MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
transfer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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  // If the local mesh is empty, skip all computations
295  if (nel_ho == 0) { return; }
296 
297  const CoarseFineTransformations &cf_tr = mesh_lor->GetRefinementTransforms();
298 
299  int nref_max = 0;
300  Array<Geometry::Type> geoms;
301  mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
302  for (int ig = 0; ig < geoms.Size(); ++ig)
303  {
304  Geometry::Type geom = geoms[ig];
305  nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
306  }
307 
308  BuildHo2Lor(nel_ho, nel_lor, cf_tr);
309 
310  offsets.SetSize(nel_ho+1);
311  offsets[0] = 0;
312  for (int iho = 0; iho < nel_ho; ++iho)
313  {
314  int nref = ho2lor.RowSize(iho);
315  const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
316  const FiniteElement &fe_lor = *fes_lor.GetFE(ho2lor.GetRow(iho)[0]);
317  offsets[iho+1] = offsets[iho] + fe_ho.GetDof()*fe_lor.GetDof()*nref;
318  }
319  // R will contain the restriction (L^2 projection operator) defined on each
320  // coarse HO element (and corresponding patch of LOR elements)
321  R.SetSize(offsets[nel_ho]);
322  // P will contain the corresponding prolongation operator
323  P.SetSize(offsets[nel_ho]);
324 
326  IsoparametricTransformation &emb_tr = ip_tr.Transf;
327 
328  for (int iho = 0; iho < nel_ho; ++iho)
329  {
330  Array<int> lor_els;
331  ho2lor.GetRow(iho, lor_els);
332  int nref = ho2lor.RowSize(iho);
333 
334  Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
335  const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
336  const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
337  int ndof_ho = fe_ho.GetDof();
338  int ndof_lor = fe_lor.GetDof();
339 
340  emb_tr.SetIdentityTransformation(geom);
341  const DenseTensor &pmats = cf_tr.point_matrices[geom];
342 
343  DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
344  DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
345 
346  DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
347  DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
348 
349  MassIntegrator mi;
350  DenseMatrix M_lor_el(ndof_lor, ndof_lor);
351  DenseMatrixInverse Minv_lor_el(&M_lor_el);
352  DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
353  DenseMatrix M_mixed_el(ndof_lor, ndof_ho);
354 
355  Minv_lor = 0.0;
356  M_lor = 0.0;
357 
358  DenseMatrix RtMlor(ndof_ho, ndof_lor*nref);
359  DenseMatrix RtMlorR(ndof_ho, ndof_ho);
360  DenseMatrixInverse RtMlorR_inv(&RtMlorR);
361 
362  for (int iref = 0; iref < nref; ++iref)
363  {
364  // Assemble the low-order refined mass matrix and invert locally
365  int ilor = lor_els[iref];
367  mi.AssembleElementMatrix(fe_lor, *el_tr, M_lor_el);
368  M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
369  Minv_lor_el.Factor();
370  Minv_lor_el.GetInverseMatrix(M_lor_el);
371  // Insert into the diagonal of the patch LOR mass matrix
372  Minv_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
373 
374  // Now assemble the block-row of the mixed mass matrix associated
375  // with integrating HO functions against LOR functions on the LOR
376  // sub-element.
377 
378  // Create the transformation that embeds the fine low-order element
379  // within the coarse high-order element in reference space
380  emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
381 
382  ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, M_mixed_el);
383 
384  M_mixed.CopyMN(M_mixed_el, iref*ndof_lor, 0);
385  }
386  mfem::Mult(Minv_lor, M_mixed, R_iho);
387 
388  mfem::MultAtB(R_iho, M_lor, RtMlor);
389  mfem::Mult(RtMlor, R_iho, RtMlorR);
390  RtMlorR_inv.Factor();
391  RtMlorR_inv.Mult(RtMlor, P_iho);
392  }
393 }
394 
396  const Vector &x, Vector &y) const
397 {
398  int vdim = fes_ho.GetVDim();
399  Array<int> vdofs;
400  DenseMatrix xel_mat, yel_mat;
401  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
402  {
403  int nref = ho2lor.RowSize(iho);
404  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
405  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
406  xel_mat.SetSize(ndof_ho, vdim);
407  yel_mat.SetSize(ndof_lor*nref, vdim);
408  DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
409 
410  fes_ho.GetElementVDofs(iho, vdofs);
411  x.GetSubVector(vdofs, xel_mat.GetData());
412  mfem::Mult(R_iho, xel_mat, yel_mat);
413  // Place result correctly into the low-order vector
414  for (int iref = 0; iref < nref; ++iref)
415  {
416  int ilor = ho2lor.GetRow(iho)[iref];
417  for (int vd=0; vd<vdim; ++vd)
418  {
419  fes_lor.GetElementDofs(ilor, vdofs);
420  fes_lor.DofsToVDofs(vd, vdofs);
421  y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
422  }
423  }
424  }
425 }
426 
428  const Vector &x, Vector &y) const
429 {
430  int vdim = fes_ho.GetVDim();
431  Array<int> vdofs;
432  DenseMatrix xel_mat, yel_mat;
433  y = 0.0;
434  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
435  {
436  int nref = ho2lor.RowSize(iho);
437  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
438  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
439  xel_mat.SetSize(ndof_lor*nref, vdim);
440  yel_mat.SetSize(ndof_ho, vdim);
441  DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
442 
443  // Extract the LOR DOFs
444  for (int iref=0; iref<nref; ++iref)
445  {
446  int ilor = ho2lor.GetRow(iho)[iref];
447  for (int vd=0; vd<vdim; ++vd)
448  {
449  fes_lor.GetElementDofs(ilor, vdofs);
450  fes_lor.DofsToVDofs(vd, vdofs);
451  x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
452  }
453  }
454  // Multiply locally by the transpose
455  mfem::MultAtB(R_iho, xel_mat, yel_mat);
456  // Place the result in the HO vector
457  fes_ho.GetElementVDofs(iho, vdofs);
458  y.AddElementVector(vdofs, yel_mat.GetData());
459  }
460 }
461 
463  const Vector &x, Vector &y) const
464 {
465  int vdim = fes_ho.GetVDim();
466  Array<int> vdofs;
467  DenseMatrix xel_mat,yel_mat;
468  y = 0.0;
469  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
470  {
471  int nref = ho2lor.RowSize(iho);
472  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
473  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
474  xel_mat.SetSize(ndof_lor*nref, vdim);
475  yel_mat.SetSize(ndof_ho, vdim);
476  DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
477 
478  // Extract the LOR DOFs
479  for (int iref = 0; iref < nref; ++iref)
480  {
481  int ilor = ho2lor.GetRow(iho)[iref];
482  for (int vd = 0; vd < vdim; ++vd)
483  {
484  fes_lor.GetElementDofs(ilor, vdofs);
485  fes_lor.DofsToVDofs(vd, vdofs);
486  x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
487  }
488  }
489  // Locally prolongate
490  mfem::Mult(P_iho, xel_mat, yel_mat);
491  // Place the result in the HO vector
492  fes_ho.GetElementVDofs(iho, vdofs);
493  y.AddElementVector(vdofs, yel_mat.GetData());
494  }
495 }
496 
498  const Vector &x, Vector &y) const
499 {
500  int vdim = fes_ho.GetVDim();
501  Array<int> vdofs;
502  DenseMatrix xel_mat,yel_mat;
503  for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
504  {
505  int nref = ho2lor.RowSize(iho);
506  int ndof_ho = fes_ho.GetFE(iho)->GetDof();
507  int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
508  xel_mat.SetSize(ndof_ho, vdim);
509  yel_mat.SetSize(ndof_lor*nref, vdim);
510  DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
511 
512  fes_ho.GetElementVDofs(iho, vdofs);
513  x.GetSubVector(vdofs, xel_mat.GetData());
514  mfem::MultAtB(P_iho, xel_mat, yel_mat);
515 
516  // Place result correctly into the low-order vector
517  for (int iref = 0; iref < nref; ++iref)
518  {
519  int ilor = ho2lor.GetRow(iho)[iref];
520  for (int vd=0; vd<vdim; ++vd)
521  {
522  fes_lor.GetElementDofs(ilor, vdofs);
523  fes_lor.DofsToVDofs(vd, vdofs);
524  y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
525  }
526  }
527  }
528 }
529 
531  const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_)
532  : L2Projection(fes_ho_, fes_lor_)
533 {
534  Mesh* mesh_ho = fes_ho.GetMesh();
535  Mesh* mesh_lor = fes_lor.GetMesh();
536  int nel_ho = mesh_ho->GetNE();
537  int nel_lor = mesh_lor->GetNE();
538  int ndof_lor = fes_lor.GetNDofs();
539 
540  // If the local mesh is empty, skip all computations
541  if (nel_ho == 0) { return; }
542 
543  const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
544 
545  int nref_max = 0;
546  Array<Geometry::Type> geoms;
547  mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
548  for (int ig = 0; ig < geoms.Size(); ++ig)
549  {
550  Geometry::Type geom = geoms[ig];
551  nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
552  }
553 
554  BuildHo2Lor(nel_ho, nel_lor, cf_tr);
555 
556  // ML_inv contains the inverse lumped (row sum) mass matrix. Note that the
557  // method will also work with a full (consistent) mass matrix, though this is
558  // not implemented here. L refers to the low-order refined mesh
559  Vector ML_inv(ndof_lor);
560  ML_inv = 0.0;
561 
562  // Compute ML_inv
563  for (int iho = 0; iho < nel_ho; ++iho)
564  {
565  Array<int> lor_els;
566  ho2lor.GetRow(iho, lor_els);
567  int nref = ho2lor.RowSize(iho);
568 
569  Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
570  const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
571  int nedof_lor = fe_lor.GetDof();
572 
573  // Instead of using a MassIntegrator, manually loop over integration
574  // points so we can row sum and store the diagonal as a Vector.
575  Vector ML_el(nedof_lor);
576  Vector shape_lor(nedof_lor);
577  Array<int> dofs_lor(nedof_lor);
578 
579  for (int iref = 0; iref < nref; ++iref)
580  {
581  int ilor = lor_els[iref];
583 
584  int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
585  const IntegrationRule* ir = &IntRules.Get(geom, order);
586  ML_el = 0.0;
587  for (int i = 0; i < ir->GetNPoints(); ++i)
588  {
589  const IntegrationPoint& ip_lor = ir->IntPoint(i);
590  fe_lor.CalcShape(ip_lor, shape_lor);
591  el_tr->SetIntPoint(&ip_lor);
592  ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
593  }
594  fes_lor.GetElementDofs(ilor, dofs_lor);
595  ML_inv.AddElementVector(dofs_lor, ML_el);
596  }
597  }
598  // DOF by DOF inverse of non-zero entries
599  for (int i = 0; i < ndof_lor; ++i)
600  {
601  ML_inv[i] = 1.0 / ML_inv[i];
602  }
603 
604  // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate
605  AllocR();
606  // Allocate M_LH (same sparsity pattern as R)
607  // L refers to the low-order refined mesh (DOFs correspond to rows)
608  // H refers to the higher-order mesh (DOFs correspond to columns)
609  M_LH = SparseMatrix(R.GetI(), R.GetJ(), NULL,
610  R.Height(), R.Width(), false, true, true);
611 
613  IsoparametricTransformation& emb_tr = ip_tr.Transf;
614 
615  // Compute M_LH and R
616  for (int iho = 0; iho < nel_ho; ++iho)
617  {
618  Array<int> lor_els;
619  ho2lor.GetRow(iho, lor_els);
620  int nref = ho2lor.RowSize(iho);
621 
622  Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
623  const FiniteElement& fe_ho = *fes_ho.GetFE(iho);
624  const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
625 
626  emb_tr.SetIdentityTransformation(geom);
627  const DenseTensor& pmats = cf_tr.point_matrices[geom];
628 
629  int nedof_ho = fe_ho.GetDof();
630  int nedof_lor = fe_lor.GetDof();
631  DenseMatrix M_LH_el(nedof_lor, nedof_ho);
632  DenseMatrix R_el(nedof_lor, nedof_ho);
633 
634  for (int iref = 0; iref < nref; ++iref)
635  {
636  int ilor = lor_els[iref];
638 
639  // Create the transformation that embeds the fine low-order element
640  // within the coarse high-order element in reference space
641  emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
642 
643  ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, M_LH_el);
644 
645  Array<int> dofs_lor(nedof_lor);
646  fes_lor.GetElementDofs(ilor, dofs_lor);
647  Vector R_row;
648  for (int i = 0; i < nedof_lor; ++i)
649  {
650  M_LH_el.GetRow(i, R_row);
651  R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
652  }
653  Array<int> dofs_ho(nedof_ho);
654  fes_ho.GetElementDofs(iho, dofs_ho);
655  M_LH.AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
656  R.AddSubMatrix(dofs_lor, dofs_ho, R_el);
657  }
658  }
659 
660  // Create PCG solver
661  RTxM_LH = TransposeMult(R, M_LH);
662  pcg.SetPrintLevel(0);
663  pcg.SetMaxIter(1000);
664  // initial values for relative and absolute tolerance
665  SetRelTol(1e-13);
666  SetAbsTol(1e-13);
667  Ds = DSmoother(*RTxM_LH);
668  pcg.SetPreconditioner(Ds);
669  pcg.SetOperator(*RTxM_LH);
670 }
671 
673 {
674  delete RTxM_LH;
675 }
676 
678  const Vector& x, Vector& y) const
679 {
680  int vdim = fes_ho.GetVDim();
681  const int ndof_ho = fes_ho.GetNDofs();
682  const int ndof_lor = fes_lor.GetNDofs();
683  Array<int> dofs_ho(ndof_ho);
684  Array<int> dofs_lor(ndof_lor);
685  Vector x_dim(ndof_ho);
686  Vector y_dim(ndof_lor);
687 
688  for (int d = 0; d < vdim; ++d)
689  {
690  fes_ho.GetVDofs(d, dofs_ho);
691  fes_lor.GetVDofs(d, dofs_lor);
692  x.GetSubVector(dofs_ho, x_dim);
693  R.Mult(x_dim, y_dim);
694  y.SetSubVector(dofs_lor, y_dim);
695  }
696 }
697 
699  const Vector& x, Vector& y) const
700 {
701  int vdim = fes_ho.GetVDim();
702  const int ndof_ho = fes_ho.GetNDofs();
703  const int ndof_lor = fes_lor.GetNDofs();
704  Array<int> dofs_ho(ndof_ho);
705  Array<int> dofs_lor(ndof_lor);
706  Vector x_dim(ndof_lor);
707  Vector y_dim(ndof_ho);
708 
709  for (int d = 0; d < vdim; ++d)
710  {
711  fes_ho.GetVDofs(d, dofs_ho);
712  fes_lor.GetVDofs(d, dofs_lor);
713  x.GetSubVector(dofs_lor, x_dim);
714  R.MultTranspose(x_dim, y_dim);
715  y.SetSubVector(dofs_ho, y_dim);
716  }
717 }
718 
720  const Vector& x, Vector& y) const
721 {
722  int vdim = fes_ho.GetVDim();
723  const int ndof_ho = fes_ho.GetNDofs();
724  const int ndof_lor = fes_lor.GetNDofs();
725  Array<int> dofs_ho(ndof_ho);
726  Array<int> dofs_lor(ndof_lor);
727  Vector x_dim(ndof_lor);
728  Vector y_dim(ndof_ho);
729  Vector xbar(ndof_ho);
730 
731  for (int d = 0; d < vdim; ++d)
732  {
733  fes_lor.GetVDofs(d, dofs_lor);
734  x.GetSubVector(dofs_lor, x_dim);
735  // Compute y = P x = (R^T M_LH)^(-1) M_LH^T x = (R^T M_LH)^(-1) xbar
736  M_LH.MultTranspose(x_dim, xbar);
737  y_dim = 0.0;
738  pcg.Mult(xbar, y_dim);
739  fes_ho.GetVDofs(d, dofs_ho);
740  y.SetSubVector(dofs_ho, y_dim);
741  }
742 }
743 
745  const Vector& x, Vector& y) const
746 {
747  int vdim = fes_ho.GetVDim();
748  const int ndof_ho = fes_ho.GetNDofs();
749  const int ndof_lor = fes_lor.GetNDofs();
750  Array<int> dofs_ho(ndof_ho);
751  Array<int> dofs_lor(ndof_lor);
752  Vector x_dim(ndof_ho);
753  Vector y_dim(ndof_lor);
754  Vector xbar(ndof_ho);
755 
756  for (int d = 0; d < vdim; ++d)
757  {
758  fes_ho.GetVDofs(d, dofs_ho);
759  x.GetSubVector(dofs_ho, x_dim);
760  // Compute y = P^T x = M_LH (R^T M_LH)^(-1) x = M_LH xbar
761  xbar = 0.0;
762  pcg.Mult(x_dim, xbar);
763  M_LH.Mult(xbar, y_dim);
764  fes_lor.GetVDofs(d, dofs_lor);
765  y.SetSubVector(dofs_lor, y_dim);
766  }
767 }
768 
770 {
771  pcg.SetRelTol(p_rtol_);
772 }
773 
775 {
776  pcg.SetAbsTol(p_atol_);
777 }
778 
779 void L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR()
780 {
781  const Table& elem_dof_ho = fes_ho.GetElementToDofTable();
782  const Table& elem_dof_lor = fes_lor.GetElementToDofTable();
783  const int ndof_ho = fes_ho.GetNDofs();
784  const int ndof_lor = fes_lor.GetNDofs();
785 
786  Table dof_elem_lor;
787  Transpose(elem_dof_lor, dof_elem_lor, ndof_lor);
788 
789  Mesh* mesh_lor = fes_lor.GetMesh();
790  const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
791 
792  // mfem::Mult but uses ho2lor to map HO elements to LOR elements
793  const int* elem_dof_hoI = elem_dof_ho.GetI();
794  const int* elem_dof_hoJ = elem_dof_ho.GetJ();
795  const int* dof_elem_lorI = dof_elem_lor.GetI();
796  const int* dof_elem_lorJ = dof_elem_lor.GetJ();
797 
798  Array<int> I(ndof_lor + 1);
799 
800  // figure out the size of J
801  Array<int> dof_used_ho;
802  dof_used_ho.SetSize(ndof_ho, -1);
803 
804  int sizeJ = 0;
805  for (int ilor = 0; ilor < ndof_lor; ++ilor)
806  {
807  for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
808  {
809  int el_lor = dof_elem_lorJ[jlor];
810  int iho = cf_tr.embeddings[el_lor].parent;
811  for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
812  {
813  int dof_ho = elem_dof_hoJ[jho];
814  if (dof_used_ho[dof_ho] != ilor)
815  {
816  dof_used_ho[dof_ho] = ilor;
817  ++sizeJ;
818  }
819  }
820  }
821  }
822 
823  // initialize dof_ho_dof_lor
824  Table dof_lor_dof_ho;
825  dof_lor_dof_ho.SetDims(ndof_lor, sizeJ);
826 
827  for (int i = 0; i < ndof_ho; ++i)
828  {
829  dof_used_ho[i] = -1;
830  }
831 
832  // set values of J
833  int* dof_dofI = dof_lor_dof_ho.GetI();
834  int* dof_dofJ = dof_lor_dof_ho.GetJ();
835  sizeJ = 0;
836  for (int ilor = 0; ilor < ndof_lor; ++ilor)
837  {
838  dof_dofI[ilor] = sizeJ;
839  for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
840  {
841  int el_lor = dof_elem_lorJ[jlor];
842  int iho = cf_tr.embeddings[el_lor].parent;
843  for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
844  {
845  int dof_ho = elem_dof_hoJ[jho];
846  if (dof_used_ho[dof_ho] != ilor)
847  {
848  dof_used_ho[dof_ho] = ilor;
849  dof_dofJ[sizeJ] = dof_ho;
850  ++sizeJ;
851  }
852  }
853  }
854  }
855 
856  dof_lor_dof_ho.SortRows();
857  double* data = Memory<double>(dof_dofI[ndof_lor]);
858 
859  R = SparseMatrix(dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho,
860  true, true, true);
861  R = 0.0;
862 
863  dof_lor_dof_ho.LoseData();
864 }
865 
867 {
868  delete F;
869  delete B;
870 }
871 
873 {
874  if (!F) { BuildF(); }
875  return *F;
876 }
877 
879 {
880  if (!B)
881  {
882  if (!F) { BuildF(); }
883  B = new L2Prolongation(*F);
884  }
885  return *B;
886 }
887 
888 void L2ProjectionGridTransfer::BuildF()
889 {
890  if (!force_l2_space &&
892  {
893  F = new L2ProjectionH1Space(dom_fes, ran_fes);
894  }
895  else
896  {
897  F = new L2ProjectionL2Space(dom_fes, ran_fes);
898  }
899 }
900 
901 
903  const FiniteElementSpace& hFESpace_)
904  : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize())
905 {
906  bool isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
907  if (lFESpace_.FEColl() == hFESpace_.FEColl() && !isvar_order)
908  {
910  hFESpace_.GetTransferOperator(lFESpace_, P);
911  P.SetOperatorOwner(false);
912  opr = P.Ptr();
913  }
914  else if (lFESpace_.GetMesh()->GetNE() > 0
915  && hFESpace_.GetMesh()->GetNE() > 0
916  && lFESpace_.GetVDim() == 1
917  && hFESpace_.GetVDim() == 1
918  && dynamic_cast<const TensorBasisElement*>(lFESpace_.GetFE(0))
919  && dynamic_cast<const TensorBasisElement*>(hFESpace_.GetFE(0))
920  && !isvar_order
921  && (hFESpace_.FEColl()->GetContType() ==
923  hFESpace_.FEColl()->GetContType() ==
925  {
926  opr = new TensorProductPRefinementTransferOperator(lFESpace_, hFESpace_);
927  }
928  else
929  {
930  opr = new PRefinementTransferOperator(lFESpace_, hFESpace_);
931  }
932 }
933 
935 
936 void TransferOperator::Mult(const Vector& x, Vector& y) const
937 {
938  opr->Mult(x, y);
939 }
940 
942 {
943  opr->MultTranspose(x, y);
944 }
945 
946 
948  const FiniteElementSpace& lFESpace_, const FiniteElementSpace& hFESpace_)
949  : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
950  hFESpace(hFESpace_)
951 {
952  isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
953 }
954 
956 
958 {
959  Mesh* mesh = hFESpace.GetMesh();
960  Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
961  DenseMatrix loc_prol;
962  Vector subY, subX;
963 
964  Geometry::Type cached_geom = Geometry::INVALID;
965  const FiniteElement* h_fe = NULL;
966  const FiniteElement* l_fe = NULL;
968 
969  int vdim = lFESpace.GetVDim();
970 
971  for (int i = 0; i < mesh->GetNE(); i++)
972  {
973  DofTransformation * doftrans_h = hFESpace.GetElementDofs(i, h_dofs);
974  DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
975 
976  const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
977  if (geom != cached_geom || isvar_order)
978  {
979  h_fe = hFESpace.GetFE(i);
980  l_fe = lFESpace.GetFE(i);
982  h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
983  subY.SetSize(loc_prol.Height());
984  cached_geom = geom;
985  }
986 
987  for (int vd = 0; vd < vdim; vd++)
988  {
989  l_dofs.Copy(l_vdofs);
990  lFESpace.DofsToVDofs(vd, l_vdofs);
991  h_dofs.Copy(h_vdofs);
992  hFESpace.DofsToVDofs(vd, h_vdofs);
993  x.GetSubVector(l_vdofs, subX);
994  if (doftrans_l)
995  {
996  doftrans_l->InvTransformPrimal(subX);
997  }
998  loc_prol.Mult(subX, subY);
999  if (doftrans_h)
1000  {
1001  doftrans_h->TransformPrimal(subY);
1002  }
1003  y.SetSubVector(h_vdofs, subY);
1004  }
1005  }
1006 }
1007 
1009  Vector& y) const
1010 {
1011  y = 0.0;
1012 
1013  Mesh* mesh = hFESpace.GetMesh();
1014  Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
1015  DenseMatrix loc_prol;
1016  Vector subY, subX;
1017 
1018  Array<char> processed(hFESpace.GetVSize());
1019  processed = 0;
1020 
1021  Geometry::Type cached_geom = Geometry::INVALID;
1022  const FiniteElement* h_fe = NULL;
1023  const FiniteElement* l_fe = NULL;
1025 
1026  int vdim = lFESpace.GetVDim();
1027 
1028  for (int i = 0; i < mesh->GetNE(); i++)
1029  {
1030  DofTransformation * doftrans_h = hFESpace.GetElementDofs(i, h_dofs);
1031  DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
1032 
1033  const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
1034  if (geom != cached_geom || isvar_order)
1035  {
1036  h_fe = hFESpace.GetFE(i);
1037  l_fe = lFESpace.GetFE(i);
1039  h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
1040  loc_prol.Transpose();
1041  subY.SetSize(loc_prol.Height());
1042  cached_geom = geom;
1043  }
1044 
1045  for (int vd = 0; vd < vdim; vd++)
1046  {
1047  l_dofs.Copy(l_vdofs);
1048  lFESpace.DofsToVDofs(vd, l_vdofs);
1049  h_dofs.Copy(h_vdofs);
1050  hFESpace.DofsToVDofs(vd, h_vdofs);
1051 
1052  x.GetSubVector(h_vdofs, subX);
1053  if (doftrans_h)
1054  {
1055  doftrans_h->InvTransformDual(subX);
1056  }
1057  for (int p = 0; p < h_dofs.Size(); ++p)
1058  {
1059  if (processed[lFESpace.DecodeDof(h_dofs[p])])
1060  {
1061  subX[p] = 0.0;
1062  }
1063  }
1064 
1065  loc_prol.Mult(subX, subY);
1066  if (doftrans_l)
1067  {
1068  doftrans_l->TransformDual(subY);
1069  }
1070  y.AddElementVector(l_vdofs, subY);
1071  }
1072 
1073  for (int p = 0; p < h_dofs.Size(); ++p)
1074  {
1075  processed[lFESpace.DecodeDof(h_dofs[p])] = 1;
1076  }
1077  }
1078 }
1079 
1080 
1083  const FiniteElementSpace& lFESpace_,
1084  const FiniteElementSpace& hFESpace_)
1085  : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
1086  hFESpace(hFESpace_)
1087 {
1088  // Assuming the same element type
1089  Mesh* mesh = lFESpace.GetMesh();
1090  dim = mesh->Dimension();
1091  if (mesh->GetNE() == 0)
1092  {
1093  return;
1094  }
1095  const FiniteElement& el = *lFESpace.GetFE(0);
1096 
1097  const TensorBasisElement* ltel =
1098  dynamic_cast<const TensorBasisElement*>(&el);
1099  MFEM_VERIFY(ltel, "Low order FE space must be tensor product space");
1100 
1101  const TensorBasisElement* htel =
1102  dynamic_cast<const TensorBasisElement*>(hFESpace.GetFE(0));
1103  MFEM_VERIFY(htel, "High order FE space must be tensor product space");
1104  const Array<int>& hdofmap = htel->GetDofMap();
1105 
1106  const IntegrationRule& ir = hFESpace.GetFE(0)->GetNodes();
1107  IntegrationRule irLex = ir;
1108 
1109  // The quadrature points, or equivalently, the dofs of the high order space
1110  // must be sorted in lexicographical order
1111  for (int i = 0; i < ir.GetNPoints(); ++i)
1112  {
1113  int j = hdofmap[i] >=0 ? hdofmap[i] : -1 - hdofmap[i];
1114  irLex.IntPoint(i) = ir.IntPoint(j);
1115  }
1116 
1117  NE = lFESpace.GetNE();
1118  const DofToQuad& maps = el.GetDofToQuad(irLex, DofToQuad::TENSOR);
1119 
1120  D1D = maps.ndof;
1121  Q1D = maps.nqpt;
1122  B = maps.B;
1123  Bt = maps.Bt;
1124 
1125  elem_restrict_lex_l =
1127 
1128  MFEM_VERIFY(elem_restrict_lex_l,
1129  "Low order ElementRestriction not available");
1130 
1131  elem_restrict_lex_h =
1133 
1134  MFEM_VERIFY(elem_restrict_lex_h,
1135  "High order ElementRestriction not available");
1136 
1137  localL.SetSize(elem_restrict_lex_l->Height(), Device::GetMemoryType());
1138  localH.SetSize(elem_restrict_lex_h->Height(), Device::GetMemoryType());
1139  localL.UseDevice(true);
1140  localH.UseDevice(true);
1141 
1142  MFEM_VERIFY(dynamic_cast<const ElementRestriction*>(elem_restrict_lex_h),
1143  "High order element restriction is of unsupported type");
1144 
1145  mask.SetSize(localH.Size(), Device::GetMemoryType());
1146  static_cast<const ElementRestriction*>(elem_restrict_lex_h)
1147  ->BooleanMask(mask);
1148  mask.UseDevice(true);
1149 }
1150 
1151 namespace TransferKernels
1152 {
1153 void Prolongation2D(const int NE, const int D1D, const int Q1D,
1154  const Vector& localL, Vector& localH,
1155  const Array<double>& B, const Vector& mask)
1156 {
1157  auto x_ = Reshape(localL.Read(), D1D, D1D, NE);
1158  auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, NE);
1159  auto B_ = Reshape(B.Read(), Q1D, D1D);
1160  auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
1161 
1162  localH = 0.0;
1163 
1164  MFEM_FORALL(e, NE,
1165  {
1166  for (int dy = 0; dy < D1D; ++dy)
1167  {
1168  double sol_x[MAX_Q1D];
1169  for (int qy = 0; qy < Q1D; ++qy)
1170  {
1171  sol_x[qy] = 0.0;
1172  }
1173  for (int dx = 0; dx < D1D; ++dx)
1174  {
1175  const double s = x_(dx, dy, e);
1176  for (int qx = 0; qx < Q1D; ++qx)
1177  {
1178  sol_x[qx] += B_(qx, dx) * s;
1179  }
1180  }
1181  for (int qy = 0; qy < Q1D; ++qy)
1182  {
1183  const double d2q = B_(qy, dy);
1184  for (int qx = 0; qx < Q1D; ++qx)
1185  {
1186  y_(qx, qy, e) += d2q * sol_x[qx];
1187  }
1188  }
1189  }
1190  for (int qy = 0; qy < Q1D; ++qy)
1191  {
1192  for (int qx = 0; qx < Q1D; ++qx)
1193  {
1194  y_(qx, qy, e) *= m_(qx, qy, e);
1195  }
1196  }
1197  });
1198 }
1199 
1200 void Prolongation3D(const int NE, const int D1D, const int Q1D,
1201  const Vector& localL, Vector& localH,
1202  const Array<double>& B, const Vector& mask)
1203 {
1204  auto x_ = Reshape(localL.Read(), D1D, D1D, D1D, NE);
1205  auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, Q1D, NE);
1206  auto B_ = Reshape(B.Read(), Q1D, D1D);
1207  auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
1208 
1209  localH = 0.0;
1210 
1211  MFEM_FORALL(e, NE,
1212  {
1213  for (int dz = 0; dz < D1D; ++dz)
1214  {
1215  double sol_xy[MAX_Q1D][MAX_Q1D];
1216  for (int qy = 0; qy < Q1D; ++qy)
1217  {
1218  for (int qx = 0; qx < Q1D; ++qx)
1219  {
1220  sol_xy[qy][qx] = 0.0;
1221  }
1222  }
1223  for (int dy = 0; dy < D1D; ++dy)
1224  {
1225  double sol_x[MAX_Q1D];
1226  for (int qx = 0; qx < Q1D; ++qx)
1227  {
1228  sol_x[qx] = 0;
1229  }
1230  for (int dx = 0; dx < D1D; ++dx)
1231  {
1232  const double s = x_(dx, dy, dz, e);
1233  for (int qx = 0; qx < Q1D; ++qx)
1234  {
1235  sol_x[qx] += B_(qx, dx) * s;
1236  }
1237  }
1238  for (int qy = 0; qy < Q1D; ++qy)
1239  {
1240  const double wy = B_(qy, dy);
1241  for (int qx = 0; qx < Q1D; ++qx)
1242  {
1243  sol_xy[qy][qx] += wy * sol_x[qx];
1244  }
1245  }
1246  }
1247  for (int qz = 0; qz < Q1D; ++qz)
1248  {
1249  const double wz = B_(qz, dz);
1250  for (int qy = 0; qy < Q1D; ++qy)
1251  {
1252  for (int qx = 0; qx < Q1D; ++qx)
1253  {
1254  y_(qx, qy, qz, e) += wz * sol_xy[qy][qx];
1255  }
1256  }
1257  }
1258  }
1259  for (int qz = 0; qz < Q1D; ++qz)
1260  {
1261  for (int qy = 0; qy < Q1D; ++qy)
1262  {
1263  for (int qx = 0; qx < Q1D; ++qx)
1264  {
1265  y_(qx, qy, qz, e) *= m_(qx, qy, qz, e);
1266  }
1267  }
1268  }
1269  });
1270 }
1271 
1272 void Restriction2D(const int NE, const int D1D, const int Q1D,
1273  const Vector& localH, Vector& localL,
1274  const Array<double>& Bt, const Vector& mask)
1275 {
1276  auto x_ = Reshape(localH.Read(), Q1D, Q1D, NE);
1277  auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, NE);
1278  auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
1279  auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
1280 
1281  localL = 0.0;
1282 
1283  MFEM_FORALL(e, NE,
1284  {
1285  for (int qy = 0; qy < Q1D; ++qy)
1286  {
1287  double sol_x[MAX_D1D];
1288  for (int dx = 0; dx < D1D; ++dx)
1289  {
1290  sol_x[dx] = 0.0;
1291  }
1292  for (int qx = 0; qx < Q1D; ++qx)
1293  {
1294  const double s = m_(qx, qy, e) * x_(qx, qy, e);
1295  for (int dx = 0; dx < D1D; ++dx)
1296  {
1297  sol_x[dx] += Bt_(dx, qx) * s;
1298  }
1299  }
1300  for (int dy = 0; dy < D1D; ++dy)
1301  {
1302  const double q2d = Bt_(dy, qy);
1303  for (int dx = 0; dx < D1D; ++dx)
1304  {
1305  y_(dx, dy, e) += q2d * sol_x[dx];
1306  }
1307  }
1308  }
1309  });
1310 }
1311 void Restriction3D(const int NE, const int D1D, const int Q1D,
1312  const Vector& localH, Vector& localL,
1313  const Array<double>& Bt, const Vector& mask)
1314 {
1315  auto x_ = Reshape(localH.Read(), Q1D, Q1D, Q1D, NE);
1316  auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, D1D, NE);
1317  auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
1318  auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
1319 
1320  localL = 0.0;
1321 
1322  MFEM_FORALL(e, NE,
1323  {
1324  for (int qz = 0; qz < Q1D; ++qz)
1325  {
1326  double sol_xy[MAX_D1D][MAX_D1D];
1327  for (int dy = 0; dy < D1D; ++dy)
1328  {
1329  for (int dx = 0; dx < D1D; ++dx)
1330  {
1331  sol_xy[dy][dx] = 0;
1332  }
1333  }
1334  for (int qy = 0; qy < Q1D; ++qy)
1335  {
1336  double sol_x[MAX_D1D];
1337  for (int dx = 0; dx < D1D; ++dx)
1338  {
1339  sol_x[dx] = 0;
1340  }
1341  for (int qx = 0; qx < Q1D; ++qx)
1342  {
1343  const double s = m_(qx, qy, qz, e) * x_(qx, qy, qz, e);
1344  for (int dx = 0; dx < D1D; ++dx)
1345  {
1346  sol_x[dx] += Bt_(dx, qx) * s;
1347  }
1348  }
1349  for (int dy = 0; dy < D1D; ++dy)
1350  {
1351  const double wy = Bt_(dy, qy);
1352  for (int dx = 0; dx < D1D; ++dx)
1353  {
1354  sol_xy[dy][dx] += wy * sol_x[dx];
1355  }
1356  }
1357  }
1358  for (int dz = 0; dz < D1D; ++dz)
1359  {
1360  const double wz = Bt_(dz, qz);
1361  for (int dy = 0; dy < D1D; ++dy)
1362  {
1363  for (int dx = 0; dx < D1D; ++dx)
1364  {
1365  y_(dx, dy, dz, e) += wz * sol_xy[dy][dx];
1366  }
1367  }
1368  }
1369  }
1370  });
1371 }
1372 } // namespace TransferKernels
1373 
1374 
1377 {
1378 }
1379 
1381  Vector& y) const
1382 {
1383  if (lFESpace.GetMesh()->GetNE() == 0)
1384  {
1385  return;
1386  }
1387 
1388  elem_restrict_lex_l->Mult(x, localL);
1389  if (dim == 2)
1390  {
1391  TransferKernels::Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
1392  }
1393  else if (dim == 3)
1394  {
1395  TransferKernels::Prolongation3D(NE, D1D, Q1D, localL, localH, B, mask);
1396  }
1397  else
1398  {
1399  MFEM_ABORT("TensorProductPRefinementTransferOperator::Mult not "
1400  "implemented for dim = "
1401  << dim);
1402  }
1403  elem_restrict_lex_h->MultTranspose(localH, y);
1404 }
1405 
1407  Vector& y) const
1408 {
1409  if (lFESpace.GetMesh()->GetNE() == 0)
1410  {
1411  return;
1412  }
1413 
1414  elem_restrict_lex_h->Mult(x, localH);
1415  if (dim == 2)
1416  {
1417  TransferKernels::Restriction2D(NE, D1D, Q1D, localH, localL, Bt, mask);
1418  }
1419  else if (dim == 3)
1420  {
1421  TransferKernels::Restriction3D(NE, D1D, Q1D, localH, localL, Bt, mask);
1422  }
1423  else
1424  {
1425  MFEM_ABORT("TensorProductPRefinementTransferOperator::MultTranspose not "
1426  "implemented for dim = "
1427  << dim);
1428  }
1429  elem_restrict_lex_l->MultTranspose(localL, y);
1430 }
1431 
1432 
1434  const FiniteElementSpace& hFESpace_)
1435  : Operator(hFESpace_.GetTrueVSize(), lFESpace_.GetTrueVSize()),
1436  lFESpace(lFESpace_),
1437  hFESpace(hFESpace_)
1438 {
1439  localTransferOperator = new TransferOperator(lFESpace_, hFESpace_);
1440 
1441  P = lFESpace.GetProlongationMatrix();
1442  R = hFESpace.IsVariableOrder() ? hFESpace.GetHpRestrictionMatrix() :
1443  hFESpace.GetRestrictionMatrix();
1444 
1445  // P and R can be both null
1446  // P can be null and R not null
1447  // If P is not null it is assumed that R is not null as well
1448  if (P) { MFEM_VERIFY(R, "Both P and R have to be not NULL") }
1449 
1450  if (P)
1451  {
1452  tmpL.SetSize(lFESpace_.GetVSize());
1453  tmpH.SetSize(hFESpace_.GetVSize());
1454  }
1455  // P can be null and R not null
1456  else if (R)
1457  {
1458  tmpH.SetSize(hFESpace_.GetVSize());
1459  }
1460 }
1461 
1463 {
1464  delete localTransferOperator;
1465 }
1466 
1467 void TrueTransferOperator::Mult(const Vector& x, Vector& y) const
1468 {
1469  if (P)
1470  {
1471  P->Mult(x, tmpL);
1472  localTransferOperator->Mult(tmpL, tmpH);
1473  R->Mult(tmpH, y);
1474  }
1475  else if (R)
1476  {
1477  localTransferOperator->Mult(x, tmpH);
1478  R->Mult(tmpH, y);
1479  }
1480  else
1481  {
1482  localTransferOperator->Mult(x, y);
1483  }
1484 }
1485 
1487 {
1488  if (P)
1489  {
1490  R->MultTranspose(x, tmpH);
1491  localTransferOperator->MultTranspose(tmpH, tmpL);
1492  P->MultTranspose(tmpL, y);
1493  }
1494  else if (R)
1495  {
1496  R->MultTranspose(x, tmpH);
1497  localTransferOperator->MultTranspose(tmpH, y);
1498  }
1499  else
1500  {
1501  localTransferOperator->MultTranspose(x, y);
1502  }
1503 }
1504 
1505 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
const Operator & MakeTrueOperator(FiniteElementSpace &fes_in, FiniteElementSpace &fes_out, const Operator &oper, OperatorHandle &t_oper)
Definition: transfer.cpp:34
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:573
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:104
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3127
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const Table *coarse_elem_fos, const DenseTensor localP[]) const
Definition: fespace.cpp:1387
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
int * GetJ()
Definition: table.hpp:114
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:1406
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1232
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:463
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:923
virtual void Prolongate(const Vector &x, Vector &y) const
Definition: transfer.cpp:462
virtual const SparseMatrix * GetHpRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:499
BilinearFormIntegrator * mass_integ
Ownership depends on own_mass_integ.
Definition: transfer.hpp:124
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:162
L2ProjectionH1Space(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
Definition: transfer.cpp:530
void SetRow(int r, const double *row)
Definition: densemat.cpp:2139
virtual ~TensorProductPRefinementTransferOperator()
Destructor.
Definition: transfer.cpp:1376
Field is discontinuous across element interfaces.
Definition: fe_coll.hpp:48
static const int NumGeom
Definition: geom.hpp:42
FiniteElementSpace & dom_fes
Domain FE space.
Definition: transfer.hpp:32
virtual int GetContType() const =0
Matrix-free transfer operator between finite element spaces on the same mesh.
Definition: transfer.hpp:385
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
static int DecodeDof(int dof)
Helpers to remove encoded sign from a DOF.
Definition: fespace.hpp:262
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
virtual void Mult(const Vector &x, Vector &y) const
Definition: transfer.cpp:677
int SizeK() const
Definition: densemat.hpp:1001
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
int * GetJ()
Return the array J.
Definition: sparsemat.hpp:227
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:1486
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:856
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
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
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:5919
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe_base.hpp:174
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
virtual ~TransferOperator()
Destructor.
Definition: transfer.cpp:934
Operator::Type oper_type
Desired Operator::Type for the construction of all operators defined by the underlying transfer algor...
Definition: transfer.hpp:38
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe_base.hpp:170
int * GetI()
Return the array I.
Definition: sparsemat.hpp:222
~TrueTransferOperator()
Destructor.
Definition: transfer.cpp:1462
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
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
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:936
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:349
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Definition: densemat.cpp:3890
virtual const Operator & ForwardOperator()
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
Definition: transfer.cpp:872
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3878
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition: hypre.cpp:2849
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1069
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:118
void BuildHo2Lor(int nel_ho, int nel_lor, const CoarseFineTransformations &cf_tr)
Definition: transfer.cpp:238
virtual ~PRefinementTransferOperator()
Destructor.
Definition: transfer.cpp:955
virtual void Mult(const Vector &x, Vector &y) const
Definition: transfer.cpp:395
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:94
virtual void ProlongateTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:497
Derefinement operator, used by the friend class InterpolationGridTransfer.
Definition: fespace.hpp:324
FiniteElementSpace & ran_fes
Range FE space.
Definition: transfer.hpp:33
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1239
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
ID for class SparseMatrix.
Definition: operator.hpp:259
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
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 -&gt; (A*B)(x) = A(B(x)).
Definition: operator.hpp:747
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:258
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:941
DenseTensor point_matrices[Geometry::NumGeom]
Definition: ncmesh.hpp:77
Data type sparse matrix.
Definition: sparsemat.hpp:50
const int MAX_Q1D
Definition: forall.hpp:29
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
OperatorHandle F
Forward, coarse-to-fine, operator.
Definition: transfer.hpp:127
virtual void MultTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:698
GridTransfer(FiniteElementSpace &dom_fes_, FiniteElementSpace &ran_fes_)
Definition: transfer.cpp:19
void GetRow(int r, Vector &row) const
Definition: densemat.cpp:1276
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition: device.hpp:277
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
TensorProductPRefinementTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator from lFESpace to hFESpace which have different FE collections.
Definition: transfer.cpp:1082
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 * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:321
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:479
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
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:1272
Array< Embedding > embeddings
Fine element positions in their parents.
Definition: ncmesh.hpp:73
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:304
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
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
int Dimension() const
Definition: mesh.hpp:1006
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:639
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
Definition: sparsemat.cpp:2718
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
virtual void Prolongate(const Vector &x, Vector &y) const
Definition: transfer.cpp:719
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2678
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:3244
OperatorHandle B
Backward, fine-to-coarse, operator.
Definition: transfer.hpp:128
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
Definition: fespace.cpp:1451
List of mesh geometries stored as Array&lt;Geometry::Type&gt;.
Definition: mesh.hpp:1092
virtual const Operator & BackwardOperator()
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
Definition: transfer.cpp:878
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
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:647
Array< double > Bt
Transpose of B.
Definition: fe_base.hpp:191
HYPRE_BigInt * GetTrueDofOffsets() const
Definition: pfespace.hpp:282
virtual void InvTransformPrimal(double *v) const =0
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:35
virtual void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Definition: transfer.cpp:1008
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1381
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.
Matrix-free transfer operator between finite element spaces.
Definition: transfer.hpp:354
Operator that converts FiniteElementSpace L-vectors to E-vectors.
Definition: restriction.hpp:36
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
virtual void MultTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:427
Matrix-free transfer operator between finite element spaces on the same mesh exploiting the tensor pr...
Definition: transfer.hpp:416
virtual void TransformPrimal(double *v) const =0
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
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:1380
L2Prolongation * B
Backward, fine-to-coarse, operator.
Definition: transfer.hpp:334
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:546
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:1311
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3924
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:136
bool own_mass_integ
Ownership flag for mass_integ.
Definition: transfer.hpp:125
const CoarseFineTransformations & GetRefinementTransforms()
Definition: mesh.cpp:9928
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:267
L2Projection * F
Forward, coarse-to-fine, operator.
Definition: transfer.hpp:333
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:185
const Array< int > & GetDofMap() const
Get an Array&lt;int&gt; that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:1182
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:366
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:465
Class for integration point with weight.
Definition: intrules.hpp:25
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
bool Parallel() const
Definition: transfer.hpp:46
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
GridFunction interpolation operator applicable after mesh refinement.
Definition: fespace.hpp:298
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:957
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
Definition: sparsemat.cpp:3763
virtual void ProlongateTranspose(const Vector &x, Vector &y) const
Definition: transfer.cpp:744
const Table * GetElementToFaceOrientationTable() const
Definition: fespace.hpp:746
const FiniteElementSpace & fes_lor
Definition: transfer.hpp:187
General triple product operator x -&gt; A*B*C*x, with ownership of the factors.
Definition: operator.hpp:811
PRefinementTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator from lFESpace to hFESpace which have different FE collections.
Definition: transfer.cpp:947
const int MAX_D1D
Definition: forall.hpp:28
TrueTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator working on true degrees of freedom from lFESpace to hFESpace...
Definition: transfer.cpp:1433
IsoparametricTransformation Transf
Definition: eltrans.hpp:464
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:2781
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:1467
Lexicographic ordering for tensor-product FiniteElements.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition: fespace.cpp:1259
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:479
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
ID for class HypreParMatrix.
Definition: operator.hpp:260
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
int * GetI()
Definition: table.hpp:113
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:1553
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:750
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:70
RefCoord s[3]
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:118
Field is continuous across element interfaces.
Definition: fe_coll.hpp:45
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:1911
int RowSize(int i) const
Definition: table.hpp:108
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
virtual void InvTransformDual(double *v) const =0
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:495
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:953
L2ProjectionL2Space(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
Definition: transfer.cpp:285
virtual const Operator * GetRestrictionOperator() const
An abstract operator that performs the same action as GetRestrictionMatrix.
Definition: fespace.hpp:491
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
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:1153
virtual void TransformDual(double *v) const =0
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:449
TransferOperator(const FiniteElementSpace &lFESpace, const FiniteElementSpace &hFESpace)
Constructs a transfer operator from lFESpace to hFESpace.
Definition: transfer.cpp:902
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
void 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:1200
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:215