MFEM  v4.6.0
Finite element discretization library
transfer.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 "transfer.hpp"
13 #include "bilinearform.hpp"
14 #include "../general/forall.hpp"
15 
16 namespace mfem
17 {
18 
20  FiniteElementSpace &ran_fes_)
21  : dom_fes(dom_fes_), ran_fes(ran_fes_),
22  oper_type(Operator::ANY_TYPE),
23  fw_t_oper(), bw_t_oper()
24 {
25 #ifdef MFEM_USE_MPI
26  const bool par_dom = dynamic_cast<ParFiniteElementSpace*>(&dom_fes);
27  const bool par_ran = dynamic_cast<ParFiniteElementSpace*>(&ran_fes);
28  MFEM_VERIFY(par_dom == par_ran, "the domain and range FE spaces must both"
29  " be either serial or parallel");
30  parallel = par_dom;
31 #endif
32 }
33 
35  FiniteElementSpace &fes_in, FiniteElementSpace &fes_out,
36  const Operator &oper, OperatorHandle &t_oper)
37 {
38  if (t_oper.Ptr())
39  {
40  return *t_oper.Ptr();
41  }
42 
43  if (!Parallel())
44  {
45  const SparseMatrix *in_cP = fes_in.GetConformingProlongation();
46  const SparseMatrix *out_cR = fes_out.GetConformingRestriction();
48  {
49  const SparseMatrix *mat = dynamic_cast<const SparseMatrix *>(&oper);
50  MFEM_VERIFY(mat != NULL, "Operator is not a SparseMatrix");
51  if (!out_cR)
52  {
53  t_oper.Reset(const_cast<SparseMatrix*>(mat), false);
54  }
55  else
56  {
57  t_oper.Reset(mfem::Mult(*out_cR, *mat));
58  }
59  if (in_cP)
60  {
61  t_oper.Reset(mfem::Mult(*t_oper.As<SparseMatrix>(), *in_cP));
62  }
63  }
64  else if (oper_type == Operator::ANY_TYPE)
65  {
66  const int RP_case = bool(out_cR) + 2*bool(in_cP);
67  switch (RP_case)
68  {
69  case 0:
70  t_oper.Reset(const_cast<Operator*>(&oper), false);
71  break;
72  case 1:
73  t_oper.Reset(
74  new ProductOperator(out_cR, &oper, false, false));
75  break;
76  case 2:
77  t_oper.Reset(
78  new ProductOperator(&oper, in_cP, false, false));
79  break;
80  case 3:
81  t_oper.Reset(
83  out_cR, &oper, in_cP, false, false, false));
84  break;
85  }
86  }
87  else
88  {
89  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
90  }
91  }
92  else // Parallel() == true
93  {
94 #ifdef MFEM_USE_MPI
96  {
97  const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
98  const ParFiniteElementSpace *pfes_in =
99  dynamic_cast<const ParFiniteElementSpace *>(&fes_in);
100  const ParFiniteElementSpace *pfes_out =
101  dynamic_cast<const ParFiniteElementSpace *>(&fes_out);
102  const SparseMatrix *sp_mat = dynamic_cast<const SparseMatrix *>(&oper);
103  const HypreParMatrix *hy_mat;
104  if (sp_mat)
105  {
106  SparseMatrix *RA = mfem::Mult(*out_R, *sp_mat);
107  t_oper.Reset(pfes_in->Dof_TrueDof_Matrix()->
108  LeftDiagMult(*RA, pfes_out->GetTrueDofOffsets()));
109  delete RA;
110  }
111  else if ((hy_mat = dynamic_cast<const HypreParMatrix *>(&oper)))
112  {
113  HypreParMatrix *RA =
114  hy_mat->LeftDiagMult(*out_R, pfes_out->GetTrueDofOffsets());
115  t_oper.Reset(mfem::ParMult(RA, pfes_in->Dof_TrueDof_Matrix()));
116  delete RA;
117  }
118  else
119  {
120  MFEM_ABORT("unknown Operator type");
121  }
122  }
123  else if (oper_type == Operator::ANY_TYPE)
124  {
125  const Operator *out_R = fes_out.GetRestrictionOperator();
126  t_oper.Reset(new TripleProductOperator(
127  out_R, &oper, fes_in.GetProlongationMatrix(),
128  false, false, false));
129  }
130  else
131  {
132  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
133  }
134 #endif
135  }
136 
137  return *t_oper.Ptr();
138 }
139 
140 
142 {
143  if (own_mass_integ) { delete mass_integ; }
144 }
145 
147  BilinearFormIntegrator *mass_integ_, bool own_mass_integ_)
148 {
149  if (own_mass_integ) { delete mass_integ; }
150 
151  mass_integ = mass_integ_;
152  own_mass_integ = own_mass_integ_;
153 }
154 
156 {
157  if (F.Ptr())
158  {
159  return *F.Ptr();
160  }
161 
162  // Construct F
164  {
166  }
168  {
169  Mesh::GeometryList elem_geoms(*ran_fes.GetMesh());
170 
172  for (int i = 0; i < elem_geoms.Size(); i++)
173  {
175  localP[elem_geoms[i]]);
176  }
180  }
181  else
182  {
183  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
184  }
185 
186  return *F.Ptr();
187 }
188 
190 {
191  if (B.Ptr())
192  {
193  return *B.Ptr();
194  }
195 
196  // Construct B, if not set, define a suitable mass_integ
197  if (!mass_integ && ran_fes.GetNE() > 0)
198  {
199  const FiniteElement *f_fe_0 = ran_fes.GetFE(0);
200  const int map_type = f_fe_0->GetMapType();
201  if (map_type == FiniteElement::VALUE ||
202  map_type == FiniteElement::INTEGRAL)
203  {
205  }
206  else if (map_type == FiniteElement::H_DIV ||
207  map_type == FiniteElement::H_CURL)
208  {
210  }
211  else
212  {
213  MFEM_ABORT("unknown type of FE space");
214  }
215  own_mass_integ = true;
216  }
218  {
220  &ran_fes, &dom_fes, mass_integ));
221  }
222  else
223  {
224  MFEM_ABORT("Operator::Type is not supported: " << oper_type);
225  }
226 
227  return *B.Ptr();
228 }
229 
230 
232  const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
233  : Operator(fes_lor_.GetVSize(), fes_ho_.GetVSize()),
234  fes_ho(fes_ho_),
235  fes_lor(fes_lor_)
236 { }
237 
239  int nel_ho, int nel_lor, const CoarseFineTransformations& cf_tr)
240 {
241  // Construct the mapping from HO to LOR
242  // ho2lor.GetRow(iho) will give all the LOR elements contained in iho
243  ho2lor.MakeI(nel_ho);
244  for (int ilor = 0; ilor < nel_lor; ++ilor)
245  {
246  int iho = cf_tr.embeddings[ilor].parent;
247  ho2lor.AddAColumnInRow(iho);
248  }
249  ho2lor.MakeJ();
250  for (int ilor = 0; ilor < nel_lor; ++ilor)
251  {
252  int iho = cf_tr.embeddings[ilor].parent;
253  ho2lor.AddConnection(iho, ilor);
254  }
255  ho2lor.ShiftUpI();
256 }
257 
259  Geometry::Type geom, const FiniteElement& fe_ho,
260  const FiniteElement& fe_lor, ElementTransformation* el_tr,
262  DenseMatrix& M_mixed_el) const
263 {
264  int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
265  const IntegrationRule* ir = &IntRules.Get(geom, order);
266  M_mixed_el = 0.0;
267  for (int i = 0; i < ir->GetNPoints(); i++)
268  {
269  const IntegrationPoint& ip_lor = ir->IntPoint(i);
270  IntegrationPoint ip_ho;
271  ip_tr.Transform(ip_lor, ip_ho);
272  Vector shape_lor(fe_lor.GetDof());
273  fe_lor.CalcShape(ip_lor, shape_lor);
274  Vector shape_ho(fe_ho.GetDof());
275  fe_ho.CalcShape(ip_ho, shape_ho);
276  el_tr->SetIntPoint(&ip_lor);
277  // For now we use the geometry information from the LOR space, which means
278  // we won't be mass conservative if the mesh is curved
279  double w = el_tr->Weight() * ip_lor.weight;
280  shape_lor *= w;
281  AddMultVWt(shape_lor, shape_ho, M_mixed_el);
282  }
283 }
284 
286  const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
287  : L2Projection(fes_ho_, fes_lor_)
288 {
289  Mesh *mesh_ho = fes_ho.GetMesh();
290  Mesh *mesh_lor = fes_lor.GetMesh();
291  int nel_ho = mesh_ho->GetNE();
292  int nel_lor = mesh_lor->GetNE();
293 
294  // The prolongation operation is only well-defined when the LOR space has at
295  // least as many DOFs as the high-order space.
296  const bool build_P = fes_lor.GetTrueVSize() >= fes_ho.GetTrueVSize();
297 
298  // If the local mesh is empty, skip all computations
299  if (nel_ho == 0) { return; }
300 
301  const CoarseFineTransformations &cf_tr = mesh_lor->GetRefinementTransforms();
302 
303  int nref_max = 0;
304  Array<Geometry::Type> geoms;
305  mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
306  for (int ig = 0; ig < geoms.Size(); ++ig)
307  {
308  Geometry::Type geom = geoms[ig];
309  nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
310  }
311 
312  BuildHo2Lor(nel_ho, nel_lor, cf_tr);
313 
314  offsets.SetSize(nel_ho+1);
315  offsets[0] = 0;
316  for (int iho = 0; iho < nel_ho; ++iho)
317  {
318  int nref = ho2lor.RowSize(iho);
319  const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
320  const FiniteElement &fe_lor = *fes_lor.GetFE(ho2lor.GetRow(iho)[0]);
321  offsets[iho+1] = offsets[iho] + fe_ho.GetDof()*fe_lor.GetDof()*nref;
322  }
323  // R will contain the restriction (L^2 projection operator) defined on each
324  // coarse HO element (and corresponding patch of LOR elements)
325  R.SetSize(offsets[nel_ho]);
326  if (build_P)
327  {
328  // P will contain the corresponding prolongation operator
329  P.SetSize(offsets[nel_ho]);
330  }
331 
333  IsoparametricTransformation &emb_tr = ip_tr.Transf;
334 
335  for (int iho = 0; iho < nel_ho; ++iho)
336  {
337  Array<int> lor_els;
338  ho2lor.GetRow(iho, lor_els);
339  int nref = ho2lor.RowSize(iho);
340 
341  Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
342  const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
343  const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
344  int ndof_ho = fe_ho.GetDof();
345  int ndof_lor = fe_lor.GetDof();
346 
347  emb_tr.SetIdentityTransformation(geom);
348  const DenseTensor &pmats = cf_tr.point_matrices[geom];
349 
350  DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
351 
352  DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
353  DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
354 
355  MassIntegrator mi;
356  DenseMatrix M_lor_el(ndof_lor, ndof_lor);
357  DenseMatrixInverse Minv_lor_el(&M_lor_el);
358  DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
359  DenseMatrix M_mixed_el(ndof_lor, ndof_ho);
360 
361  Minv_lor = 0.0;
362  M_lor = 0.0;
363 
364  DenseMatrix RtMlor(ndof_ho, ndof_lor*nref);
365  DenseMatrix RtMlorR(ndof_ho, ndof_ho);
366  DenseMatrixInverse RtMlorR_inv(&RtMlorR);
367 
368  for (int iref = 0; iref < nref; ++iref)
369  {
370  // Assemble the low-order refined mass matrix and invert locally
371  int ilor = lor_els[iref];
373  mi.AssembleElementMatrix(fe_lor, *el_tr, M_lor_el);
374  M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
375  Minv_lor_el.Factor();
376  Minv_lor_el.GetInverseMatrix(M_lor_el);
377  // Insert into the diagonal of the patch LOR mass matrix
378  Minv_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
379 
380  // Now assemble the block-row of the mixed mass matrix associated
381  // with integrating HO functions against LOR functions on the LOR
382  // sub-element.
383 
384  // Create the transformation that embeds the fine low-order element
385  // within the coarse high-order element in reference space
386  emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
387 
388  ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, M_mixed_el);
389 
390  M_mixed.CopyMN(M_mixed_el, iref*ndof_lor, 0);
391  }
392  mfem::Mult(Minv_lor, M_mixed, R_iho);
393 
394  if (build_P)
395  {
396  DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
397 
398  mfem::MultAtB(R_iho, M_lor, RtMlor);
399  mfem::Mult(RtMlor, R_iho, RtMlorR);
400  RtMlorR_inv.Factor();
401  RtMlorR_inv.Mult(RtMlor, P_iho);
402  }
403  }
404 }
405 
407  const Vector &x, Vector &y) const
408 {
409  int vdim = fes_ho.GetVDim();
410  Array<int> vdofs;
411  DenseMatrix xel_mat, yel_mat;
412  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
413  {
414  int nref = ho2lor.RowSize(iho);
415  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
416  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
417  xel_mat.SetSize(ndof_ho, vdim);
418  yel_mat.SetSize(ndof_lor*nref, vdim);
419  DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
420 
421  fes_ho.GetElementVDofs(iho, vdofs);
422  x.GetSubVector(vdofs, xel_mat.GetData());
423  mfem::Mult(R_iho, xel_mat, yel_mat);
424  // Place result correctly into the low-order vector
425  for (int iref = 0; iref < nref; ++iref)
426  {
427  int ilor = ho2lor.GetRow(iho)[iref];
428  for (int vd=0; vd<vdim; ++vd)
429  {
430  fes_lor.GetElementDofs(ilor, vdofs);
431  fes_lor.DofsToVDofs(vd, vdofs);
432  y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
433  }
434  }
435  }
436 }
437 
439  const Vector &x, Vector &y) const
440 {
441  int vdim = fes_ho.GetVDim();
442  Array<int> vdofs;
443  DenseMatrix xel_mat, yel_mat;
444  y = 0.0;
445  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
446  {
447  int nref = ho2lor.RowSize(iho);
448  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
449  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
450  xel_mat.SetSize(ndof_lor*nref, vdim);
451  yel_mat.SetSize(ndof_ho, vdim);
452  DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
453 
454  // Extract the LOR DOFs
455  for (int iref=0; iref<nref; ++iref)
456  {
457  int ilor = ho2lor.GetRow(iho)[iref];
458  for (int vd=0; vd<vdim; ++vd)
459  {
460  fes_lor.GetElementDofs(ilor, vdofs);
461  fes_lor.DofsToVDofs(vd, vdofs);
462  x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
463  }
464  }
465  // Multiply locally by the transpose
466  mfem::MultAtB(R_iho, xel_mat, yel_mat);
467  // Place the result in the HO vector
468  fes_ho.GetElementVDofs(iho, vdofs);
469  y.AddElementVector(vdofs, yel_mat.GetData());
470  }
471 }
472 
474  const Vector &x, Vector &y) const
475 {
476  if (fes_ho.GetNE() == 0) { return; }
477  MFEM_VERIFY(P.Size() > 0, "Prolongation not supported for these spaces.")
478  int vdim = fes_ho.GetVDim();
479  Array<int> vdofs;
480  DenseMatrix xel_mat,yel_mat;
481  y = 0.0;
482  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
483  {
484  int nref = ho2lor.RowSize(iho);
485  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
486  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
487  xel_mat.SetSize(ndof_lor*nref, vdim);
488  yel_mat.SetSize(ndof_ho, vdim);
489  DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
490 
491  // Extract the LOR DOFs
492  for (int iref = 0; iref < nref; ++iref)
493  {
494  int ilor = ho2lor.GetRow(iho)[iref];
495  for (int vd = 0; vd < vdim; ++vd)
496  {
497  fes_lor.GetElementDofs(ilor, vdofs);
498  fes_lor.DofsToVDofs(vd, vdofs);
499  x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
500  }
501  }
502  // Locally prolongate
503  mfem::Mult(P_iho, xel_mat, yel_mat);
504  // Place the result in the HO vector
505  fes_ho.GetElementVDofs(iho, vdofs);
506  y.AddElementVector(vdofs, yel_mat.GetData());
507  }
508 }
509 
511  const Vector &x, Vector &y) const
512 {
513  if (fes_ho.GetNE() == 0) { return; }
514  MFEM_VERIFY(P.Size() > 0, "Prolongation not supported for these spaces.")
515  int vdim = fes_ho.GetVDim();
516  Array<int> vdofs;
517  DenseMatrix xel_mat,yel_mat;
518  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
519  {
520  int nref = ho2lor.RowSize(iho);
521  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
522  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
523  xel_mat.SetSize(ndof_ho, vdim);
524  yel_mat.SetSize(ndof_lor*nref, vdim);
525  DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
526 
527  fes_ho.GetElementVDofs(iho, vdofs);
528  x.GetSubVector(vdofs, xel_mat.GetData());
529  mfem::MultAtB(P_iho, xel_mat, yel_mat);
530 
531  // Place result correctly into the low-order vector
532  for (int iref = 0; iref < nref; ++iref)
533  {
534  int ilor = ho2lor.GetRow(iho)[iref];
535  for (int vd=0; vd<vdim; ++vd)
536  {
537  fes_lor.GetElementDofs(ilor, vdofs);
538  fes_lor.DofsToVDofs(vd, vdofs);
539  y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
540  }
541  }
542  }
543 }
544 
546  const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_)
547  : L2Projection(fes_ho_, fes_lor_)
548 {
549  std::unique_ptr<SparseMatrix> R_mat, M_LH_mat;
550  std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH();
551 
552  FiniteElementSpace fes_ho_scalar(fes_ho.GetMesh(), fes_ho.FEColl(), 1);
553  FiniteElementSpace fes_lor_scalar(fes_lor.GetMesh(), fes_lor.FEColl(), 1);
554 
555  const SparseMatrix *P_ho = fes_ho_scalar.GetConformingProlongation();
556  const SparseMatrix *P_lor = fes_lor_scalar.GetConformingProlongation();
557 
558  if (P_ho || P_lor)
559  {
560  if (P_ho && P_lor)
561  {
562  R_mat.reset(RAP(*P_lor, *R_mat, *P_ho));
563  M_LH_mat.reset(RAP(*P_lor, *M_LH_mat, *P_ho));
564  }
565  else if (P_ho)
566  {
567  R_mat.reset(mfem::Mult(*R_mat, *P_ho));
568  M_LH_mat.reset(mfem::Mult(*M_LH_mat, *P_ho));
569  }
570  else // P_lor != nullptr
571  {
572  R_mat.reset(mfem::Mult(*P_lor, *R_mat));
573  M_LH_mat.reset(mfem::Mult(*P_lor, *M_LH_mat));
574  }
575  }
576 
577  SparseMatrix *RTxM_LH_mat = TransposeMult(*R_mat, *M_LH_mat);
578  precon.reset(new DSmoother(*RTxM_LH_mat));
579 
580  // Set ownership
581  RTxM_LH.reset(RTxM_LH_mat);
582  R = std::move(R_mat);
583  M_LH = std::move(M_LH_mat);
584 
585  SetupPCG();
586 }
587 
588 #ifdef MFEM_USE_MPI
589 
591  const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor)
592  : L2Projection(pfes_ho, pfes_lor),
593  pcg(pfes_ho.GetComm())
594 {
595  std::tie(R, M_LH) = ComputeSparseRAndM_LH();
596 
597  ParFiniteElementSpace pfes_ho_scalar(pfes_ho.GetParMesh(),
598  pfes_ho.FEColl(), 1);
599  ParFiniteElementSpace pfes_lor_scalar(pfes_lor.GetParMesh(),
600  pfes_lor.FEColl(), 1);
601 
602  HypreParMatrix R_local = HypreParMatrix(pfes_ho.GetComm(),
603  pfes_lor_scalar.GlobalVSize(),
604  pfes_ho_scalar.GlobalVSize(),
605  pfes_lor_scalar.GetDofOffsets(),
606  pfes_ho_scalar.GetDofOffsets(),
607  static_cast<SparseMatrix*>(R.get()));
608  HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho.GetComm(),
609  pfes_lor_scalar.GlobalVSize(),
610  pfes_ho_scalar.GlobalVSize(),
611  pfes_lor_scalar.GetDofOffsets(),
612  pfes_ho_scalar.GetDofOffsets(),
613  static_cast<SparseMatrix*>(M_LH.get()));
614 
615  HypreParMatrix *R_mat = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(),
616  &R_local, pfes_ho_scalar.Dof_TrueDof_Matrix());
617  HypreParMatrix *M_LH_mat = RAP(pfes_lor_scalar.Dof_TrueDof_Matrix(),
618  &M_LH_local, pfes_ho_scalar.Dof_TrueDof_Matrix());
619 
620  std::unique_ptr<HypreParMatrix> R_T(R_mat->Transpose());
621  HypreParMatrix *RTxM_LH_mat = ParMult(R_T.get(), M_LH_mat, true);
622 
623  HypreBoomerAMG *amg = new HypreBoomerAMG(*RTxM_LH_mat);
624  amg->SetPrintLevel(0);
625 
626  R.reset(R_mat);
627  M_LH.reset(M_LH_mat);
628  RTxM_LH.reset(RTxM_LH_mat);
629  precon.reset(amg);
630 
631  SetupPCG();
634 }
635 
636 #endif
637 
639 {
640  // Basic PCG solver setup
641  pcg.SetPrintLevel(0);
642  // pcg.SetPrintLevel(IterativeSolver::PrintLevel().Summary());
643  pcg.SetMaxIter(1000);
644  // initial values for relative and absolute tolerance
645  pcg.SetRelTol(1e-13);
646  pcg.SetAbsTol(1e-13);
647  pcg.SetPreconditioner(*precon);
648  pcg.SetOperator(*RTxM_LH);
649 }
650 
652  const Vector& x, Vector& y) const
653 {
654  Vector X(fes_ho.GetTrueVSize());
655  Vector X_dim(R->Width());
656 
657  Vector Y_dim(R->Height());
658  Vector Y(fes_lor.GetTrueVSize());
659 
660  Array<int> vdofs_list;
661 
662  GetTDofs(fes_ho, x, X);
663 
664  for (int d = 0; d < fes_ho.GetVDim(); ++d)
665  {
666  TDofsListByVDim(fes_ho, d, vdofs_list);
667  X.GetSubVector(vdofs_list, X_dim);
668  R->Mult(X_dim, Y_dim);
669  TDofsListByVDim(fes_lor, d, vdofs_list);
670  Y.SetSubVector(vdofs_list, Y_dim);
671  }
672 
673  SetFromTDofs(fes_lor, Y, y);
674 }
675 
677  const Vector& x, Vector& y) const
678 {
679  Vector X(fes_lor.GetTrueVSize());
680  Vector X_dim(R->Height());
681 
682  Vector Y_dim(R->Width());
683  Vector Y(fes_ho.GetTrueVSize());
684 
685  Array<int> vdofs_list;
686 
687  GetTDofsTranspose(fes_lor, x, X);
688 
689  for (int d = 0; d < fes_ho.GetVDim(); ++d)
690  {
691  TDofsListByVDim(fes_lor, d, vdofs_list);
692  X.GetSubVector(vdofs_list, X_dim);
693  R->MultTranspose(X_dim, Y_dim);
694  TDofsListByVDim(fes_ho, d, vdofs_list);
695  Y.SetSubVector(vdofs_list, Y_dim);
696  }
697 
698  SetFromTDofsTranspose(fes_ho, Y, y);
699 }
700 
702  const Vector& x, Vector& y) const
703 {
704  Vector X(fes_lor.GetTrueVSize());
705  Vector X_dim(M_LH->Height());
706  Vector Xbar(pcg.Width());
707 
708  Vector Y_dim(pcg.Height());
709  Vector Y(fes_ho.GetTrueVSize());
710 
711  Array<int> vdofs_list;
712 
713  GetTDofs(fes_lor, x, X);
714 
715  for (int d = 0; d < fes_ho.GetVDim(); ++d)
716  {
717  TDofsListByVDim(fes_lor, d, vdofs_list);
718  X.GetSubVector(vdofs_list, X_dim);
719  // Compute y = P x = (R^T M_LH)^(-1) M_LH^T X = (R^T M_LH)^(-1) Xbar
720  M_LH->MultTranspose(X_dim, Xbar);
721  Y_dim = 0.0;
722  pcg.Mult(Xbar, Y_dim);
723  TDofsListByVDim(fes_ho, d, vdofs_list);
724  Y.SetSubVector(vdofs_list, Y_dim);
725  }
726 
727  SetFromTDofs(fes_ho, Y, y);
728 }
729 
731  const Vector& x, Vector& y) const
732 {
733  Vector X(fes_ho.GetTrueVSize());
734  Vector X_dim(pcg.Width());
735  Vector Xbar(pcg.Height());
736 
737  Vector Y_dim(M_LH->Height());
738  Vector Y(fes_lor.GetTrueVSize());
739 
740  Array<int> vdofs_list;
741 
742  GetTDofsTranspose(fes_ho, x, X);
743 
744  for (int d = 0; d < fes_ho.GetVDim(); ++d)
745  {
746  TDofsListByVDim(fes_ho, d, vdofs_list);
747  X.GetSubVector(vdofs_list, X_dim);
748  // Compute y = P^T x = M_LH (R^T M_LH)^(-1) X = M_LH Xbar
749  Xbar = 0.0;
750  pcg.Mult(X_dim, Xbar);
751  M_LH->Mult(Xbar, Y_dim);
752  TDofsListByVDim(fes_lor, d, vdofs_list);
753  Y.SetSubVector(vdofs_list, Y_dim);
754  }
755 
756  SetFromTDofsTranspose(fes_lor, Y, y);
757 }
758 
760 {
761  pcg.SetRelTol(p_rtol_);
762 }
763 
765 {
766  pcg.SetAbsTol(p_atol_);
767 }
768 
769 std::pair<
770 std::unique_ptr<SparseMatrix>,
771 std::unique_ptr<SparseMatrix>>
773 {
774  std::pair<std::unique_ptr<SparseMatrix>,
775  std::unique_ptr<SparseMatrix>> r_and_mlh;
776 
777  Mesh* mesh_ho = fes_ho.GetMesh();
778  Mesh* mesh_lor = fes_lor.GetMesh();
779  int nel_ho = mesh_ho->GetNE();
780  int nel_lor = mesh_lor->GetNE();
781  int ndof_lor = fes_lor.GetNDofs();
782 
783  // If the local mesh is empty, skip all computations
784  if (nel_ho == 0) { return {nullptr, nullptr}; }
785 
786  const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
787 
788  int nref_max = 0;
789  Array<Geometry::Type> geoms;
790  mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
791  for (int ig = 0; ig < geoms.Size(); ++ig)
792  {
793  Geometry::Type geom = geoms[ig];
794  nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
795  }
796 
797  BuildHo2Lor(nel_ho, nel_lor, cf_tr);
798 
799  // ML_inv contains the inverse lumped (row sum) mass matrix. Note that the
800  // method will also work with a full (consistent) mass matrix, though this is
801  // not implemented here. L refers to the low-order refined mesh
802  Vector ML_inv(ndof_lor);
803  ML_inv = 0.0;
804 
805  // Compute ML_inv
806  for (int iho = 0; iho < nel_ho; ++iho)
807  {
808  Array<int> lor_els;
809  ho2lor.GetRow(iho, lor_els);
810  int nref = ho2lor.RowSize(iho);
811 
812  Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
813  const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
814  int nedof_lor = fe_lor.GetDof();
815 
816  // Instead of using a MassIntegrator, manually loop over integration
817  // points so we can row sum and store the diagonal as a Vector.
818  Vector ML_el(nedof_lor);
819  Vector shape_lor(nedof_lor);
820  Array<int> dofs_lor(nedof_lor);
821 
822  for (int iref = 0; iref < nref; ++iref)
823  {
824  int ilor = lor_els[iref];
825  ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
826 
827  int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
828  const IntegrationRule* ir = &IntRules.Get(geom, order);
829  ML_el = 0.0;
830  for (int i = 0; i < ir->GetNPoints(); ++i)
831  {
832  const IntegrationPoint& ip_lor = ir->IntPoint(i);
833  fe_lor.CalcShape(ip_lor, shape_lor);
834  el_tr->SetIntPoint(&ip_lor);
835  ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
836  }
837  fes_lor.GetElementDofs(ilor, dofs_lor);
838  ML_inv.AddElementVector(dofs_lor, ML_el);
839  }
840  }
841  // DOF by DOF inverse of non-zero entries
842  LumpedMassInverse(ML_inv);
843 
844  // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate
845  r_and_mlh.first = AllocR();
846  // Allocate M_LH (same sparsity pattern as R)
847  // L refers to the low-order refined mesh (DOFs correspond to rows)
848  // H refers to the higher-order mesh (DOFs correspond to columns)
849  Memory<int> I(r_and_mlh.first->Height() + 1);
850  for (int icol = 0; icol < r_and_mlh.first->Height() + 1; ++icol)
851  {
852  I[icol] = r_and_mlh.first->GetI()[icol];
853  }
854  Memory<int> J(r_and_mlh.first->NumNonZeroElems());
855  for (int jcol = 0; jcol < r_and_mlh.first->NumNonZeroElems(); ++jcol)
856  {
857  J[jcol] = r_and_mlh.first->GetJ()[jcol];
858  }
859  r_and_mlh.second = std::unique_ptr<SparseMatrix>(
860  new SparseMatrix(I, J, NULL, r_and_mlh.first->Height(),
861  r_and_mlh.first->Width(), true, true, true));
862 
864  IsoparametricTransformation& emb_tr = ip_tr.Transf;
865 
866  // Compute M_LH and R
867  for (int iho = 0; iho < nel_ho; ++iho)
868  {
869  Array<int> lor_els;
870  ho2lor.GetRow(iho, lor_els);
871  int nref = ho2lor.RowSize(iho);
872 
873  Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
874  const FiniteElement& fe_ho = *fes_ho.GetFE(iho);
875  const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
876 
877  emb_tr.SetIdentityTransformation(geom);
878  const DenseTensor& pmats = cf_tr.point_matrices[geom];
879 
880  int nedof_ho = fe_ho.GetDof();
881  int nedof_lor = fe_lor.GetDof();
882  DenseMatrix M_LH_el(nedof_lor, nedof_ho);
883  DenseMatrix R_el(nedof_lor, nedof_ho);
884 
885  for (int iref = 0; iref < nref; ++iref)
886  {
887  int ilor = lor_els[iref];
888  ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
889 
890  // Create the transformation that embeds the fine low-order element
891  // within the coarse high-order element in reference space
892  emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
893 
894  ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, M_LH_el);
895 
896  Array<int> dofs_lor(nedof_lor);
897  fes_lor.GetElementDofs(ilor, dofs_lor);
898  Vector R_row;
899  for (int i = 0; i < nedof_lor; ++i)
900  {
901  M_LH_el.GetRow(i, R_row);
902  R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
903  }
904  Array<int> dofs_ho(nedof_ho);
905  fes_ho.GetElementDofs(iho, dofs_ho);
906  r_and_mlh.second->AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
907  r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
908  }
909  }
910 
911  return r_and_mlh;
912 }
913 
915  const FiniteElementSpace& fes, const Vector& x, Vector& X) const
916 {
917  const Operator* res = fes.GetRestrictionOperator();
918  if (res)
919  {
920  res->Mult(x, X);
921  }
922  else
923  {
924  X = x;
925  }
926 }
927 
929  const FiniteElementSpace& fes, const Vector &X, Vector& x) const
930 {
931  const Operator* P = fes.GetProlongationMatrix();
932  if (P)
933  {
934  P->Mult(X, x);
935  }
936  else
937  {
938  x = X;
939  }
940 }
941 
943  const FiniteElementSpace& fes, const Vector& x, Vector& X) const
944 {
945  const Operator* P = fes.GetProlongationMatrix();
946  if (P)
947  {
948  P->MultTranspose(x, X);
949  }
950  else
951  {
952  X = x;
953  }
954 }
955 
957  const FiniteElementSpace& fes, const Vector &X, Vector& x) const
958 {
959  const Operator *R_op = fes.GetRestrictionOperator();
960  if (R_op)
961  {
962  R_op->MultTranspose(X, x);
963  }
964  else
965  {
966  x = X;
967  }
968 }
969 
971  const FiniteElementSpace& fes, int vdim, Array<int>& vdofs_list) const
972 {
973  const SparseMatrix *R_mat = fes.GetRestrictionMatrix();
974  if (R_mat)
975  {
976  Array<int> x_vdofs_list(fes.GetNDofs());
977  Array<int> x_vdofs_marker(fes.GetVSize());
978  Array<int> X_vdofs_marker(fes.GetTrueVSize());
979  fes.GetVDofs(vdim, x_vdofs_list);
980  FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker);
981  R_mat->BooleanMult(x_vdofs_marker, X_vdofs_marker);
982  FiniteElementSpace::MarkerToList(X_vdofs_marker, vdofs_list);
983  }
984  else
985  {
986  vdofs_list.SetSize(fes.GetNDofs());
987  fes.GetVDofs(vdim, vdofs_list);
988  }
989 }
990 
992  Vector& ML_inv) const
993 {
994  Vector ML_inv_full(fes_lor.GetVSize());
995  // set ML_inv on dofs for vdim = 0
996  Array<int> vdofs_list(fes_lor.GetNDofs());
997  fes_lor.GetVDofs(0, vdofs_list);
998  ML_inv_full.SetSubVector(vdofs_list, ML_inv);
999 
1000  Vector ML_inv_true(fes_lor.GetTrueVSize());
1001  const Operator *P = fes_lor.GetProlongationMatrix();
1002  if (P) { P->MultTranspose(ML_inv_full, ML_inv_true); }
1003  else { ML_inv_true = ML_inv_full; }
1004 
1005  for (int i = 0; i < ML_inv_true.Size(); ++i)
1006  {
1007  ML_inv_true[i] = 1.0 / ML_inv_true[i];
1008  }
1009 
1010  if (P) { P->Mult(ML_inv_true, ML_inv_full); }
1011  else { ML_inv_full = ML_inv_true; }
1012 
1013  ML_inv_full.GetSubVector(vdofs_list, ML_inv);
1014 }
1015 
1016 std::unique_ptr<SparseMatrix>
1018 {
1019  const Table& elem_dof_ho = fes_ho.GetElementToDofTable();
1020  const Table& elem_dof_lor = fes_lor.GetElementToDofTable();
1021  const int ndof_ho = fes_ho.GetNDofs();
1022  const int ndof_lor = fes_lor.GetNDofs();
1023 
1024  Table dof_elem_lor;
1025  Transpose(elem_dof_lor, dof_elem_lor, ndof_lor);
1026 
1027  Mesh* mesh_lor = fes_lor.GetMesh();
1028  const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
1029 
1030  // mfem::Mult but uses ho2lor to map HO elements to LOR elements
1031  const int* elem_dof_hoI = elem_dof_ho.GetI();
1032  const int* elem_dof_hoJ = elem_dof_ho.GetJ();
1033  const int* dof_elem_lorI = dof_elem_lor.GetI();
1034  const int* dof_elem_lorJ = dof_elem_lor.GetJ();
1035 
1036  Array<int> I(ndof_lor + 1);
1037 
1038  // figure out the size of J
1039  Array<int> dof_used_ho;
1040  dof_used_ho.SetSize(ndof_ho, -1);
1041 
1042  int sizeJ = 0;
1043  for (int ilor = 0; ilor < ndof_lor; ++ilor)
1044  {
1045  for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
1046  {
1047  int el_lor = dof_elem_lorJ[jlor];
1048  int iho = cf_tr.embeddings[el_lor].parent;
1049  for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
1050  {
1051  int dof_ho = elem_dof_hoJ[jho];
1052  if (dof_used_ho[dof_ho] != ilor)
1053  {
1054  dof_used_ho[dof_ho] = ilor;
1055  ++sizeJ;
1056  }
1057  }
1058  }
1059  }
1060 
1061  // initialize dof_ho_dof_lor
1062  Table dof_lor_dof_ho;
1063  dof_lor_dof_ho.SetDims(ndof_lor, sizeJ);
1064 
1065  for (int i = 0; i < ndof_ho; ++i)
1066  {
1067  dof_used_ho[i] = -1;
1068  }
1069 
1070  // set values of J
1071  int* dof_dofI = dof_lor_dof_ho.GetI();
1072  int* dof_dofJ = dof_lor_dof_ho.GetJ();
1073  sizeJ = 0;
1074  for (int ilor = 0; ilor < ndof_lor; ++ilor)
1075  {
1076  dof_dofI[ilor] = sizeJ;
1077  for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
1078  {
1079  int el_lor = dof_elem_lorJ[jlor];
1080  int iho = cf_tr.embeddings[el_lor].parent;
1081  for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
1082  {
1083  int dof_ho = elem_dof_hoJ[jho];
1084  if (dof_used_ho[dof_ho] != ilor)
1085  {
1086  dof_used_ho[dof_ho] = ilor;
1087  dof_dofJ[sizeJ] = dof_ho;
1088  ++sizeJ;
1089  }
1090  }
1091  }
1092  }
1093 
1094  dof_lor_dof_ho.SortRows();
1095  double* data = Memory<double>(dof_dofI[ndof_lor]);
1096 
1097  std::unique_ptr<SparseMatrix> R_local(new SparseMatrix(
1098  dof_dofI, dof_dofJ, data, ndof_lor,
1099  ndof_ho, true, true, true));
1100  (*R_local) = 0.0;
1101 
1102  dof_lor_dof_ho.LoseData();
1103 
1104  return R_local;
1105 }
1106 
1108 {
1109  delete F;
1110  delete B;
1111 }
1112 
1114 {
1115  if (!F) { BuildF(); }
1116  return *F;
1117 }
1118 
1120 {
1121  if (!B)
1122  {
1123  if (!F) { BuildF(); }
1124  B = new L2Prolongation(*F);
1125  }
1126  return *B;
1127 }
1128 
1129 void L2ProjectionGridTransfer::BuildF()
1130 {
1131  if (!force_l2_space &&
1133  {
1134  if (!Parallel())
1135  {
1136  F = new L2ProjectionH1Space(dom_fes, ran_fes);
1137  }
1138  else
1139  {
1140 #ifdef MFEM_USE_MPI
1141  const mfem::ParFiniteElementSpace& dom_pfes =
1142  static_cast<mfem::ParFiniteElementSpace&>(dom_fes);
1143  const mfem::ParFiniteElementSpace& ran_pfes =
1144  static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
1145  F = new L2ProjectionH1Space(dom_pfes, ran_pfes);
1146 #endif
1147  }
1148  }
1149  else
1150  {
1151  F = new L2ProjectionL2Space(dom_fes, ran_fes);
1152  }
1153 }
1154 
1156 {
1157  return ran_fes.GetTrueVSize() >= dom_fes.GetTrueVSize();
1158 }
1159 
1160 
1162  const FiniteElementSpace& hFESpace_)
1163  : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize())
1164 {
1165  bool isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
1166  if (lFESpace_.FEColl() == hFESpace_.FEColl() && !isvar_order)
1167  {
1169  hFESpace_.GetTransferOperator(lFESpace_, P);
1170  P.SetOperatorOwner(false);
1171  opr = P.Ptr();
1172  }
1173  else if (lFESpace_.GetMesh()->GetNE() > 0
1174  && hFESpace_.GetMesh()->GetNE() > 0
1175  && lFESpace_.GetVDim() == 1
1176  && hFESpace_.GetVDim() == 1
1177  && dynamic_cast<const TensorBasisElement*>(lFESpace_.GetFE(0))
1178  && dynamic_cast<const TensorBasisElement*>(hFESpace_.GetFE(0))
1179  && !isvar_order
1180  && (hFESpace_.FEColl()->GetContType() ==
1182  hFESpace_.FEColl()->GetContType() ==
1184  {
1185  opr = new TensorProductPRefinementTransferOperator(lFESpace_, hFESpace_);
1186  }
1187  else
1188  {
1189  opr = new PRefinementTransferOperator(lFESpace_, hFESpace_);
1190  }
1191 }
1192 
1194 
1195 void TransferOperator::Mult(const Vector& x, Vector& y) const
1196 {
1197  opr->Mult(x, y);
1198 }
1199 
1201 {
1202  opr->MultTranspose(x, y);
1203 }
1204 
1205 
1207  const FiniteElementSpace& lFESpace_, const FiniteElementSpace& hFESpace_)
1208  : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
1209  hFESpace(hFESpace_)
1210 {
1211  isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
1212 }
1213 
1215 
1217 {
1218  Mesh* mesh = hFESpace.GetMesh();
1219  Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
1220  DenseMatrix loc_prol;
1221  Vector subY, subX;
1222 
1223  Geometry::Type cached_geom = Geometry::INVALID;
1224  const FiniteElement* h_fe = NULL;
1225  const FiniteElement* l_fe = NULL;
1227 
1228  int vdim = lFESpace.GetVDim();
1229 
1230  for (int i = 0; i < mesh->GetNE(); i++)
1231  {
1232  DofTransformation * doftrans_h = hFESpace.GetElementDofs(i, h_dofs);
1233  DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
1234 
1235  const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
1236  if (geom != cached_geom || isvar_order)
1237  {
1238  h_fe = hFESpace.GetFE(i);
1239  l_fe = lFESpace.GetFE(i);
1240  T.SetIdentityTransformation(h_fe->GetGeomType());
1241  h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
1242  subY.SetSize(loc_prol.Height());
1243  cached_geom = geom;
1244  }
1245 
1246  for (int vd = 0; vd < vdim; vd++)
1247  {
1248  l_dofs.Copy(l_vdofs);
1249  lFESpace.DofsToVDofs(vd, l_vdofs);
1250  h_dofs.Copy(h_vdofs);
1251  hFESpace.DofsToVDofs(vd, h_vdofs);
1252  x.GetSubVector(l_vdofs, subX);
1253  if (doftrans_l)
1254  {
1255  doftrans_l->InvTransformPrimal(subX);
1256  }
1257  loc_prol.Mult(subX, subY);
1258  if (doftrans_h)
1259  {
1260  doftrans_h->TransformPrimal(subY);
1261  }
1262  y.SetSubVector(h_vdofs, subY);
1263  }
1264  }
1265 }
1266 
1268  Vector& y) const
1269 {
1270  y = 0.0;
1271 
1272  Mesh* mesh = hFESpace.GetMesh();
1273  Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
1274  DenseMatrix loc_prol;
1275  Vector subY, subX;
1276 
1277  Array<char> processed(hFESpace.GetVSize());
1278  processed = 0;
1279 
1280  Geometry::Type cached_geom = Geometry::INVALID;
1281  const FiniteElement* h_fe = NULL;
1282  const FiniteElement* l_fe = NULL;
1284 
1285  int vdim = lFESpace.GetVDim();
1286 
1287  for (int i = 0; i < mesh->GetNE(); i++)
1288  {
1289  DofTransformation * doftrans_h = hFESpace.GetElementDofs(i, h_dofs);
1290  DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
1291 
1292  const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
1293  if (geom != cached_geom || isvar_order)
1294  {
1295  h_fe = hFESpace.GetFE(i);
1296  l_fe = lFESpace.GetFE(i);
1297  T.SetIdentityTransformation(h_fe->GetGeomType());
1298  h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
1299  loc_prol.Transpose();
1300  subY.SetSize(loc_prol.Height());
1301  cached_geom = geom;
1302  }
1303 
1304  for (int vd = 0; vd < vdim; vd++)
1305  {
1306  l_dofs.Copy(l_vdofs);
1307  lFESpace.DofsToVDofs(vd, l_vdofs);
1308  h_dofs.Copy(h_vdofs);
1309  hFESpace.DofsToVDofs(vd, h_vdofs);
1310 
1311  x.GetSubVector(h_vdofs, subX);
1312  if (doftrans_h)
1313  {
1314  doftrans_h->InvTransformDual(subX);
1315  }
1316  for (int p = 0; p < h_dofs.Size(); ++p)
1317  {
1318  if (processed[lFESpace.DecodeDof(h_dofs[p])])
1319  {
1320  subX[p] = 0.0;
1321  }
1322  }
1323 
1324  loc_prol.Mult(subX, subY);
1325  if (doftrans_l)
1326  {
1327  doftrans_l->TransformDual(subY);
1328  }
1329  y.AddElementVector(l_vdofs, subY);
1330  }
1331 
1332  for (int p = 0; p < h_dofs.Size(); ++p)
1333  {
1334  processed[lFESpace.DecodeDof(h_dofs[p])] = 1;
1335  }
1336  }
1337 }
1338 
1339 
1342  const FiniteElementSpace& lFESpace_,
1343  const FiniteElementSpace& hFESpace_)
1344  : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
1345  hFESpace(hFESpace_)
1346 {
1347  // Assuming the same element type
1348  Mesh* mesh = lFESpace.GetMesh();
1349  dim = mesh->Dimension();
1350  if (mesh->GetNE() == 0)
1351  {
1352  return;
1353  }
1354  const FiniteElement& el = *lFESpace.GetFE(0);
1355 
1356  const TensorBasisElement* ltel =
1357  dynamic_cast<const TensorBasisElement*>(&el);
1358  MFEM_VERIFY(ltel, "Low order FE space must be tensor product space");
1359 
1360  const TensorBasisElement* htel =
1361  dynamic_cast<const TensorBasisElement*>(hFESpace.GetFE(0));
1362  MFEM_VERIFY(htel, "High order FE space must be tensor product space");
1363  const Array<int>& hdofmap = htel->GetDofMap();
1364 
1365  const IntegrationRule& ir = hFESpace.GetFE(0)->GetNodes();
1366  IntegrationRule irLex = ir;
1367 
1368  // The quadrature points, or equivalently, the dofs of the high order space
1369  // must be sorted in lexicographical order
1370  for (int i = 0; i < ir.GetNPoints(); ++i)
1371  {
1372  int j = hdofmap[i] >=0 ? hdofmap[i] : -1 - hdofmap[i];
1373  irLex.IntPoint(i) = ir.IntPoint(j);
1374  }
1375 
1376  NE = lFESpace.GetNE();
1377  const DofToQuad& maps = el.GetDofToQuad(irLex, DofToQuad::TENSOR);
1378 
1379  D1D = maps.ndof;
1380  Q1D = maps.nqpt;
1381  B = maps.B;
1382  Bt = maps.Bt;
1383 
1384  elem_restrict_lex_l =
1386 
1387  MFEM_VERIFY(elem_restrict_lex_l,
1388  "Low order ElementRestriction not available");
1389 
1390  elem_restrict_lex_h =
1392 
1393  MFEM_VERIFY(elem_restrict_lex_h,
1394  "High order ElementRestriction not available");
1395 
1396  localL.SetSize(elem_restrict_lex_l->Height(), Device::GetMemoryType());
1397  localH.SetSize(elem_restrict_lex_h->Height(), Device::GetMemoryType());
1398  localL.UseDevice(true);
1399  localH.UseDevice(true);
1400 
1401  MFEM_VERIFY(dynamic_cast<const ElementRestriction*>(elem_restrict_lex_h),
1402  "High order element restriction is of unsupported type");
1403 
1404  mask.SetSize(localH.Size(), Device::GetMemoryType());
1405  static_cast<const ElementRestriction*>(elem_restrict_lex_h)
1406  ->BooleanMask(mask);
1407  mask.UseDevice(true);
1408 }
1409 
1410 namespace TransferKernels
1411 {
1412 void Prolongation2D(const int NE, const int D1D, const int Q1D,
1413  const Vector& localL, Vector& localH,
1414  const Array<double>& B, const Vector& mask)
1415 {
1416  auto x_ = Reshape(localL.Read(), D1D, D1D, NE);
1417  auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, NE);
1418  auto B_ = Reshape(B.Read(), Q1D, D1D);
1419  auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
1420 
1421  localH = 0.0;
1422 
1423  mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1424  {
1425  for (int dy = 0; dy < D1D; ++dy)
1426  {
1427  double sol_x[DofQuadLimits::MAX_Q1D];
1428  for (int qy = 0; qy < Q1D; ++qy)
1429  {
1430  sol_x[qy] = 0.0;
1431  }
1432  for (int dx = 0; dx < D1D; ++dx)
1433  {
1434  const double s = x_(dx, dy, e);
1435  for (int qx = 0; qx < Q1D; ++qx)
1436  {
1437  sol_x[qx] += B_(qx, dx) * s;
1438  }
1439  }
1440  for (int qy = 0; qy < Q1D; ++qy)
1441  {
1442  const double d2q = B_(qy, dy);
1443  for (int qx = 0; qx < Q1D; ++qx)
1444  {
1445  y_(qx, qy, e) += d2q * sol_x[qx];
1446  }
1447  }
1448  }
1449  for (int qy = 0; qy < Q1D; ++qy)
1450  {
1451  for (int qx = 0; qx < Q1D; ++qx)
1452  {
1453  y_(qx, qy, e) *= m_(qx, qy, e);
1454  }
1455  }
1456  });
1457 }
1458 
1459 void Prolongation3D(const int NE, const int D1D, const int Q1D,
1460  const Vector& localL, Vector& localH,
1461  const Array<double>& B, const Vector& mask)
1462 {
1463  auto x_ = Reshape(localL.Read(), D1D, D1D, D1D, NE);
1464  auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, Q1D, NE);
1465  auto B_ = Reshape(B.Read(), Q1D, D1D);
1466  auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
1467 
1468  localH = 0.0;
1469 
1470  mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1471  {
1472  for (int dz = 0; dz < D1D; ++dz)
1473  {
1474  double sol_xy[DofQuadLimits::MAX_Q1D][DofQuadLimits::MAX_Q1D];
1475  for (int qy = 0; qy < Q1D; ++qy)
1476  {
1477  for (int qx = 0; qx < Q1D; ++qx)
1478  {
1479  sol_xy[qy][qx] = 0.0;
1480  }
1481  }
1482  for (int dy = 0; dy < D1D; ++dy)
1483  {
1484  double sol_x[DofQuadLimits::MAX_Q1D];
1485  for (int qx = 0; qx < Q1D; ++qx)
1486  {
1487  sol_x[qx] = 0;
1488  }
1489  for (int dx = 0; dx < D1D; ++dx)
1490  {
1491  const double s = x_(dx, dy, dz, e);
1492  for (int qx = 0; qx < Q1D; ++qx)
1493  {
1494  sol_x[qx] += B_(qx, dx) * s;
1495  }
1496  }
1497  for (int qy = 0; qy < Q1D; ++qy)
1498  {
1499  const double wy = B_(qy, dy);
1500  for (int qx = 0; qx < Q1D; ++qx)
1501  {
1502  sol_xy[qy][qx] += wy * sol_x[qx];
1503  }
1504  }
1505  }
1506  for (int qz = 0; qz < Q1D; ++qz)
1507  {
1508  const double wz = B_(qz, dz);
1509  for (int qy = 0; qy < Q1D; ++qy)
1510  {
1511  for (int qx = 0; qx < Q1D; ++qx)
1512  {
1513  y_(qx, qy, qz, e) += wz * sol_xy[qy][qx];
1514  }
1515  }
1516  }
1517  }
1518  for (int qz = 0; qz < Q1D; ++qz)
1519  {
1520  for (int qy = 0; qy < Q1D; ++qy)
1521  {
1522  for (int qx = 0; qx < Q1D; ++qx)
1523  {
1524  y_(qx, qy, qz, e) *= m_(qx, qy, qz, e);
1525  }
1526  }
1527  }
1528  });
1529 }
1530 
1531 void Restriction2D(const int NE, const int D1D, const int Q1D,
1532  const Vector& localH, Vector& localL,
1533  const Array<double>& Bt, const Vector& mask)
1534 {
1535  auto x_ = Reshape(localH.Read(), Q1D, Q1D, NE);
1536  auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, NE);
1537  auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
1538  auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
1539 
1540  localL = 0.0;
1541 
1542  mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1543  {
1544  for (int qy = 0; qy < Q1D; ++qy)
1545  {
1546  double sol_x[DofQuadLimits::MAX_D1D];
1547  for (int dx = 0; dx < D1D; ++dx)
1548  {
1549  sol_x[dx] = 0.0;
1550  }
1551  for (int qx = 0; qx < Q1D; ++qx)
1552  {
1553  const double s = m_(qx, qy, e) * x_(qx, qy, e);
1554  for (int dx = 0; dx < D1D; ++dx)
1555  {
1556  sol_x[dx] += Bt_(dx, qx) * s;
1557  }
1558  }
1559  for (int dy = 0; dy < D1D; ++dy)
1560  {
1561  const double q2d = Bt_(dy, qy);
1562  for (int dx = 0; dx < D1D; ++dx)
1563  {
1564  y_(dx, dy, e) += q2d * sol_x[dx];
1565  }
1566  }
1567  }
1568  });
1569 }
1570 void Restriction3D(const int NE, const int D1D, const int Q1D,
1571  const Vector& localH, Vector& localL,
1572  const Array<double>& Bt, const Vector& mask)
1573 {
1574  auto x_ = Reshape(localH.Read(), Q1D, Q1D, Q1D, NE);
1575  auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, D1D, NE);
1576  auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
1577  auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
1578 
1579  localL = 0.0;
1580 
1581  mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
1582  {
1583  for (int qz = 0; qz < Q1D; ++qz)
1584  {
1585  double sol_xy[DofQuadLimits::MAX_D1D][DofQuadLimits::MAX_D1D];
1586  for (int dy = 0; dy < D1D; ++dy)
1587  {
1588  for (int dx = 0; dx < D1D; ++dx)
1589  {
1590  sol_xy[dy][dx] = 0;
1591  }
1592  }
1593  for (int qy = 0; qy < Q1D; ++qy)
1594  {
1595  double sol_x[DofQuadLimits::MAX_D1D];
1596  for (int dx = 0; dx < D1D; ++dx)
1597  {
1598  sol_x[dx] = 0;
1599  }
1600  for (int qx = 0; qx < Q1D; ++qx)
1601  {
1602  const double s = m_(qx, qy, qz, e) * x_(qx, qy, qz, e);
1603  for (int dx = 0; dx < D1D; ++dx)
1604  {
1605  sol_x[dx] += Bt_(dx, qx) * s;
1606  }
1607  }
1608  for (int dy = 0; dy < D1D; ++dy)
1609  {
1610  const double wy = Bt_(dy, qy);
1611  for (int dx = 0; dx < D1D; ++dx)
1612  {
1613  sol_xy[dy][dx] += wy * sol_x[dx];
1614  }
1615  }
1616  }
1617  for (int dz = 0; dz < D1D; ++dz)
1618  {
1619  const double wz = Bt_(dz, qz);
1620  for (int dy = 0; dy < D1D; ++dy)
1621  {
1622  for (int dx = 0; dx < D1D; ++dx)
1623  {
1624  y_(dx, dy, dz, e) += wz * sol_xy[dy][dx];
1625  }
1626  }
1627  }
1628  }
1629  });
1630 }
1631 } // namespace TransferKernels
1632 
1633 
1636 {
1637 }
1638 
1640  Vector& y) const
1641 {
1642  if (lFESpace.GetMesh()->GetNE() == 0)
1643  {
1644  return;
1645  }
1646 
1647  elem_restrict_lex_l->Mult(x, localL);
1648  if (dim == 2)
1649  {
1650  TransferKernels::Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
1651  }
1652  else if (dim == 3)
1653  {
1654  TransferKernels::Prolongation3D(NE, D1D, Q1D, localL, localH, B, mask);
1655  }
1656  else
1657  {
1658  MFEM_ABORT("TensorProductPRefinementTransferOperator::Mult not "
1659  "implemented for dim = "
1660  << dim);
1661  }
1662  elem_restrict_lex_h->MultTranspose(localH, y);
1663 }
1664 
1666  Vector& y) const
1667 {
1668  if (lFESpace.GetMesh()->GetNE() == 0)
1669  {
1670  return;
1671  }
1672 
1673  elem_restrict_lex_h->Mult(x, localH);
1674  if (dim == 2)
1675  {
1676  TransferKernels::Restriction2D(NE, D1D, Q1D, localH, localL, Bt, mask);
1677  }
1678  else if (dim == 3)
1679  {
1680  TransferKernels::Restriction3D(NE, D1D, Q1D, localH, localL, Bt, mask);
1681  }
1682  else
1683  {
1684  MFEM_ABORT("TensorProductPRefinementTransferOperator::MultTranspose not "
1685  "implemented for dim = "
1686  << dim);
1687  }
1688  elem_restrict_lex_l->MultTranspose(localL, y);
1689 }
1690 
1691 
1693  const FiniteElementSpace& hFESpace_)
1694  : Operator(hFESpace_.GetTrueVSize(), lFESpace_.GetTrueVSize()),
1695  lFESpace(lFESpace_),
1696  hFESpace(hFESpace_)
1697 {
1698  localTransferOperator = new TransferOperator(lFESpace_, hFESpace_);
1699 
1700  P = lFESpace.GetProlongationMatrix();
1701  R = hFESpace.IsVariableOrder() ? hFESpace.GetHpRestrictionMatrix() :
1702  hFESpace.GetRestrictionMatrix();
1703 
1704  // P and R can be both null
1705  // P can be null and R not null
1706  // If P is not null it is assumed that R is not null as well
1707  if (P) { MFEM_VERIFY(R, "Both P and R have to be not NULL") }
1708 
1709  if (P)
1710  {
1711  tmpL.SetSize(lFESpace_.GetVSize());
1712  tmpH.SetSize(hFESpace_.GetVSize());
1713  }
1714  // P can be null and R not null
1715  else if (R)
1716  {
1717  tmpH.SetSize(hFESpace_.GetVSize());
1718  }
1719 }
1720 
1722 {
1723  delete localTransferOperator;
1724 }
1725 
1726 void TrueTransferOperator::Mult(const Vector& x, Vector& y) const
1727 {
1728  if (P)
1729  {
1730  P->Mult(x, tmpL);
1731  localTransferOperator->Mult(tmpL, tmpH);
1732  R->Mult(tmpH, y);
1733  }
1734  else if (R)
1735  {
1736  localTransferOperator->Mult(x, tmpH);
1737  R->Mult(tmpH, y);
1738  }
1739  else
1740  {
1741  localTransferOperator->Mult(x, y);
1742  }
1743 }
1744 
1746 {
1747  if (P)
1748  {
1749  R->MultTranspose(x, tmpH);
1750  localTransferOperator->MultTranspose(tmpH, tmpL);
1751  P->MultTranspose(tmpL, y);
1752  }
1753  else if (R)
1754  {
1755  R->MultTranspose(x, tmpH);
1756  localTransferOperator->MultTranspose(tmpH, y);
1757  }
1758  else
1759  {
1760  localTransferOperator->MultTranspose(x, y);
1761  }
1762 }
1763 
1764 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
Abstract class for all finite elements.
Definition: fe_base.hpp:233
std::pair< std::unique_ptr< SparseMatrix >, std::unique_ptr< SparseMatrix > > ComputeSparseRAndM_LH()
Computes on-rank R and M_LH matrices.
Definition: transfer.cpp:772
const Operator & MakeTrueOperator(FiniteElementSpace &fes_in, FiniteElementSpace &fes_out, const Operator &oper, OperatorHandle &t_oper)
Definition: transfer.cpp:34
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: operator.hpp:93
virtual void Prolongate(const Vector &x, Vector &y) const
Definition: transfer.cpp:473
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3146
int * GetJ()
Definition: table.hpp:114
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:577
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:1665
Data type for scaled Jacobi-type smoother of sparse matrix.
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
BilinearFormIntegrator * mass_integ
Ownership depends on own_mass_integ.
Definition: transfer.hpp:126
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:161
L2ProjectionH1Space(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
Definition: transfer.cpp:545
void SetRow(int r, const double *row)
Definition: densemat.cpp:2177
virtual ~TensorProductPRefinementTransferOperator()
Destructor.
Definition: transfer.cpp:1635
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
Field is discontinuous across element interfaces.
Definition: fe_coll.hpp:48
virtual void Mult(const Vector &x, Vector &y) const
Definition: transfer.cpp:651
static const int NumGeom
Definition: geom.hpp:42
FiniteElementSpace & dom_fes
Domain FE space.
Definition: transfer.hpp:32
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
virtual int GetContType() const =0
Matrix-free transfer operator between finite element spaces on the same mesh.
Definition: transfer.hpp:428
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3943
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition: table.cpp:199
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3909
static int DecodeDof(int dof)
Helper to return the DOF associated with a sign encoded DOF.
Definition: fespace.hpp:996
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1275
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
void SetDims(int rows, int nnz)
Definition: table.cpp:140
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:1745
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
void TransformDual(double *v) const
Definition: doftrans.hpp:189
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual const Operator * GetRestrictionOperator() const
An abstract operator that performs the same action as GetRestrictionMatrix.
Definition: fespace.hpp:612
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe_base.hpp:173
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual ~TransferOperator()
Destructor.
Definition: transfer.cpp:1193
Operator::Type oper_type
Desired Operator::Type for the construction of all operators defined by the underlying transfer algor...
Definition: transfer.hpp:38
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:706
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe_base.hpp:169
~TrueTransferOperator()
Destructor.
Definition: transfer.cpp:1721
virtual void Mult(const Vector &x, Vector &y) const override
Interpolation or prolongation of a vector x corresponding to the coarse space to the vector y corresp...
Definition: transfer.cpp:1195
void LoseData()
Call this if data has been stolen.
Definition: table.hpp:180
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual const Operator & ForwardOperator()
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
Definition: transfer.cpp:1113
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3897
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
void GetVDofs(int vd, Array< int > &dofs, int ndofs=-1) const
Returns the indices of all of the VDofs for the specified dimension &#39;vd&#39;.
Definition: fespace.cpp:195
void TDofsListByVDim(const FiniteElementSpace &fes, int vdim, Array< int > &vdofs_list) const
Fills the vdofs_list array with a list of vdofs for a given vdim and a given finite element space...
Definition: transfer.cpp:970
virtual void Prolongate(const Vector &x, Vector &y) const
Definition: transfer.cpp:701
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition: hypre.cpp:2853
virtual void ProlongateTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:510
int RowSize(int i) const
Definition: table.hpp:108
std::unique_ptr< SparseMatrix > AllocR()
Computes sparsity pattern and initializes R matrix.
Definition: transfer.cpp:1017
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 BuildHo2Lor(int nel_ho, int nel_lor, const CoarseFineTransformations &cf_tr)
Definition: transfer.cpp:238
virtual ~PRefinementTransferOperator()
Destructor.
Definition: transfer.cpp:1214
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition: fespace.cpp:621
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
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
Derefinement operator, used by the friend class InterpolationGridTransfer.
Definition: fespace.hpp:438
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1302
FiniteElementSpace & ran_fes
Range FE space.
Definition: transfer.hpp:33
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
ID for class SparseMatrix.
Definition: operator.hpp:286
const Table * GetElementToFaceOrientationTable() const
Definition: fespace.hpp:1093
virtual const DofToQuad & GetDofToQuad(const IntegrationRule &ir, DofToQuad::Mode mode) const
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:365
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
void SetFromTDofsTranspose(const FiniteElementSpace &fes, const Vector &X, Vector &x) const
Sets dual field coefficients given a vector of dual field coefficients on the tdofs and a finite elem...
Definition: transfer.cpp:956
virtual const Operator & BackwardOperator()
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
Definition: transfer.cpp:189
void SetPointMat(const DenseMatrix &pm)
Set the underlying point matrix describing the transformation.
Definition: eltrans.hpp:401
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
General product operator: x -> (A*B)(x) = A(B(x)).
Definition: operator.hpp:774
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:285
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:1200
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
DenseTensor point_matrices[Geometry::NumGeom]
Definition: ncmesh.hpp:77
Data type sparse matrix.
Definition: sparsemat.hpp:50
void ElemMixedMass(Geometry::Type geom, const FiniteElement &fe_ho, const FiniteElement &fe_lor, ElementTransformation *el_tr, IntegrationPointTransformation &ip_tr, DenseMatrix &M_mixed_el) const
Definition: transfer.cpp:258
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
OperatorHandle F
Forward, coarse-to-fine, operator.
Definition: transfer.hpp:129
GridTransfer(FiniteElementSpace &dom_fes_, FiniteElementSpace &ran_fes_)
Definition: transfer.cpp:19
void InvTransformPrimal(double *v) const
Definition: doftrans.hpp:182
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:616
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition: device.hpp:277
void SetFromTDofs(const FiniteElementSpace &fes, const Vector &X, Vector &x) const
Sets dof values given a vector of tdofs and a finite element space.
Definition: transfer.cpp:928
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
TensorProductPRefinementTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator from lFESpace to hFESpace which have different FE collections.
Definition: transfer.cpp:1341
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_BigInt *row_starts=NULL) const
Multiply the HypreParMatrix on the left by a block-diagonal parallel matrix D and return the result a...
Definition: hypre.cpp:1915
void Restriction2D(const int NE, const int D1D, const int Q1D, const Vector &localH, Vector &localL, const Array< double > &Bt, const Vector &mask)
Definition: transfer.cpp:1531
virtual void SetAbsTol(double p_atol_)
Sets absolute tolerance in preconditioned conjugate gradient solver.
Definition: transfer.cpp:764
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1268
Array< Embedding > embeddings
Fine element positions in their parents.
Definition: ncmesh.hpp:73
L2Projection(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
Definition: transfer.cpp:231
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition: fe_base.cpp:119
virtual void Mult(const Vector &x, Vector &y) const
Definition: transfer.cpp:406
void LumpedMassInverse(Vector &ML_inv) const
Returns the inverse of an on-rank lumped mass matrix.
Definition: transfer.cpp:991
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
Definition: densemat.cpp:3232
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:172
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition: vector.cpp:671
OperatorHandle B
Backward, fine-to-coarse, operator.
Definition: transfer.hpp:130
List of mesh geometries stored as Array<Geometry::Type>.
Definition: mesh.hpp:1273
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
virtual const Operator & BackwardOperator()
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
Definition: transfer.cpp:1119
void SetMassIntegrator(BilinearFormIntegrator *mass_integ_, bool own_mass_integ_=true)
Assign a mass integrator to be used in the construction of the backward, fine-to-coarse, transfer operator.
Definition: transfer.cpp:146
Array< double > Bt
Transpose of B.
Definition: fe_base.hpp:190
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:712
virtual void MultTranspose(const Vector &x, Vector &y) const
Multiply a vector with the transposed matrix. y = At * x.
Definition: sparsemat.cpp:908
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:1267
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1419
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition: eltrans.cpp:383
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
void forall(int N, lambda &&body)
Definition: forall.hpp:742
Matrix-free transfer operator between finite element spaces.
Definition: transfer.hpp:397
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:37
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const Table *coarse_elem_fos, const DenseTensor localP[]) const
Definition: fespace.cpp:1430
Matrix-free transfer operator between finite element spaces on the same mesh exploiting the tensor pr...
Definition: transfer.hpp:459
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1314
void TransformPrimal(double *v) const
Definition: doftrans.hpp:163
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
virtual void Mult(const Vector &x, Vector &y) const override
Interpolation or prolongation of a vector x corresponding to the coarse space to the vector y corresp...
Definition: transfer.cpp:1639
L2Prolongation * B
Backward, fine-to-coarse, operator.
Definition: transfer.hpp:375
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
void SetupPCG()
Sets up the PCG solver (sets parameters, operator, and preconditioner)
Definition: transfer.cpp:638
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:546
void GetTDofsTranspose(const FiniteElementSpace &fes, const Vector &x, Vector &X) const
Recovers a vector of dual field coefficients on the tdofs given a vector of dual coefficients and a f...
Definition: transfer.cpp:942
void Restriction3D(const int NE, const int D1D, const int Q1D, const Vector &localH, Vector &localL, const Array< double > &Bt, const Vector &mask)
Definition: transfer.cpp:1570
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
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:136
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;. The returned indices are offsets into an ldo...
Definition: fespace.cpp:2731
bool own_mass_integ
Ownership flag for mass_integ.
Definition: transfer.hpp:127
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:1611
const CoarseFineTransformations & GetRefinementTransforms()
Definition: mesh.cpp:10344
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:263
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
Definition: mesh.cpp:6285
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1026
virtual void MultTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:438
L2Projection * F
Forward, coarse-to-fine, operator.
Definition: transfer.hpp:374
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:184
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:593
virtual const Operator & ForwardOperator()
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
Definition: transfer.cpp:155
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:469
Class for integration point with weight.
Definition: intrules.hpp:31
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
GridFunction interpolation operator applicable after mesh refinement.
Definition: fespace.hpp:412
virtual void Mult(const Vector &x, Vector &y) const override
Interpolation or prolongation of a vector x corresponding to the coarse space to the vector y corresp...
Definition: transfer.cpp:1216
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
Definition: sparsemat.cpp:3767
const FiniteElementSpace & fes_lor
Definition: transfer.hpp:195
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition: operator.hpp:838
PRefinementTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator from lFESpace to hFESpace which have different FE collections.
Definition: transfer.cpp:1206
TrueTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator working on true degrees of freedom from lFESpace to hFESpace...
Definition: transfer.cpp:1692
IsoparametricTransformation Transf
Definition: eltrans.hpp:464
virtual void Mult(const Vector &x, Vector &y) const override
Interpolation or prolongation of a true dof vector x to a true dof vector y.
Definition: transfer.cpp:1726
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:864
void GetTDofs(const FiniteElementSpace &fes, const Vector &x, Vector &X) const
Recovers vector of tdofs given a vector of dofs and a finite element space.
Definition: transfer.cpp:914
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:349
Lexicographic ordering for tensor-product FiniteElements.
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
bool Parallel() const
Definition: transfer.hpp:46
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void InvTransformDual(double *v) const
Definition: doftrans.hpp:195
MPI_Comm GetComm() const
Definition: pfespace.hpp:269
static void ListToMarker(const Array< int > &list, int marker_size, Array< int > &marker, int mark_val=-1)
Convert an array of indices (list) to a Boolean marker array where all indices in the list are marked...
Definition: fespace.cpp:640
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
Vector data type.
Definition: vector.hpp:58
ID for class HypreParMatrix.
Definition: operator.hpp:287
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
int SizeK() const
Definition: densemat.hpp:1144
int * GetI()
Definition: table.hpp:113
virtual void ProlongateTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:730
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
Definition: densemat.cpp:1591
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:70
RefCoord s[3]
Field is continuous across element interfaces.
Definition: fe_coll.hpp:45
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:317
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
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
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:389
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:1096
L2ProjectionL2Space(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
Definition: transfer.cpp:285
HYPRE_BigInt * GetTrueDofOffsets() const
Definition: pfespace.hpp:278
void Prolongation2D(const int NE, const int D1D, const int Q1D, const Vector &localL, Vector &localH, const Array< double > &B, const Vector &mask)
Definition: transfer.cpp:1412
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
virtual const SparseMatrix * GetHpRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:620
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
Definition: fespace.cpp:1494
virtual void SetRelTol(double p_rtol_)
Sets relative tolerance in preconditioned conjugate gradient solver.
Definition: transfer.cpp:759
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
Definition: fespace.cpp:215
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
virtual bool SupportsBackwardsOperator() const
Definition: transfer.cpp:1155
TransferOperator(const FiniteElementSpace &lFESpace, const FiniteElementSpace &hFESpace)
Constructs a transfer operator from lFESpace to hFESpace.
Definition: transfer.cpp:1161
void GetTransferOperator(const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
Construct and return an Operator that can be used to transfer GridFunction data from coarse_fes...
Definition: fespace.cpp:3312
virtual void MultTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:676
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
void Prolongation3D(const int NE, const int D1D, const int Q1D, const Vector &localL, Vector &localH, const Array< double > &B, const Vector &mask)
Definition: transfer.cpp:1459
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:769