MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
transfer.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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 "pbilinearform.hpp"
15#include "../general/forall.hpp"
16
17namespace mfem
18{
19
21 FiniteElementSpace &ran_fes_)
22 : dom_fes(dom_fes_), ran_fes(ran_fes_),
23 oper_type(Operator::ANY_TYPE),
24 fw_t_oper(), bw_t_oper(), use_ea(false), d_mt(Device::GetHostMemoryType())
25{
26#ifdef MFEM_USE_MPI
27 const bool par_dom = dynamic_cast<ParFiniteElementSpace*>(&dom_fes);
28 const bool par_ran = dynamic_cast<ParFiniteElementSpace*>(&ran_fes);
29 MFEM_VERIFY(par_dom == par_ran, "the domain and range FE spaces must both"
30 " be either serial or parallel");
31 parallel = par_dom;
32#endif
33}
34
36 FiniteElementSpace &fes_in, FiniteElementSpace &fes_out,
37 const Operator &oper, OperatorHandle &t_oper)
38{
39 if (t_oper.Ptr())
40 {
41 return *t_oper.Ptr();
42 }
43
44 if (!Parallel())
45 {
46 const SparseMatrix *in_cP = fes_in.GetConformingProlongation();
47 const SparseMatrix *out_cR = fes_out.GetConformingRestriction();
49 {
50 const SparseMatrix *mat = dynamic_cast<const SparseMatrix *>(&oper);
51 MFEM_VERIFY(mat != NULL, "Operator is not a SparseMatrix");
52 if (!out_cR)
53 {
54 t_oper.Reset(const_cast<SparseMatrix*>(mat), false);
55 }
56 else
57 {
58 t_oper.Reset(mfem::Mult(*out_cR, *mat));
59 }
60 if (in_cP)
61 {
62 t_oper.Reset(mfem::Mult(*t_oper.As<SparseMatrix>(), *in_cP));
63 }
64 }
65 else if (oper_type == Operator::ANY_TYPE)
66 {
67 const int RP_case = bool(out_cR) + 2*bool(in_cP);
68 switch (RP_case)
69 {
70 case 0:
71 t_oper.Reset(const_cast<Operator*>(&oper), false);
72 break;
73 case 1:
74 t_oper.Reset(
75 new ProductOperator(out_cR, &oper, false, false));
76 break;
77 case 2:
78 t_oper.Reset(
79 new ProductOperator(&oper, in_cP, false, false));
80 break;
81 case 3:
82 t_oper.Reset(
84 out_cR, &oper, in_cP, false, false, false));
85 break;
86 }
87 }
88 else
89 {
90 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
91 }
92 }
93 else // Parallel() == true
94 {
95#ifdef MFEM_USE_MPI
97 {
98 const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
99 const ParFiniteElementSpace *pfes_in =
100 dynamic_cast<const ParFiniteElementSpace *>(&fes_in);
101 const ParFiniteElementSpace *pfes_out =
102 dynamic_cast<const ParFiniteElementSpace *>(&fes_out);
103 const SparseMatrix *sp_mat = dynamic_cast<const SparseMatrix *>(&oper);
104 const HypreParMatrix *hy_mat;
105 if (sp_mat)
106 {
107 SparseMatrix *RA = mfem::Mult(*out_R, *sp_mat);
108 t_oper.Reset(pfes_in->Dof_TrueDof_Matrix()->
109 LeftDiagMult(*RA, pfes_out->GetTrueDofOffsets()));
110 delete RA;
111 }
112 else if ((hy_mat = dynamic_cast<const HypreParMatrix *>(&oper)))
113 {
114 HypreParMatrix *RA =
115 hy_mat->LeftDiagMult(*out_R, pfes_out->GetTrueDofOffsets());
116 t_oper.Reset(mfem::ParMult(RA, pfes_in->Dof_TrueDof_Matrix()));
117 delete RA;
118 }
119 else
120 {
121 MFEM_ABORT("unknown Operator type");
122 }
123 }
124 else if (oper_type == Operator::ANY_TYPE)
125 {
126 const Operator *out_R = fes_out.GetRestrictionOperator();
127 t_oper.Reset(new TripleProductOperator(
128 out_R, &oper, fes_in.GetProlongationMatrix(),
129 false, false, false));
130 }
131 else
132 {
133 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
134 }
135#endif
136 }
137
138 return *t_oper.Ptr();
139}
140
141
146
148 BilinearFormIntegrator *mass_integ_, bool own_mass_integ_)
149{
150 if (own_mass_integ) { delete mass_integ; }
151
152 mass_integ = mass_integ_;
153 own_mass_integ = own_mass_integ_;
154}
155
157{
158 if (F.Ptr())
159 {
160 return *F.Ptr();
161 }
162
163 // Construct F
165 {
167 }
169 {
170 Mesh::GeometryList elem_geoms(*ran_fes.GetMesh());
171
173 for (int i = 0; i < elem_geoms.Size(); i++)
174 {
176 localP[elem_geoms[i]]);
177 }
181 }
182 else
183 {
184 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
185 }
186
187 return *F.Ptr();
188}
189
191{
192 if (B.Ptr())
193 {
194 return *B.Ptr();
195 }
196
197 // Construct B, if not set, define a suitable mass_integ
198 if (!mass_integ)
199 {
200 const FiniteElement *f_fe_0 = ran_fes.GetTypicalFE();
201 const int map_type = f_fe_0->GetMapType();
202 if (map_type == FiniteElement::VALUE ||
203 map_type == FiniteElement::INTEGRAL)
204 {
206 }
207 else if (map_type == FiniteElement::H_DIV ||
208 map_type == FiniteElement::H_CURL)
209 {
211 }
212 else
213 {
214 MFEM_ABORT("unknown type of FE space");
215 }
216 own_mass_integ = true;
217 }
219 {
222 }
223 else
224 {
225 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
226 }
227
228 return *B.Ptr();
229}
230
231
233 const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
234 MemoryType d_mt_)
235 : Operator(fes_lor_.GetVSize(), fes_ho_.GetVSize()),
236 fes_ho(fes_ho_), fes_lor(fes_lor_), d_mt(d_mt_)
237{ }
238
240 int nel_ho, int nel_lor, const CoarseFineTransformations& cf_tr)
241{
242 // Construct the mapping from HO to LOR
243 // ho2lor.GetRow(iho) will give all the LOR elements contained in iho
244 ho2lor.MakeI(nel_ho);
245 for (int ilor = 0; ilor < nel_lor; ++ilor)
246 {
247 int iho = cf_tr.embeddings[ilor].parent;
248 ho2lor.AddAColumnInRow(iho);
249 }
250 ho2lor.MakeJ();
251 for (int ilor = 0; ilor < nel_lor; ++ilor)
252 {
253 int iho = cf_tr.embeddings[ilor].parent;
254 ho2lor.AddConnection(iho, ilor);
255 }
256 ho2lor.ShiftUpI();
257}
258
260 Geometry::Type geom, const FiniteElement& fe_ho,
261 const FiniteElement& fe_lor, ElementTransformation* tr_ho,
262 ElementTransformation* tr_lor,
264 DenseMatrix& M_mixed_el) const
265{
266 int order = fe_lor.GetOrder() + fe_ho.GetOrder() + tr_lor->OrderW();
267 const IntegrationRule* ir = &IntRules.Get(geom, order);
268 M_mixed_el = 0.0;
269 for (int i = 0; i < ir->GetNPoints(); i++)
270 {
271 const IntegrationPoint& ip_lor = ir->IntPoint(i);
272 IntegrationPoint ip_ho;
273 ip_tr.Transform(ip_lor, ip_ho);
274 Vector shape_lor(fe_lor.GetDof());
275 fe_lor.CalcShape(ip_lor, shape_lor);
276 Vector shape_ho(fe_ho.GetDof());
277 tr_ho->SetIntPoint(&ip_ho);
278 fe_ho.CalcPhysShape(*tr_ho, shape_ho);
279 tr_lor->SetIntPoint(&ip_lor);
280 // For now we use the geometry information from the LOR space, which means
281 // we won't be mass conservative if the mesh is curved
282 real_t w = ip_lor.weight;
283 if (fe_lor.GetMapType() == FiniteElement::VALUE)
284 {
285 w *= tr_lor->Weight();
286 }
287 shape_lor *= w;
288 AddMultVWt(shape_lor, shape_ho, M_mixed_el);
289 }
290}
291
293 Geometry::Type geom, const FiniteElement& fe_ho,
294 const FiniteElement& fe_lor, ElementTransformation* el_tr,
296 DenseMatrix& B_L, DenseMatrix& B_H) const
297{
298 int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
299 const IntegrationRule* ir = &IntRules.Get(geom, order);
300
301 for (int i = 0; i < ir->GetNPoints(); i++)
302 {
303 const IntegrationPoint& ip_lor = ir->IntPoint(i);
304 IntegrationPoint ip_ho;
305
306 // maps integration point ip_lor -> ip_ho
307 ip_tr.Transform(ip_lor, ip_ho);
308 Vector shape_lor(fe_lor.GetDof());
309 fe_lor.CalcShape(ip_lor, shape_lor);
310 Vector shape_ho(fe_ho.GetDof());
311 fe_ho.CalcShape(ip_ho, shape_ho);
312
313 for (int j=0; j<shape_lor.Size(); ++j)
314 {
315 B_L(i, j) = shape_lor(j);
316 }
317
318 for (int j=0; j<shape_ho.Size(); ++j)
319 {
320 B_H(i, j) = shape_ho(j);
321 }
322 }
323
324}
325
327 const FiniteElementSpace& fes_ho_ea,
328 const FiniteElementSpace& fes_lor_ea,
329 Vector &M_LH, MemoryType d_mt_)
330{
331 Mesh* mesh_ho = fes_ho_ea.GetMesh();
332 Mesh* mesh_lor = fes_lor_ea.GetMesh();
333 int nel_ho = mesh_ho->GetNE();
334 int nel_lor = mesh_lor->GetNE();
335
336 const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
337
338 int nref_max = 0;
340 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
341 for (int ig = 0; ig < geoms.Size(); ++ig)
342 {
343 Geometry::Type geom = geoms[ig];
344 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
345 }
346
347 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
348
350 IsoparametricTransformation &emb_tr = ip_tr.Transf;
351
352 // Gather basis functions (B_L, B_HO) and data at quadrature points
353 DenseTensor B_L, B_H, D;
354 {
355 // Assume all HO elements are LOR in the same way
356 const int iho = 0;
357 {
358 Array<int> lor_els;
359 ho2lor.GetRow(iho, lor_els);
360 int nref = ho2lor.RowSize(iho);
361
362 Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
363 const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
364 const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
365
366 // Allocate space for DenseTensors
367 ElementTransformation *el_tr = fes_lor_ea.GetElementTransformation(0);
368 int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
369 const IntegrationRule* ir_ea = &IntRules.Get(geom, order);
370 int qPts = ir_ea->GetNPoints();
371
372 // Containers for the basis functions sampled at quadrature points
373 B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt);
374 B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt);
375 D.SetSize(qPts, nref, nel_ho, d_mt);
376
377 const GeometricFactors *geo_facts =
379
380 MFEM_ASSERT(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
381
382 // Setup data at quadrature points
383 // TODO add support for user coefficient
384 const auto W = Reshape(ir_ea->GetWeights().Read(), qPts);
385 const auto J = Reshape(geo_facts->detJ.Read(), qPts, nel_lor);
386 const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
387
388 mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
389 {
390 const int q = tid % qPts;
391 const int iref = (tid / qPts) % nref;
392 const int iho = (tid / (qPts * nref)) % nel_ho;
393
394 const int lo_el_id = iref + nref*iho;
395 const real_t detJ = J(q, lo_el_id);
396
397 d_D(q, iref, iho) = W(q) * detJ;
398
399 });
400
401 emb_tr.SetIdentityTransformation(geom);
402 const DenseTensor &pmats = cf_tr.point_matrices[geom];
403
404 // Collect the basis functions
405 for (int iref = 0; iref < nref; ++iref)
406 {
407 int ilor = lor_els[iref];
408 // Now assemble the block-row of the mixed mass matrix associated
409 // with integrating HO functions against LOR functions on the LOR
410 // sub-element.
411
412 // Create the transformation that embeds the fine low-order element
413 // within the coarse high-order element in reference space
414 emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
415
416 DenseMatrix &b_lo = B_L(ilor);
417 DenseMatrix &b_ho = B_H(ilor);
418
419 ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, b_lo, b_ho);
420
421 } // loop over subcells of ho element
422 // end of quadrature point setup
423 }
424
425 } // completed setup of basis function and quadrature point
426
427 // Assemble mixed mass matrix
428 {
429 int iho = 0;
430 Array<int> lor_els;
431 ho2lor.GetRow(iho, lor_els);
432 int nref = ho2lor.RowSize(iho);
433
434 const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
435 const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
436 const int ndof_ho = fe_ho.GetDof();
437 const int ndof_lor = fe_lor.GetDof();
438
439 const int qPts = D.SizeI();
440
441 M_LH.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt);
442
443 // Rows x columns
444 // Recall MFEM is column major
445 // rows x columns is inverted - matrix is ndof_lor x ndof_ho
446 auto v_M_LH = Reshape(M_LH.Write(), ndof_lor, ndof_ho, nref,
447 nel_ho);
448
449 const int fe_ho_ndof = fe_ho.GetDof();
450 const int fe_lor_ndof = fe_lor.GetDof();
451
452 auto d_B_L = Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
453 auto d_B_H = Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
454 auto d_D = Reshape(D.Read(), qPts, nref, nel_ho);
455
456 mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (int idx)
457 {
458 const int bh = idx % fe_ho_ndof;
459 const int iref = (idx / fe_ho_ndof) % nref;
460 const int iho = idx / fe_ho_ndof / nref;
461 // (B_lo_dofs x Q) x (Q x B_ho_dofs)
462 for (int bl = 0; bl < fe_lor_ndof; ++bl)
463 {
464 real_t dot = 0.0;
465 for (int qi=0; qi<qPts; ++qi)
466 {
467 dot += d_B_L(qi, bl, iref) * d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
468 }
469 // column major storage
470 v_M_LH(bl, bh, iref, iho) = dot;
471 }
472 });
473 } // end of mixed assembly mass matrix
474}
475
477(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
478 const bool use_ea_, MemoryType d_mt_)
479 : L2Projection(fes_ho_, fes_lor_, d_mt_),
480 use_ea(use_ea_)
481{
482 if (use_ea)
483 {
485 return;
486 }
487
488 Mesh *mesh_ho = fes_ho.GetMesh();
489 Mesh *mesh_lor = fes_lor.GetMesh();
490 int nel_ho = mesh_ho->GetNE();
491 int nel_lor = mesh_lor->GetNE();
492
493 // The prolongation operation is only well-defined when the LOR space has at
494 // least as many DOFs as the high-order space.
495 const bool build_P = fes_lor.GetTrueVSize() >= fes_ho.GetTrueVSize();
496
497 // If the local mesh is empty, skip all computations
498 if (nel_ho == 0) { return; }
499
500 const CoarseFineTransformations &cf_tr = mesh_lor->GetRefinementTransforms();
501
502 int nref_max = 0;
504 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
505 for (int ig = 0; ig < geoms.Size(); ++ig)
506 {
507 Geometry::Type geom = geoms[ig];
508 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
509 }
510
511 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
512
513 offsets.SetSize(nel_ho+1);
514 offsets[0] = 0;
515 for (int iho = 0; iho < nel_ho; ++iho)
516 {
517 int nref = ho2lor.RowSize(iho);
518 const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
519 const FiniteElement &fe_lor = *fes_lor.GetFE(ho2lor.GetRow(iho)[0]);
520 offsets[iho+1] = offsets[iho] + fe_ho.GetDof()*fe_lor.GetDof()*nref;
521 }
522 // R will contain the restriction (L^2 projection operator) defined on each
523 // coarse HO element (and corresponding patch of LOR elements)
524 R.SetSize(offsets[nel_ho]);
525 if (build_P)
526 {
527 // P will contain the corresponding prolongation operator
528 P.SetSize(offsets[nel_ho]);
529 }
530
532 IsoparametricTransformation &emb_tr = ip_tr.Transf;
533
534 for (int iho = 0; iho < nel_ho; ++iho)
535 {
536 Array<int> lor_els;
537 ho2lor.GetRow(iho, lor_els);
538 int nref = ho2lor.RowSize(iho);
539
540 Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
541 const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
542 const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
543 int ndof_ho = fe_ho.GetDof();
544 int ndof_lor = fe_lor.GetDof();
545
547
548 emb_tr.SetIdentityTransformation(geom);
549 const DenseTensor &pmats = cf_tr.point_matrices[geom];
550
551 DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
552
553 DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
554 DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
555
557 DenseMatrix M_lor_el(ndof_lor, ndof_lor);
558 DenseMatrixInverse Minv_lor_el(&M_lor_el);
559 DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
560 DenseMatrix M_mixed_el(ndof_lor, ndof_ho);
561
562 Minv_lor = 0.0;
563 M_lor = 0.0;
564
565 DenseMatrix RtMlor(ndof_ho, ndof_lor*nref);
566 DenseMatrix RtMlorR(ndof_ho, ndof_ho);
567 DenseMatrixInverse RtMlorR_inv(&RtMlorR);
568
569 for (int iref = 0; iref < nref; ++iref)
570 {
571 // Assemble the low-order refined mass matrix and invert locally
572 int ilor = lor_els[iref];
574 mi.AssembleElementMatrix(fe_lor, *tr_lor, M_lor_el);
575 M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
576 Minv_lor_el.Factor();
577 Minv_lor_el.GetInverseMatrix(M_lor_el);
578 // Insert into the diagonal of the patch LOR mass matrix
579 Minv_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
580
581 // Now assemble the block-row of the mixed mass matrix associated
582 // with integrating HO functions against LOR functions on the LOR
583 // sub-element.
584
585 // Create the transformation that embeds the fine low-order element
586 // within the coarse high-order element in reference space
587 emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
588
589 ElemMixedMass(geom, fe_ho, fe_lor, tr_ho, tr_lor, ip_tr, M_mixed_el);
590
591 M_mixed.CopyMN(M_mixed_el, iref*ndof_lor, 0);
592 }
593 mfem::Mult(Minv_lor, M_mixed, R_iho);
594
595 if (build_P)
596 {
597 DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
598
599 mfem::MultAtB(R_iho, M_lor, RtMlor);
600 mfem::Mult(RtMlor, R_iho, RtMlorR);
601 RtMlorR_inv.Factor();
602 RtMlorR_inv.Mult(RtMlor, P_iho);
603 }
604 }
605
606}
607
608
610{
611 Mesh *mesh_ho = fes_ho.GetMesh();
612 Mesh *mesh_lor = fes_lor.GetMesh();
613 int nel_ho = mesh_ho->GetNE();
614 int nel_lor = mesh_lor->GetNE();
615
616 // The prolongation operation is only well-defined when the LOR space has at
617 // least as many DOFs as the high-order space.
618 const bool build_P = fes_lor.GetTrueVSize() >= fes_ho.GetTrueVSize();
619
620 // If the local mesh is empty, skip all computations
621 if (nel_ho == 0) { return; }
622
623 const CoarseFineTransformations &cf_tr = mesh_lor->GetRefinementTransforms();
624
625 int nref_max = 0;
627 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
628 for (int ig = 0; ig < geoms.Size(); ++ig)
629 {
630 Geometry::Type geom = geoms[ig];
631 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
632 }
633
634 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
635
636 offsets.SetSize(nel_ho+1);
637 offsets[0] = 0;
638 for (int iho = 0; iho < nel_ho; ++iho)
639 {
640 int nref = ho2lor.RowSize(iho);
641 const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
642 const FiniteElement &fe_lor = *fes_lor.GetFE(ho2lor.GetRow(iho)[0]);
643 offsets[iho+1] = offsets[iho] + fe_ho.GetDof()*fe_lor.GetDof()*nref;
644 }
645
646 // R will contain the restriction (L^2 projection operator) defined on each
647 // coarse HO element (and corresponding patch of LOR elements)
648 R.SetSize(offsets[nel_ho]);
649
650 if (build_P)
651 {
652 // P will contain the corresponding prolongation operator
653 P.SetSize(offsets[nel_ho]);
654 }
655
656 // Assemble mixed mass matrix
657 Vector M_mixed_all;
658 MixedMassEA(fes_ho, fes_lor, M_mixed_all, d_mt);
659
660
661 // R = inv(M_L) * M_mixed
662 // Need to compute M_L
663 // Note: Using user-inputted M_LH IntegrationRule ir
664 // (higher order than needed) in order to re-use coeff
666
667 Vector M_ea_lor;
668 int ndof_lor;
669 int ndof_ho;
670 int nref;
671 {
672 int iho = 0;
673 Array<int> lor_els;
674 ho2lor.GetRow(iho, lor_els);
675 nref = ho2lor.RowSize(iho);
676
677 const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
678 const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
679 ndof_ho = fe_ho.GetDof();
680 ndof_lor = fe_lor.GetDof();
681
682 M_ea_lor.SetSize(ndof_lor*ndof_lor*nel_lor, d_mt);
683 }
684
685 const bool add = false;
686 mi.AssembleEA(fes_lor, M_ea_lor, add);
687
688 DenseTensor Minv_ear_lor;
689 Minv_ear_lor.SetSize(ndof_lor, ndof_lor, nel_lor, d_mt);
690 Minv_ear_lor.GetMemory().CopyFrom(M_ea_lor.GetMemory(), M_ea_lor.Size());
691
692 BatchedLinAlg::Invert(Minv_ear_lor);
693 {
694 // Recall mfem is column major
695 // ndof_lor x ndof_ho
696 auto v_M_mixed_all = Reshape(M_mixed_all.Read(), ndof_lor, ndof_ho, nref,
697 nel_ho);
698
699 // matrix is symmetric
700 auto v_Minv_ear_lor = Reshape(Minv_ear_lor.Read(), ndof_lor, ndof_lor,
701 nel_lor);
702
703 // ndof_lor x ndof_ho
704 auto v_R = Reshape(R.Write(), ndof_lor, nref, ndof_ho, nel_ho);
705
706 MFEM_VERIFY(nel_lor==nel_ho*nref, "nel_lor != nel_ho*nref");
707
708 // (ndofs_lor x ndofs_lor) x (ndofs_lor x ndof_ho)
709 mfem::forall(ndof_lor * nref * ndof_ho * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
710 {
711
712 const int i = tid % ndof_lor;
713 const int iref = (tid / ndof_lor) % nref;
714 const int j = (tid / (ndof_lor * nref) ) % ndof_ho;
715 const int iho = (tid / (ndof_lor * nref * ndof_ho)) % nel_ho;
716
717 const int lor_idx = iref + iho * nref;
718
719 //matrices are stored in the transpose position
720 real_t dot = 0.0;
721 for (int k=0; k<ndof_lor; ++k)
722 {
723 dot += v_Minv_ear_lor(i, k, lor_idx) * v_M_mixed_all(k, j, iref, iho);
724 }
725 v_R(i, iref, j, iho) = dot;
726
727 });
728 }
729
730 if (build_P)
731 {
732 // P = inv(R^T M_L R) * R^T M_L
733
734 // M_lor is size of ndof_lor x ndof_lor
735 // R is size of (ndof_lor x nref x ndof_ho)
736 auto v_M_ea_lor = Reshape(M_ea_lor.Read(), ndof_lor, ndof_lor, nel_lor);
737 auto v_R = Reshape(R.Read(), ndof_lor, nref, ndof_ho, nel_ho);
738 // R^T M_LO is of size nref x ndof_lor
739
740 // Compute R^T M_L
741 Vector RtM_L(ndof_ho*nref*ndof_lor*nel_ho, d_mt);
742 auto v_RtM_L = Reshape(RtM_L.Write(), ndof_ho, ndof_lor, nref, nel_ho);
743
744 mfem::forall(ndof_lor * nref * ndof_ho * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
745 {
746
747 const int jlo = tid % ndof_lor;
748 const int iref = (tid / ndof_lor) % nref;
749 const int iho = (tid / (ndof_lor * nref)) % ndof_ho;
750 const int e = (tid / (ndof_lor * nref * ndof_ho)) % nel_ho;
751
752 const int lor_idx = iref + e * nref;
753
754 real_t dot = 0.0;
755 for (int t=0; t<ndof_lor; ++t)
756 {
757 dot += v_R(t, iref, iho, e) * v_M_ea_lor(t, jlo, lor_idx);
758 }
759
760 v_RtM_L(iho, jlo, iref, e) = dot;
761
762 });
763
764 // Resulting matrix should be: ndof_ho x ndof_ho
765 // R^T M_L x R
766 DenseTensor RtM_L_dt;
767 RtM_L_dt.NewMemoryAndSize(RtM_L.GetMemory(), ndof_ho, ndof_lor*nref,
768 nel_ho, false);
769 Vector R_vec;
770 R_vec.NewMemoryAndSize(R.GetMemory(), R.Size(), false);
771 Vector RtM_LR(ndof_ho * ndof_ho * nel_ho, d_mt);
772 BatchedLinAlg::Mult(RtM_L_dt, R_vec, RtM_LR);
773 // Ensure that changes to the alias R_vec are propagated to the base, R
774 R_vec.GetMemory().SyncAlias(R.GetMemory(), P.Size());
775
776 // Compute the inverse of InvRtM_LR
777 DenseTensor InvRtM_LR;
778 InvRtM_LR.NewMemoryAndSize(RtM_LR.GetMemory(), ndof_ho, ndof_ho, nel_ho, false);
779 BatchedLinAlg::Invert(InvRtM_LR);
780
781 // Form P
782 // P should be of dimension (ndof_ho x ndof_ho) x (ndof_ho x nref*ndof_lor)
783 // P ndof_ho x nref*ndof_lor
784 Vector P_vec;
785 P_vec.NewMemoryAndSize(P.GetMemory(), P.Size(), false);
786 BatchedLinAlg::Mult(InvRtM_LR, RtM_L, P_vec);
787 // Ensure that changes to the alias P_vec are propagated to the base, P
788 P_vec.GetMemory().SyncAlias(P.GetMemory(), P.Size());
789 }
790}
791
793 const Vector &x, Vector &y) const
794{
795
796 if (use_ea)
797 {
798 return EAMult(x,y);
799 }
800
801
802 int vdim = fes_ho.GetVDim();
803 Array<int> vdofs;
804 DenseMatrix xel_mat, yel_mat;
805 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
806 {
807 int nref = ho2lor.RowSize(iho);
808 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
809 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
810 xel_mat.SetSize(ndof_ho, vdim);
811 yel_mat.SetSize(ndof_lor*nref, vdim);
812 DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
813
814 fes_ho.GetElementVDofs(iho, vdofs);
815 x.GetSubVector(vdofs, xel_mat.GetData());
816 mfem::Mult(R_iho, xel_mat, yel_mat);
817 // Place result correctly into the low-order vector
818 for (int iref = 0; iref < nref; ++iref)
819 {
820 int ilor = ho2lor.GetRow(iho)[iref];
821 for (int vd=0; vd<vdim; ++vd)
822 {
823 fes_lor.GetElementDofs(ilor, vdofs);
824 fes_lor.DofsToVDofs(vd, vdofs);
825 y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
826 }
827 }
828 }
829}
830
832 const Vector &x, Vector &y) const
833{
834 const int iho = 0;
835 const int nref = ho2lor.RowSize(iho);
836 const int ndof_ho = fes_ho.GetFE(iho)->GetDof();
837 const int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
838 const int nel_ho = fes_ho.GetMesh()->GetNE();
839
840 DenseTensor R_dt;
841 R_dt.NewMemoryAndSize(R.GetMemory(), ndof_lor*nref, ndof_ho, nel_ho, false);
842 BatchedLinAlg::Mult(R_dt, x, y);
843}
844
846 const Vector &x, Vector &y) const
847{
848
849 if (use_ea)
850 {
851 return EAMultTranspose(x,y);
852 }
853
854 int vdim = fes_ho.GetVDim();
855 Array<int> vdofs;
856 DenseMatrix xel_mat, yel_mat;
857 y = 0.0;
858 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
859 {
860 int nref = ho2lor.RowSize(iho);
861 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
862 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
863 xel_mat.SetSize(ndof_lor*nref, vdim);
864 yel_mat.SetSize(ndof_ho, vdim);
865 DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
866
867 // Extract the LOR DOFs
868 for (int iref=0; iref<nref; ++iref)
869 {
870 int ilor = ho2lor.GetRow(iho)[iref];
871 for (int vd=0; vd<vdim; ++vd)
872 {
873 fes_lor.GetElementDofs(ilor, vdofs);
874 fes_lor.DofsToVDofs(vd, vdofs);
875 x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
876 }
877 }
878 // Multiply locally by the transpose
879 mfem::MultAtB(R_iho, xel_mat, yel_mat);
880 // Place the result in the HO vector
881 fes_ho.GetElementVDofs(iho, vdofs);
882 y.AddElementVector(vdofs, yel_mat.GetData());
883 }
884
885}
886
888 const Vector &x, Vector &y) const
889{
890 const int iho = 0;
891 const int nref = ho2lor.RowSize(iho);
892 const int ndof_ho = fes_ho.GetFE(iho)->GetDof();
893 const int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
894 const int nel_ho = fes_ho.GetMesh()->GetNE();
895
896 DenseTensor R_dt;
897 R_dt.NewMemoryAndSize(R.GetMemory(), ndof_lor*nref, ndof_ho, nel_ho, false);
899}
900
902 const Vector &x, Vector &y) const
903{
904
905 if (fes_ho.GetNE() == 0) { return; }
906
907 if (use_ea)
908 {
909 return EAProlongate(x,y);
910 }
911
912 MFEM_VERIFY(P.Size() > 0, "Prolongation not supported for these spaces.")
913 int vdim = fes_ho.GetVDim();
914 Array<int> vdofs;
915 DenseMatrix xel_mat,yel_mat;
916 y = 0.0;
917 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
918 {
919 int nref = ho2lor.RowSize(iho);
920 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
921 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
922 xel_mat.SetSize(ndof_lor*nref, vdim);
923 yel_mat.SetSize(ndof_ho, vdim);
924 DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
925
926 // Extract the LOR DOFs
927 for (int iref = 0; iref < nref; ++iref)
928 {
929 int ilor = ho2lor.GetRow(iho)[iref];
930 for (int vd = 0; vd < vdim; ++vd)
931 {
932 fes_lor.GetElementDofs(ilor, vdofs);
933 fes_lor.DofsToVDofs(vd, vdofs);
934 x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
935 }
936 }
937 // Locally prolongate
938 mfem::Mult(P_iho, xel_mat, yel_mat);
939 // Place the result in the HO vector
940 fes_ho.GetElementVDofs(iho, vdofs);
941 y.AddElementVector(vdofs, yel_mat.GetData());
942 }
943
944}
945
947 const Vector &x, Vector &y) const
948{
949 const int iho = 0;
950 const int nref = ho2lor.RowSize(iho);
951 const int ndof_ho = fes_ho.GetFE(iho)->GetDof();
952 const int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
953 const int nel_ho = fes_ho.GetMesh()->GetNE();
954
955 DenseTensor P_dt;
956 P_dt.NewMemoryAndSize(P.GetMemory(), ndof_ho, ndof_lor * nref, nel_ho, false);
957 BatchedLinAlg::Mult(P_dt, x, y);
958}
959
961 const Vector &x, Vector &y) const
962{
963
964 if (use_ea)
965 {
966 return EAProlongateTranspose(x,y);
967 }
968
969
970 if (fes_ho.GetNE() == 0) { return; }
971 MFEM_VERIFY(P.Size() > 0, "Prolongation not supported for these spaces.")
972 int vdim = fes_ho.GetVDim();
973 Array<int> vdofs;
974 DenseMatrix xel_mat,yel_mat;
975 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
976 {
977 int nref = ho2lor.RowSize(iho);
978 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
979 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
980 xel_mat.SetSize(ndof_ho, vdim);
981 yel_mat.SetSize(ndof_lor*nref, vdim);
982 DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
983
984 fes_ho.GetElementVDofs(iho, vdofs);
985 x.GetSubVector(vdofs, xel_mat.GetData());
986 mfem::MultAtB(P_iho, xel_mat, yel_mat);
987
988 // Place result correctly into the low-order vector
989 for (int iref = 0; iref < nref; ++iref)
990 {
991 int ilor = ho2lor.GetRow(iho)[iref];
992 for (int vd=0; vd<vdim; ++vd)
993 {
994 fes_lor.GetElementDofs(ilor, vdofs);
995 fes_lor.DofsToVDofs(vd, vdofs);
996 y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
997 }
998 }
999 }
1000
1001}
1002
1004 const Vector &x, Vector &y) const
1005{
1006 const int iho = 0;
1007 const int nref = ho2lor.RowSize(iho);
1008 const int ndof_ho = fes_ho.GetFE(iho)->GetDof();
1009 const int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
1010 const int nel_ho = fes_ho.GetMesh()->GetNE();
1011
1012 DenseTensor P_dt;
1013 P_dt.NewMemoryAndSize(P.GetMemory(), ndof_ho, ndof_lor * nref, nel_ho, false);
1014 BatchedLinAlg::MultTranspose(P_dt, x, y);
1015}
1016
1018 const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_,
1019 const bool use_ea_, MemoryType d_mt_)
1020 : L2Projection(fes_ho_, fes_lor_, d_mt_),
1021 use_ea(use_ea_)
1022{
1023
1024 // need scalar to keep dimensions matching (operators are built to apply
1025 // individually on each vdim)
1026 // needed in both matrix and element based versions
1028 fes_ho.FEColl(), 1));
1030 fes_lor.FEColl(), 1));
1031
1032 if (use_ea)
1033 {
1035 return;
1036 }
1037
1038 std::unique_ptr<SparseMatrix> R_mat, M_LH_mat;
1039
1040 std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH();
1041
1042 const SparseMatrix *P_ho = fes_ho_scalar->GetConformingProlongation();
1043 const SparseMatrix *P_lor = fes_lor_scalar->GetConformingProlongation();
1044
1045 if (P_ho || P_lor)
1046 {
1047 if (P_ho && P_lor)
1048 {
1049 R_mat.reset(RAP(*P_lor, *R_mat, *P_ho));
1050 M_LH_mat.reset(RAP(*P_lor, *M_LH_mat, *P_ho));
1051 }
1052 else if (P_ho)
1053 {
1054 R_mat.reset(mfem::Mult(*R_mat, *P_ho));
1055 M_LH_mat.reset(mfem::Mult(*M_LH_mat, *P_ho));
1056 }
1057 else // P_lor != nullptr
1058 {
1059 R_mat.reset(mfem::Mult(*P_lor, *R_mat));
1060 M_LH_mat.reset(mfem::Mult(*P_lor, *M_LH_mat));
1061 }
1062 }
1063
1064 SparseMatrix *RTxM_LH_mat = TransposeMult(*R_mat, *M_LH_mat);
1065 precon.reset(new DSmoother(*RTxM_LH_mat));
1066
1067 // Set ownership
1068 RTxM_LH.reset(RTxM_LH_mat);
1069 R = std::move(R_mat);
1070 M_LH = std::move(M_LH_mat);
1071
1072 SetupPCG();
1073}
1074
1075#ifdef MFEM_USE_MPI
1076
1078 const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor,
1079 const bool use_ea_, MemoryType d_mt_)
1080 : L2Projection(pfes_ho, pfes_lor, d_mt_),
1081 use_ea(use_ea_), pcg(pfes_ho.GetComm())
1082{
1083
1084 // need scalar to keep dimensions matching (operators are built to apply
1085 // individually on each vdim)
1086 // needed in both matrix and element based versions
1088 pfes_ho.FEColl(), 1));
1090 pfes_lor.FEColl(), 1));
1091
1092 if (use_ea)
1093 {
1094 EAL2ProjectionH1Space(pfes_ho, pfes_lor);
1095 return;
1096 }
1097
1098 std::tie(R, M_LH) = ComputeSparseRAndM_LH();
1099
1100
1101 HypreParMatrix R_local = HypreParMatrix(pfes_ho.GetComm(),
1102 pfes_lor_scalar->GlobalVSize(),
1103 pfes_ho_scalar->GlobalVSize(),
1104 pfes_lor_scalar->GetDofOffsets(),
1105 pfes_ho_scalar->GetDofOffsets(),
1106 static_cast<SparseMatrix*>(R.get()));
1107 HypreParMatrix M_LH_local = HypreParMatrix(pfes_ho.GetComm(),
1108 pfes_lor_scalar->GlobalVSize(),
1109 pfes_ho_scalar->GlobalVSize(),
1110 pfes_lor_scalar->GetDofOffsets(),
1111 pfes_ho_scalar->GetDofOffsets(),
1112 static_cast<SparseMatrix*>(M_LH.get()));
1113
1114 HypreParMatrix *R_mat = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(),
1115 &R_local, pfes_ho_scalar->Dof_TrueDof_Matrix());
1116 HypreParMatrix *M_LH_mat = RAP(pfes_lor_scalar->Dof_TrueDof_Matrix(),
1117 &M_LH_local, pfes_ho_scalar->Dof_TrueDof_Matrix());
1118
1119 std::unique_ptr<HypreParMatrix> R_T(R_mat->Transpose());
1120 HypreParMatrix *RTxM_LH_mat = ParMult(R_T.get(), M_LH_mat, true);
1121
1122 HypreBoomerAMG *amg = new HypreBoomerAMG(*RTxM_LH_mat);
1123 amg->SetPrintLevel(0);
1124
1125 R.reset(R_mat);
1126 M_LH.reset(M_LH_mat);
1127 RTxM_LH.reset(RTxM_LH_mat);
1128 precon.reset(amg);
1129
1130 SetupPCG();
1133}
1134
1135#endif
1136
1138{
1139 // Basic PCG solver setup
1140 pcg.SetPrintLevel(0);
1141 // pcg.SetPrintLevel(IterativeSolver::PrintLevel().Summary());
1142 pcg.SetMaxIter(1000);
1143 // initial values for relative and absolute tolerance
1144 pcg.SetRelTol(1e-13);
1145 pcg.SetAbsTol(1e-13);
1146 pcg.SetPreconditioner(*precon);
1147 pcg.SetOperator(*RTxM_LH);
1148}
1149
1151{
1152 Mesh* mesh_ho = fes_ho.GetMesh();
1153 Mesh* mesh_lor = fes_lor.GetMesh();
1154 int nel_ho = mesh_ho->GetNE();
1155 int nel_lor = mesh_lor->GetNE();
1156 int ndof_ho = fes_ho.GetNDofs();
1157 int ndof_lor = fes_lor.GetNDofs();
1158
1159 // If the local mesh is empty, skip all computations
1160 if (nel_ho == 0)
1161 {
1162 return;
1163 }
1164
1165 const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
1166
1167 int nref_max = 0;
1169 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
1170 for (int ig = 0; ig < geoms.Size(); ++ig)
1171 {
1172 Geometry::Type geom = geoms[ig];
1173 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
1174 }
1175
1176 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
1177
1178 // lumped M_H and inv lumped M_L
1179
1180 // M_H contains the lumped (row sum) high order mass matrix. This is built for
1181 // preconditioning the inverse needed to build the prolongation operator P
1182 Vector M_H(ndof_ho);
1183 M_H = 0.0;
1184 // ML_inv_ea contains the inverse lumped (row sum) mass matrix. Note that the
1185 // method will also work with a full (consistent) mass matrix, though this is
1186 // not implemented here. L refers to the low-order refined mesh
1187 ML_inv_ea.SetSize(ndof_lor);
1188 ML_inv_ea = 0.0;
1189
1190 BilinearForm Mho(fes_ho_scalar.get());
1193 Mho.Assemble();
1194
1195 // Processor local lumped Mass
1196 Vector ones_ho(Mho.Width()); ones_ho = 1.0;
1197 M_H = 0.0;
1198 Mho.Mult(ones_ho, M_H);
1199
1200 BilinearForm Mlor(fes_lor_scalar.get());
1203 Mlor.Assemble();
1204
1205 Vector ones_lor(Mlor.Width()); ones_lor = 1.0;
1206 Mlor.Mult(ones_lor, ML_inv_ea);
1207
1208 // DOF by DOF inverse of non-zero entries
1209 LumpedMassInverse(ML_inv_ea);
1210
1211 // mixed mass M_LH
1212 MixedMassEA(fes_ho, fes_lor, M_LH_ea, d_mt);
1213
1214 // Set ownership
1215 M_LH_local_op = new H1SpaceMixedMassOperator(fes_ho_scalar.get(),
1216 fes_lor_scalar.get(),
1217 &ho2lor,
1218 &M_LH_ea);
1219
1220 ML_inv_vea.reset(new H1SpaceLumpedMassOperator(fes_ho_scalar.get(),
1221 fes_lor_scalar.get(),
1222 ML_inv_ea));
1223 M_LH.reset(M_LH_local_op);
1224 R.reset(new ProductOperator(ML_inv_vea.get(), M_LH.get(), false,
1225 false));
1226
1227 Array<int> ess_tdof_list; // leave empty
1228 precon.reset(new OperatorJacobiSmoother(M_H, ess_tdof_list));
1229
1230 TransposeOperator* RT = new TransposeOperator(R.get());
1231 RTxM_LH.reset(new ProductOperator(RT, M_LH.get(), true, false));
1232
1233 SetupPCG();
1234}
1235
1236#ifdef MFEM_USE_MPI
1238(const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor)
1239{
1240 Mesh* mesh_ho = pfes_ho.GetParMesh();
1241 Mesh* mesh_lor = pfes_lor.GetParMesh();
1242 int nel_ho = mesh_ho->GetNE();
1243 int nel_lor = mesh_lor->GetNE();
1244 int ndof_ho = pfes_ho.GetNDofs();
1245 int ndof_lor = pfes_lor.GetNDofs();
1246
1247
1248 // If the local mesh is empty, skip all computations
1249 if (nel_ho == 0)
1250 {
1251 return;
1252 }
1253
1254 const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
1255
1256 int nref_max = 0;
1258 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
1259 for (int ig = 0; ig < geoms.Size(); ++ig)
1260 {
1261 Geometry::Type geom = geoms[ig];
1262 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
1263 }
1264
1265 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
1266
1267 // lumped M_H and inv lumped M_L
1268
1269 // M_H contains the lumped (row sum) high order mass matrix. This is built for
1270 // preconditioning the inverse needed to build the prolongation operator P
1271 Vector M_H(ndof_ho);
1272 M_H = 0.0;
1273 // ML_inv_ea contains the inverse lumped (row sum) mass matrix. Note that the
1274 // method will also work with a full (consistent) mass matrix, though this is
1275 // not implemented here. L refers to the low-order refined mesh
1276 ML_inv_ea.SetSize(ndof_lor);
1277 ML_inv_ea = 0.0;
1278
1279 ParBilinearForm pMho(pfes_ho_scalar.get());
1282 pMho.Assemble();
1283
1284 // Processor local lumped Mass
1285 Vector ones_ho(pMho.Width()); ones_ho = 1.0;
1286 M_H = 0.0;
1287 pMho.Mult(ones_ho, M_H);
1288
1289 ParBilinearForm pMlor(pfes_lor_scalar.get());
1292 pMlor.Assemble();
1293
1294 Vector ones_lor(pMlor.Width()); ones_lor = 1.0;
1295 pMlor.Mult(ones_lor, ML_inv_ea);
1296
1297
1298 // DOF by DOF inverse of non-zero entries
1299 LumpedMassInverse(ML_inv_ea);
1300
1301 // mixed mass M_LH
1302 MixedMassEA(*pfes_ho_scalar.get(), *pfes_lor_scalar.get(), M_LH_ea, d_mt);
1303
1304 // Set ownership
1305 M_LH_local_op = new H1SpaceMixedMassOperator(pfes_ho_scalar.get(),
1306 pfes_lor_scalar.get(),
1307 &ho2lor, &M_LH_ea);
1308
1309 const Operator *P_ho = pfes_ho_scalar->GetProlongationMatrix();
1310 const Operator *P_lor = pfes_lor_scalar->GetProlongationMatrix();
1311
1312 Array<int> ess_tdof_list; // leave empty
1313
1314 if (P_ho || P_lor)
1315 {
1316 if (P_ho && P_lor)
1317 {
1318 Operator *Pt_lor = new TransposeOperator(P_lor);
1319 RML_inv.SetSize(pfes_lor_scalar->GetTrueVSize());
1320 GetTDofs(*pfes_lor_scalar, ML_inv_ea, RML_inv);
1321 ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
1322 pfes_lor_scalar.get(),
1323 RML_inv));
1324 M_LH.reset(new TripleProductOperator(Pt_lor, M_LH_local_op, P_ho, true,
1325 true, false));
1326
1327 Vector RM_H(pfes_ho_scalar->GetTrueVSize());
1328 GetTDofsTranspose(*pfes_ho_scalar, M_H, RM_H);
1329 precon.reset(new OperatorJacobiSmoother(RM_H, ess_tdof_list));
1330 }
1331 else if (P_ho)
1332 {
1333 ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
1334 pfes_lor_scalar.get(),
1335 ML_inv_ea));
1336 M_LH.reset(new ProductOperator(M_LH_local_op, P_ho, true, false));
1337
1338 Vector RM_H(pfes_ho_scalar->GetTrueVSize());
1339 GetTDofsTranspose(*pfes_ho_scalar.get(), M_H, RM_H);
1340 precon.reset(new OperatorJacobiSmoother(RM_H, ess_tdof_list));
1341 }
1342 else if (P_lor)
1343 {
1344 Operator *Pt_lor = new TransposeOperator(P_lor);
1345 RML_inv.SetSize(pfes_lor_scalar->GetTrueVSize());
1346 GetTDofsTranspose(*pfes_lor_scalar, ML_inv_ea, RML_inv);
1347 ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
1348 pfes_lor_scalar.get(),
1349 RML_inv));
1350 M_LH.reset(new ProductOperator(Pt_lor, M_LH_local_op, true, true));
1351 R.reset(new ProductOperator(ML_inv_vea.get(), M_LH.get(), false,
1352 false));
1353
1354 precon.reset(new OperatorJacobiSmoother(M_H, ess_tdof_list));
1355 }
1356 else
1357 {
1358 ML_inv_vea.reset(new H1SpaceLumpedMassOperator(pfes_ho_scalar.get(),
1359 pfes_lor_scalar.get(),
1360 ML_inv_ea));
1361 M_LH.reset(M_LH_local_op);
1362
1363 precon.reset(new OperatorJacobiSmoother(M_H, ess_tdof_list));
1364 }
1365 }
1366 R.reset(new ProductOperator(ML_inv_vea.get(), M_LH.get(), false,
1367 false));
1368
1369 TransposeOperator* RT = new TransposeOperator(R.get());
1370 RTxM_LH.reset(new ProductOperator(RT, M_LH.get(), true, false));
1371
1372 SetupPCG();
1373}
1374
1375#endif
1376
1378 const Vector& x, Vector& y) const
1379{
1380 Vector X(fes_ho.GetTrueVSize());
1381 Vector X_dim(R->Width());
1382
1383 Vector Y_dim(R->Height());
1384 Vector Y(fes_lor.GetTrueVSize());
1385
1386 Array<int> vdofs_list;
1387
1388 GetTDofs(fes_ho, x, X);
1389
1390 for (int d = 0; d < fes_ho.GetVDim(); ++d)
1391 {
1392 TDofsListByVDim(fes_ho, d, vdofs_list);
1393 X.GetSubVector(vdofs_list, X_dim);
1394 R->Mult(X_dim, Y_dim);
1395 TDofsListByVDim(fes_lor, d, vdofs_list);
1396 Y.SetSubVector(vdofs_list, Y_dim);
1397 }
1398
1399 SetFromTDofs(fes_lor, Y, y);
1400
1401}
1402
1404 const Vector& x, Vector& y) const
1405{
1406 Vector X(fes_lor.GetTrueVSize());
1407 Vector X_dim(R->Height());
1408
1409 Vector Y_dim(R->Width());
1410 Vector Y(fes_ho.GetTrueVSize());
1411
1412 Array<int> vdofs_list;
1413
1414 GetTDofsTranspose(fes_lor, x, X);
1415
1416 for (int d = 0; d < fes_ho.GetVDim(); ++d)
1417 {
1418 TDofsListByVDim(fes_lor, d, vdofs_list);
1419 X.GetSubVector(vdofs_list, X_dim);
1420 R->MultTranspose(X_dim, Y_dim);
1421 TDofsListByVDim(fes_ho, d, vdofs_list);
1422 Y.SetSubVector(vdofs_list, Y_dim);
1423 }
1424
1425 SetFromTDofsTranspose(fes_ho, Y, y);
1426
1427}
1428
1430 const Vector& x, Vector& y) const
1431{
1432
1433 Vector X(fes_lor.GetTrueVSize());
1434 Vector X_dim(M_LH->Height());
1435 Vector Xbar(pcg.Width());
1436
1437 Vector Y_dim(pcg.Height());
1438 Y_dim = 0.0;
1439 Vector Y(fes_ho.GetTrueVSize());
1440
1441 Array<int> vdofs_list;
1442
1443 GetTDofs(fes_lor, x, X);
1444
1445 for (int d = 0; d < fes_ho.GetVDim(); ++d)
1446 {
1447 TDofsListByVDim(fes_lor, d, vdofs_list);
1448 X.GetSubVector(vdofs_list, X_dim);
1449 // Compute y = P x = (R^T M_LH)^(-1) M_LH^T X = (R^T M_LH)^(-1) Xbar
1450 M_LH->MultTranspose(X_dim, Xbar);
1451 pcg.Mult(Xbar, Y_dim);
1452 TDofsListByVDim(fes_ho, d, vdofs_list);
1453 Y.SetSubVector(vdofs_list, Y_dim);
1454 }
1455
1456 SetFromTDofs(fes_ho, Y, y);
1457
1458}
1459
1461 const Vector& x, Vector& y) const
1462{
1463 Vector X(fes_ho.GetTrueVSize());
1464 Vector X_dim(pcg.Width());
1465 Vector Xbar(pcg.Height());
1466
1467 Vector Y_dim(M_LH->Height());
1468 Vector Y(fes_lor.GetTrueVSize());
1469
1470 Array<int> vdofs_list;
1471
1472 GetTDofsTranspose(fes_ho, x, X);
1473
1474 for (int d = 0; d < fes_ho.GetVDim(); ++d)
1475 {
1476 TDofsListByVDim(fes_ho, d, vdofs_list);
1477 X.GetSubVector(vdofs_list, X_dim);
1478 // Compute y = P^T x = M_LH (R^T M_LH)^(-1) X = M_LH Xbar
1479 Xbar = 0.0;
1480 pcg.Mult(X_dim, Xbar);
1481 M_LH->Mult(Xbar, Y_dim);
1482 TDofsListByVDim(fes_lor, d, vdofs_list);
1483 Y.SetSubVector(vdofs_list, Y_dim);
1484 }
1485
1486 SetFromTDofsTranspose(fes_lor, Y, y);
1487
1488}
1489
1491{
1492 pcg.SetRelTol(p_rtol_);
1493}
1494
1496{
1497 pcg.SetAbsTol(p_atol_);
1498}
1499
1500std::pair<
1501std::unique_ptr<SparseMatrix>,
1502std::unique_ptr<SparseMatrix>>
1504{
1505 std::pair<std::unique_ptr<SparseMatrix>,
1506 std::unique_ptr<SparseMatrix>> r_and_mlh;
1507
1508 Mesh* mesh_ho = fes_ho.GetMesh();
1509 Mesh* mesh_lor = fes_lor.GetMesh();
1510 int nel_ho = mesh_ho->GetNE();
1511 int nel_lor = mesh_lor->GetNE();
1512 int ndof_lor = fes_lor.GetNDofs();
1513
1514 // If the local mesh is empty, skip all computations
1515 if (nel_ho == 0)
1516 {
1517 return std::make_pair(
1518 std::unique_ptr<SparseMatrix>(new SparseMatrix),
1519 std::unique_ptr<SparseMatrix>(new SparseMatrix)
1520 );
1521 }
1522
1523 const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
1524
1525 int nref_max = 0;
1527 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
1528 for (int ig = 0; ig < geoms.Size(); ++ig)
1529 {
1530 Geometry::Type geom = geoms[ig];
1531 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
1532 }
1533
1534 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
1535
1536 // ML_inv contains the inverse lumped (row sum) mass matrix. Note that the
1537 // method will also work with a full (consistent) mass matrix, though this is
1538 // not implemented here. L refers to the low-order refined mesh
1539 Vector ML_inv(ndof_lor);
1540 ML_inv = 0.0;
1541
1542 // Compute ML_inv
1543 for (int iho = 0; iho < nel_ho; ++iho)
1544 {
1545 Array<int> lor_els;
1546 ho2lor.GetRow(iho, lor_els);
1547 int nref = ho2lor.RowSize(iho);
1548
1549 Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
1550 const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
1551 int nedof_lor = fe_lor.GetDof();
1552
1553 // Instead of using a MassIntegrator, manually loop over integration
1554 // points so we can row sum and store the diagonal as a Vector.
1555 Vector ML_el(nedof_lor);
1556 Vector shape_lor(nedof_lor);
1557 Array<int> dofs_lor(nedof_lor);
1558
1559 for (int iref = 0; iref < nref; ++iref)
1560 {
1561 int ilor = lor_els[iref];
1562 ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
1563
1564 int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
1565 const IntegrationRule* ir = &IntRules.Get(geom, order);
1566 ML_el = 0.0;
1567 for (int i = 0; i < ir->GetNPoints(); ++i)
1568 {
1569 const IntegrationPoint& ip_lor = ir->IntPoint(i);
1570 fe_lor.CalcShape(ip_lor, shape_lor);
1571 el_tr->SetIntPoint(&ip_lor);
1572 ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
1573 }
1574 fes_lor.GetElementDofs(ilor, dofs_lor);
1575 ML_inv.AddElementVector(dofs_lor, ML_el);
1576 }
1577 }
1578 // DOF by DOF inverse of non-zero entries
1579 LumpedMassInverse(ML_inv);
1580
1581 // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate
1582 r_and_mlh.first = AllocR();
1583 // Allocate M_LH (same sparsity pattern as R)
1584 // L refers to the low-order refined mesh (DOFs correspond to rows)
1585 // H refers to the higher-order mesh (DOFs correspond to columns)
1586 Memory<int> I(r_and_mlh.first->Height() + 1);
1587 for (int icol = 0; icol < r_and_mlh.first->Height() + 1; ++icol)
1588 {
1589 I[icol] = r_and_mlh.first->GetI()[icol];
1590 }
1591 Memory<int> J(r_and_mlh.first->NumNonZeroElems());
1592 for (int jcol = 0; jcol < r_and_mlh.first->NumNonZeroElems(); ++jcol)
1593 {
1594 J[jcol] = r_and_mlh.first->GetJ()[jcol];
1595 }
1596 r_and_mlh.second = std::unique_ptr<SparseMatrix>(
1597 new SparseMatrix(I, J, NULL, r_and_mlh.first->Height(),
1598 r_and_mlh.first->Width(), true, true, true));
1599
1601 IsoparametricTransformation& emb_tr = ip_tr.Transf;
1602
1603 // Compute M_LH and R
1604 offsets.SetSize(nel_ho+1);
1605 offsets[0] = 0;
1606 for (int iho = 0; iho < nel_ho; ++iho)
1607 {
1608 Array<int> lor_els;
1609 ho2lor.GetRow(iho, lor_els);
1610 int nref = ho2lor.RowSize(iho);
1611
1612 Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
1613 const FiniteElement& fe_ho = *fes_ho.GetFE(iho);
1614 const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
1615 offsets[iho+1] = offsets[iho] + fe_ho.GetDof()*fe_lor.GetDof()*nref;
1616
1617 ElementTransformation *tr_ho = fes_ho.GetElementTransformation(iho);
1618
1619 emb_tr.SetIdentityTransformation(geom);
1620 const DenseTensor& pmats = cf_tr.point_matrices[geom];
1621
1622 int nedof_ho = fe_ho.GetDof();
1623 int nedof_lor = fe_lor.GetDof();
1624 DenseMatrix M_LH_el(nedof_lor, nedof_ho);
1625 DenseMatrix R_el(nedof_lor, nedof_ho);
1626
1627 for (int iref = 0; iref < nref; ++iref)
1628 {
1629 int ilor = lor_els[iref];
1630 ElementTransformation* tr_lor = fes_lor.GetElementTransformation(ilor);
1631
1632 // Create the transformation that embeds the fine low-order element
1633 // within the coarse high-order element in reference space
1634 emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
1635
1636 ElemMixedMass(geom, fe_ho, fe_lor, tr_ho, tr_lor, ip_tr, M_LH_el);
1637
1638 Array<int> dofs_lor(nedof_lor);
1639 fes_lor.GetElementDofs(ilor, dofs_lor);
1640 Vector R_row;
1641 for (int i = 0; i < nedof_lor; ++i)
1642 {
1643 M_LH_el.GetRow(i, R_row);
1644 R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
1645 }
1646 Array<int> dofs_ho(nedof_ho);
1647 fes_ho.GetElementDofs(iho, dofs_ho);
1648 r_and_mlh.second->AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
1649 r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
1650
1651 }
1652 }
1653
1654 return r_and_mlh;
1655}
1656
1658 const FiniteElementSpace& fes, const Vector& x, Vector& X) const
1659{
1660 const Operator* res = fes.GetRestrictionOperator();
1661 if (res)
1662 {
1663 res->Mult(x, X);
1664 }
1665 else
1666 {
1667 X = x;
1668 }
1669}
1670
1672 const FiniteElementSpace& fes, const Vector &X, Vector& x) const
1673{
1674 const Operator* P = fes.GetProlongationMatrix();
1675 if (P)
1676 {
1677 P->Mult(X, x);
1678 }
1679 else
1680 {
1681 x = X;
1682 }
1683}
1684
1686 const FiniteElementSpace& fes, const Vector& x, Vector& X) const
1687{
1688 const Operator* P = fes.GetProlongationMatrix();
1689 if (P)
1690 {
1691 P->MultTranspose(x, X);
1692 }
1693 else
1694 {
1695 X = x;
1696 }
1697}
1698
1700 const FiniteElementSpace& fes, const Vector &X, Vector& x) const
1701{
1702 const Operator *R_op = fes.GetRestrictionOperator();
1703 if (R_op)
1704 {
1705 R_op->MultTranspose(X, x);
1706 }
1707 else
1708 {
1709 x = X;
1710 }
1711}
1712
1714 const FiniteElementSpace& fes, int vdim, Array<int>& vdofs_list) const
1715{
1716 const SparseMatrix *R_mat = fes.GetRestrictionMatrix();
1717 if (R_mat)
1718 {
1719 Array<int> x_vdofs_list(fes.GetNDofs());
1720 Array<int> x_vdofs_marker(fes.GetVSize());
1721 Array<int> X_vdofs_marker(fes.GetTrueVSize());
1722 fes.GetVDofs(vdim, x_vdofs_list);
1723 FiniteElementSpace::ListToMarker(x_vdofs_list, fes.GetVSize(), x_vdofs_marker);
1724 R_mat->BooleanMult(x_vdofs_marker, X_vdofs_marker);
1725 FiniteElementSpace::MarkerToList(X_vdofs_marker, vdofs_list);
1726 }
1727 else
1728 {
1729 vdofs_list.SetSize(fes.GetNDofs());
1730 fes.GetVDofs(vdim, vdofs_list);
1731 }
1732}
1733
1735 Vector& ML_inv) const
1736{
1737#ifdef MFEM_USE_MPI
1738 // LumpedMassInverse may get called from serial and MPI parallel routines
1739 // since we do not know which code path is calling it we must check if
1740 // the pfes pointer is null when MPI is available.
1741 auto * fes = pfes_lor_scalar == nullptr ? fes_lor_scalar.get() :
1742 pfes_lor_scalar.get();
1743#else
1744 auto * fes = fes_lor_scalar.get();
1745#endif
1746 MFEM_ASSERT(fes != nullptr, "[p]fes_lor_scalar is nullptr");
1747
1748 Vector ML_inv_true(fes->GetTrueVSize());
1749 const Operator *P = fes->GetProlongationMatrix();
1750 if (P) { P->MultTranspose(ML_inv, ML_inv_true); }
1751 else { ML_inv_true = ML_inv; }
1752
1753 ML_inv_true.Reciprocal();
1754
1755 if (P) { P->Mult(ML_inv_true, ML_inv); }
1756 else { ML_inv = ML_inv_true; }
1757
1758}
1759
1760std::unique_ptr<SparseMatrix>
1762{
1763 const Table& elem_dof_ho = fes_ho.GetElementToDofTable();
1764 const Table& elem_dof_lor = fes_lor.GetElementToDofTable();
1765 const int ndof_ho = fes_ho.GetNDofs();
1766 const int ndof_lor = fes_lor.GetNDofs();
1767
1768 Table dof_elem_lor;
1769 Transpose(elem_dof_lor, dof_elem_lor, ndof_lor);
1770
1771 Mesh* mesh_lor = fes_lor.GetMesh();
1772 const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
1773
1774 // mfem::Mult but uses ho2lor to map HO elements to LOR elements
1775 const int* elem_dof_hoI = elem_dof_ho.GetI();
1776 const int* elem_dof_hoJ = elem_dof_ho.GetJ();
1777 const int* dof_elem_lorI = dof_elem_lor.GetI();
1778 const int* dof_elem_lorJ = dof_elem_lor.GetJ();
1779
1780 Array<int> I(ndof_lor + 1);
1781
1782 // figure out the size of J
1783 Array<int> dof_used_ho;
1784 dof_used_ho.SetSize(ndof_ho, -1);
1785
1786 int sizeJ = 0;
1787 for (int ilor = 0; ilor < ndof_lor; ++ilor)
1788 {
1789 for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
1790 {
1791 int el_lor = dof_elem_lorJ[jlor];
1792 int iho = cf_tr.embeddings[el_lor].parent;
1793 for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
1794 {
1795 int dof_ho = elem_dof_hoJ[jho];
1796 if (dof_used_ho[dof_ho] != ilor)
1797 {
1798 dof_used_ho[dof_ho] = ilor;
1799 ++sizeJ;
1800 }
1801 }
1802 }
1803 }
1804
1805 // initialize dof_ho_dof_lor
1806 Table dof_lor_dof_ho;
1807 dof_lor_dof_ho.SetDims(ndof_lor, sizeJ);
1808
1809 for (int i = 0; i < ndof_ho; ++i)
1810 {
1811 dof_used_ho[i] = -1;
1812 }
1813
1814 // set values of J
1815 int* dof_dofI = dof_lor_dof_ho.GetI();
1816 int* dof_dofJ = dof_lor_dof_ho.GetJ();
1817 sizeJ = 0;
1818 for (int ilor = 0; ilor < ndof_lor; ++ilor)
1819 {
1820 dof_dofI[ilor] = sizeJ;
1821 for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
1822 {
1823 int el_lor = dof_elem_lorJ[jlor];
1824 int iho = cf_tr.embeddings[el_lor].parent;
1825 for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
1826 {
1827 int dof_ho = elem_dof_hoJ[jho];
1828 if (dof_used_ho[dof_ho] != ilor)
1829 {
1830 dof_used_ho[dof_ho] = ilor;
1831 dof_dofJ[sizeJ] = dof_ho;
1832 ++sizeJ;
1833 }
1834 }
1835 }
1836 }
1837
1838 dof_lor_dof_ho.SortRows();
1839 real_t* data = Memory<real_t>(dof_dofI[ndof_lor]);
1840
1841 std::unique_ptr<SparseMatrix> R_local(new SparseMatrix(
1842 dof_dofI, dof_dofJ, data, ndof_lor,
1843 ndof_ho, true, true, true));
1844 (*R_local) = 0.0;
1845
1846 dof_lor_dof_ho.LoseData();
1847
1848 return R_local;
1849}
1850
1852 const FiniteElementSpace* fes_ho_, const FiniteElementSpace* fes_lor_,
1853 Table* ho2lor_, Vector* M_LH_ea_) :
1854 Operator(fes_lor_->GetElementRestriction(ElementDofOrdering::NATIVE)->Width(),
1855 fes_ho_->GetElementRestriction(ElementDofOrdering::NATIVE)->Width()),
1856 fes_ho(fes_ho_), fes_lor(fes_lor_), ho2lor(ho2lor_),
1857 M_LH_ea(M_LH_ea_)
1858{ }
1859
1861 Vector &y) const
1862{
1863 const Operator* elem_restrict_ho = fes_ho->GetElementRestriction(
1865 const Operator* elem_restrict_lor = fes_lor->GetElementRestriction(
1867
1868 const int vdim = fes_ho->GetVDim();
1869 const int iho = 0;
1870 const int nref = ho2lor->RowSize(iho);
1871 const int ndof_ho = fes_ho->GetFE(iho)->GetDof();
1872 const int ndof_lor = fes_lor->GetFE(ho2lor->GetRow(iho)[0])->GetDof();
1873 const Mesh *mesh_ho = fes_ho->GetMesh();
1874 const int nel_ho = mesh_ho->GetNE();
1875
1876 Vector tempx(elem_restrict_ho->Height());
1877 elem_restrict_ho->Mult(x, tempx);
1878
1879 Vector tempy(ndof_lor*nref*vdim*nel_ho);
1880
1881 auto v_M_mixed_ea = Reshape(M_LH_ea->Read(), ndof_lor, ndof_ho, nref,
1882 nel_ho);
1883 auto v_tempx = Reshape(tempx.Read(), ndof_ho, vdim, nel_ho);
1884 auto v_tempy = Reshape(tempy.Write(), ndof_lor, nref, vdim, nel_ho);
1885
1886
1887 mfem::forall(ndof_lor * nref * vdim * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
1888 {
1889 const int j = tid % ndof_lor;
1890 const int i = (tid / ndof_lor) % nref;
1891 const int v = (tid / (ndof_lor * nref)) % vdim;
1892 const int iho = (tid / (ndof_lor * nref * vdim)) % nel_ho;
1893
1894 real_t dot = 0.0;
1895 for (int k=0; k<ndof_ho; ++k)
1896 {
1897 dot += v_M_mixed_ea(j, k, i, iho) * v_tempx(k, v, iho);
1898 }
1899
1900 v_tempy(j, i, v, iho) = dot;
1901 });
1902
1903 elem_restrict_lor->MultTranspose(tempy, y);
1904}
1905
1907 const Vector &x, Vector &y) const
1908{
1909 const Operator* elem_restrict_ho = fes_ho->GetElementRestriction(
1911 const Operator* elem_restrict_lor = fes_lor->GetElementRestriction(
1913
1914 const int vdim = fes_ho->GetVDim();
1915 const int iho = 0;
1916 const int nref = ho2lor->RowSize(iho);
1917 const int ndof_ho = fes_ho->GetFE(iho)->GetDof();
1918 const int ndof_lor = fes_lor->GetFE(ho2lor->GetRow(iho)[0])->GetDof();
1919 const Mesh *mesh_ho = fes_ho->GetMesh();
1920 const int nel_ho = mesh_ho->GetNE();
1921
1922 Vector tempx(elem_restrict_lor->Height());
1923 elem_restrict_lor->Mult(x, tempx);
1924
1925 Vector tempy(ndof_ho*vdim*nel_ho);
1926
1927 auto v_M_mixed_ea = Reshape(M_LH_ea->Read(), ndof_lor, ndof_ho, nref,
1928 nel_ho);
1929 auto v_tempx = Reshape(tempx.Read(), ndof_lor, nref, vdim, nel_ho);
1930 auto v_tempy = Reshape(tempy.Write(), ndof_ho, vdim, nel_ho);
1931
1932 mfem::forall(ndof_ho * vdim * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
1933 {
1934 const int k = tid % ndof_ho;
1935 const int v = (tid / ndof_ho) % vdim;
1936 const int iho = (tid / (ndof_ho * vdim)) % nel_ho;
1937
1938 real_t dot = 0.0;
1939 for (int i=0; i<nref; ++i)
1940 {
1941 for (int j=0; j<ndof_lor; ++j)
1942 {
1943 dot += v_M_mixed_ea(j, k, i, iho) * v_tempx(j, i, v, iho);
1944 }
1945 v_tempy(k, v, iho) = dot;
1946 }
1947 });
1948
1949 elem_restrict_ho->MultTranspose(tempy, y);
1950}
1951
1953 const FiniteElementSpace* fes_ho_,
1954 const FiniteElementSpace* fes_lor_,
1955 Vector& ML_inv_) :
1956 Operator(ML_inv_.Size(), ML_inv_.Size()),
1957 fes_ho(fes_ho_), fes_lor(fes_lor_),
1958 ML_inv(&ML_inv_)
1959{ }
1960
1962 Vector &y) const
1963{
1964 MFEM_ASSERT(ML_inv->Size() == x.Size(), "sizes not the same");
1965 auto v_ML_inv = Reshape(ML_inv->Read(), ML_inv->Size());
1966 auto v_x = Reshape(x.Read(), x.Size());
1967 auto v_y = Reshape(y.Write(), y.Size());
1968
1969 mfem::forall(ML_inv->Size(), [=] MFEM_HOST_DEVICE(int i)
1970 { v_y(i) = v_ML_inv(i) * v_x(i); });
1971}
1972
1974 const Vector &x, Vector &y) const
1975{
1976 this->Mult(x,y); // lumped diagonal has the same Mult and MultTranspose behavior
1977}
1978
1980{
1981 delete F;
1982 delete B;
1983}
1984
1986{
1987 if (!F) { BuildF(); }
1988 return *F;
1989}
1990
1992{
1993 if (!B)
1994 {
1995 if (!F) { BuildF(); }
1996 B = new L2Prolongation(*F);
1997 }
1998 return *B;
1999}
2000
2001void L2ProjectionGridTransfer::BuildF()
2002{
2003 if (!force_l2_space &&
2005 {
2006 if (!Parallel())
2007 {
2008 F = new L2ProjectionH1Space(dom_fes, ran_fes,
2009 use_ea, d_mt);
2010 }
2011 else
2012 {
2013#ifdef MFEM_USE_MPI
2014 const mfem::ParFiniteElementSpace& dom_pfes =
2015 static_cast<mfem::ParFiniteElementSpace&>(dom_fes);
2016 const mfem::ParFiniteElementSpace& ran_pfes =
2017 static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
2018 F = new L2ProjectionH1Space(dom_pfes, ran_pfes,
2019 use_ea, d_mt);
2020#endif
2021 }
2022 }
2023 else
2024 {
2025 F = new L2ProjectionL2Space(dom_fes, ran_fes,
2026 use_ea, d_mt);
2027 }
2028}
2029
2034
2035
2037 const FiniteElementSpace& hFESpace_)
2038 : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize())
2039{
2040 bool isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
2041 if (lFESpace_.FEColl() == hFESpace_.FEColl() && !isvar_order)
2042 {
2044 hFESpace_.GetTransferOperator(lFESpace_, P);
2045 P.SetOperatorOwner(false);
2046 opr = P.Ptr();
2047 }
2048 else if (lFESpace_.GetVDim() == 1
2049 && hFESpace_.GetVDim() == 1
2050 && dynamic_cast<const TensorBasisElement*>(lFESpace_.GetTypicalFE())
2051 && dynamic_cast<const TensorBasisElement*>(hFESpace_.GetTypicalFE())
2052 && !isvar_order
2053 && (hFESpace_.FEColl()->GetContType() ==
2055 hFESpace_.FEColl()->GetContType() ==
2057 {
2058 opr = new TensorProductPRefinementTransferOperator(lFESpace_, hFESpace_);
2059 }
2060 else
2061 {
2062 opr = new PRefinementTransferOperator(lFESpace_, hFESpace_);
2063 }
2064}
2065
2067
2068void TransferOperator::Mult(const Vector& x, Vector& y) const
2069{
2070 opr->Mult(x, y);
2071}
2072
2074{
2075 opr->MultTranspose(x, y);
2076}
2077
2078
2080 const FiniteElementSpace& lFESpace_, const FiniteElementSpace& hFESpace_)
2081 : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
2082 hFESpace(hFESpace_)
2083{
2084 isvar_order = lFESpace_.IsVariableOrder() || hFESpace_.IsVariableOrder();
2085}
2086
2088{
2089 Mesh* mesh = hFESpace.GetMesh();
2090 Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
2091 DenseMatrix loc_prol;
2092 Vector subY, subX;
2093
2094 Geometry::Type cached_geom = Geometry::INVALID;
2095 const FiniteElement* h_fe = NULL;
2096 const FiniteElement* l_fe = NULL;
2098
2099 int vdim = lFESpace.GetVDim();
2100
2101 y = 0.0;
2102
2103 for (int i = 0; i < mesh->GetNE(); i++)
2104 {
2105 DofTransformation * doftrans_h = hFESpace.GetElementDofs(i, h_dofs);
2106 DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
2107
2108 const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
2109 if (geom != cached_geom || isvar_order)
2110 {
2111 h_fe = hFESpace.GetFE(i);
2112 l_fe = lFESpace.GetFE(i);
2114 h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
2115 subY.SetSize(loc_prol.Height());
2116 cached_geom = geom;
2117 }
2118
2119 for (int vd = 0; vd < vdim; vd++)
2120 {
2121 l_dofs.Copy(l_vdofs);
2122 lFESpace.DofsToVDofs(vd, l_vdofs);
2123 h_dofs.Copy(h_vdofs);
2124 hFESpace.DofsToVDofs(vd, h_vdofs);
2125 x.GetSubVector(l_vdofs, subX);
2126 if (doftrans_l)
2127 {
2128 doftrans_l->InvTransformPrimal(subX);
2129 }
2130 loc_prol.Mult(subX, subY);
2131 if (doftrans_h)
2132 {
2133 doftrans_h->TransformPrimal(subY);
2134 }
2135 y.SetSubVector(h_vdofs, subY);
2136 }
2137 }
2138}
2139
2141 Vector& y) const
2142{
2143 y = 0.0;
2144
2145 Mesh* mesh = hFESpace.GetMesh();
2146 Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
2147 DenseMatrix loc_prol;
2148 Vector subY, subX;
2149
2150 Array<char> processed(hFESpace.GetVSize());
2151 processed = 0;
2152
2153 Geometry::Type cached_geom = Geometry::INVALID;
2154 const FiniteElement* h_fe = NULL;
2155 const FiniteElement* l_fe = NULL;
2157
2158 int vdim = lFESpace.GetVDim();
2159
2160 for (int i = 0; i < mesh->GetNE(); i++)
2161 {
2162 DofTransformation * doftrans_h = hFESpace.GetElementDofs(i, h_dofs);
2163 DofTransformation * doftrans_l = lFESpace.GetElementDofs(i, l_dofs);
2164
2165 const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
2166 if (geom != cached_geom || isvar_order)
2167 {
2168 h_fe = hFESpace.GetFE(i);
2169 l_fe = lFESpace.GetFE(i);
2171 h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
2172 loc_prol.Transpose();
2173 subY.SetSize(loc_prol.Height());
2174 cached_geom = geom;
2175 }
2176
2177 for (int vd = 0; vd < vdim; vd++)
2178 {
2179 l_dofs.Copy(l_vdofs);
2180 lFESpace.DofsToVDofs(vd, l_vdofs);
2181 h_dofs.Copy(h_vdofs);
2182 hFESpace.DofsToVDofs(vd, h_vdofs);
2183
2184 x.GetSubVector(h_vdofs, subX);
2185 if (doftrans_h)
2186 {
2187 doftrans_h->InvTransformDual(subX);
2188 }
2189 for (int p = 0; p < h_dofs.Size(); ++p)
2190 {
2191 if (processed[lFESpace.DecodeDof(h_dofs[p])])
2192 {
2193 subX[p] = 0.0;
2194 }
2195 }
2196
2197 loc_prol.Mult(subX, subY);
2198 if (doftrans_l)
2199 {
2200 doftrans_l->TransformDual(subY);
2201 }
2202 y.AddElementVector(l_vdofs, subY);
2203 }
2204
2205 for (int p = 0; p < h_dofs.Size(); ++p)
2206 {
2207 processed[lFESpace.DecodeDof(h_dofs[p])] = 1;
2208 }
2209 }
2210}
2211
2212
2215 const FiniteElementSpace& lFESpace_,
2216 const FiniteElementSpace& hFESpace_)
2217 : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
2218 hFESpace(hFESpace_)
2219{
2220 // Assuming the same element type
2221 Mesh* mesh = lFESpace.GetMesh();
2222 dim = mesh->Dimension();
2223 const FiniteElement& el = *lFESpace.GetTypicalFE();
2224
2225 const TensorBasisElement* ltel =
2226 dynamic_cast<const TensorBasisElement*>(&el);
2227 MFEM_VERIFY(ltel, "Low order FE space must be tensor product space");
2228
2229 const TensorBasisElement* htel =
2230 dynamic_cast<const TensorBasisElement*>(hFESpace.GetTypicalFE());
2231 MFEM_VERIFY(htel, "High order FE space must be tensor product space");
2232 const Array<int>& hdofmap = htel->GetDofMap();
2233
2234 const IntegrationRule& ir = hFESpace.GetTypicalFE()->GetNodes();
2235 IntegrationRule irLex = ir;
2236
2237 // The quadrature points, or equivalently, the dofs of the high order space
2238 // must be sorted in lexicographical order
2239 for (int i = 0; i < ir.GetNPoints(); ++i)
2240 {
2241 int j = hdofmap[i] >=0 ? hdofmap[i] : -1 - hdofmap[i];
2242 irLex.IntPoint(i) = ir.IntPoint(j);
2243 }
2244
2245 NE = lFESpace.GetNE();
2246 const DofToQuad& maps = el.GetDofToQuad(irLex, DofToQuad::TENSOR);
2247
2248 D1D = maps.ndof;
2249 Q1D = maps.nqpt;
2250 B = maps.B;
2251 Bt = maps.Bt;
2252
2253 elem_restrict_lex_l =
2255
2256 MFEM_VERIFY(elem_restrict_lex_l,
2257 "Low order ElementRestriction not available");
2258
2259 elem_restrict_lex_h =
2261
2262 MFEM_VERIFY(elem_restrict_lex_h,
2263 "High order ElementRestriction not available");
2264
2265 localL.SetSize(elem_restrict_lex_l->Height(), Device::GetMemoryType());
2266 localH.SetSize(elem_restrict_lex_h->Height(), Device::GetMemoryType());
2267 localL.UseDevice(true);
2268 localH.UseDevice(true);
2269
2270 MFEM_VERIFY(dynamic_cast<const ElementRestriction*>(elem_restrict_lex_h),
2271 "High order element restriction is of unsupported type");
2272
2273 mask.SetSize(localH.Size(), Device::GetMemoryType());
2274 static_cast<const ElementRestriction*>(elem_restrict_lex_h)
2275 ->BooleanMask(mask);
2276 mask.UseDevice(true);
2277}
2278
2279namespace TransferKernels
2280{
2281void Prolongation2D(const int NE, const int D1D, const int Q1D,
2282 const Vector& localL, Vector& localH,
2283 const Array<real_t>& B, const Vector& mask)
2284{
2285 auto x_ = Reshape(localL.Read(), D1D, D1D, NE);
2286 auto y_ = Reshape(localH.Write(), Q1D, Q1D, NE);
2287 auto B_ = Reshape(B.Read(), Q1D, D1D);
2288 auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
2289
2290 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
2291 {
2292 for (int qy = 0; qy < Q1D; ++qy)
2293 {
2294 for (int qx = 0; qx < Q1D; ++qx)
2295 {
2296 y_(qx, qy, e) = 0.0;
2297 }
2298 }
2299
2300 for (int dy = 0; dy < D1D; ++dy)
2301 {
2302 real_t sol_x[DofQuadLimits::MAX_Q1D];
2303 for (int qy = 0; qy < Q1D; ++qy)
2304 {
2305 sol_x[qy] = 0.0;
2306 }
2307 for (int dx = 0; dx < D1D; ++dx)
2308 {
2309 const real_t s = x_(dx, dy, e);
2310 for (int qx = 0; qx < Q1D; ++qx)
2311 {
2312 sol_x[qx] += B_(qx, dx) * s;
2313 }
2314 }
2315 for (int qy = 0; qy < Q1D; ++qy)
2316 {
2317 const real_t d2q = B_(qy, dy);
2318 for (int qx = 0; qx < Q1D; ++qx)
2319 {
2320 y_(qx, qy, e) += d2q * sol_x[qx];
2321 }
2322 }
2323 }
2324 for (int qy = 0; qy < Q1D; ++qy)
2325 {
2326 for (int qx = 0; qx < Q1D; ++qx)
2327 {
2328 y_(qx, qy, e) *= m_(qx, qy, e);
2329 }
2330 }
2331 });
2332}
2333
2334void Prolongation3D(const int NE, const int D1D, const int Q1D,
2335 const Vector& localL, Vector& localH,
2336 const Array<real_t>& B, const Vector& mask)
2337{
2338 auto x_ = Reshape(localL.Read(), D1D, D1D, D1D, NE);
2339 auto y_ = Reshape(localH.Write(), Q1D, Q1D, Q1D, NE);
2340 auto B_ = Reshape(B.Read(), Q1D, D1D);
2341 auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
2342
2343 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
2344 {
2345 for (int qz = 0; qz < Q1D; ++qz)
2346 {
2347 for (int qy = 0; qy < Q1D; ++qy)
2348 {
2349 for (int qx = 0; qx < Q1D; ++qx)
2350 {
2351 y_(qx, qy, qz, e) = 0.0;
2352 }
2353 }
2354 }
2355
2356 for (int dz = 0; dz < D1D; ++dz)
2357 {
2358 real_t sol_xy[DofQuadLimits::MAX_Q1D][DofQuadLimits::MAX_Q1D];
2359 for (int qy = 0; qy < Q1D; ++qy)
2360 {
2361 for (int qx = 0; qx < Q1D; ++qx)
2362 {
2363 sol_xy[qy][qx] = 0.0;
2364 }
2365 }
2366 for (int dy = 0; dy < D1D; ++dy)
2367 {
2368 real_t sol_x[DofQuadLimits::MAX_Q1D];
2369 for (int qx = 0; qx < Q1D; ++qx)
2370 {
2371 sol_x[qx] = 0;
2372 }
2373 for (int dx = 0; dx < D1D; ++dx)
2374 {
2375 const real_t s = x_(dx, dy, dz, e);
2376 for (int qx = 0; qx < Q1D; ++qx)
2377 {
2378 sol_x[qx] += B_(qx, dx) * s;
2379 }
2380 }
2381 for (int qy = 0; qy < Q1D; ++qy)
2382 {
2383 const real_t wy = B_(qy, dy);
2384 for (int qx = 0; qx < Q1D; ++qx)
2385 {
2386 sol_xy[qy][qx] += wy * sol_x[qx];
2387 }
2388 }
2389 }
2390 for (int qz = 0; qz < Q1D; ++qz)
2391 {
2392 const real_t wz = B_(qz, dz);
2393 for (int qy = 0; qy < Q1D; ++qy)
2394 {
2395 for (int qx = 0; qx < Q1D; ++qx)
2396 {
2397 y_(qx, qy, qz, e) += wz * sol_xy[qy][qx];
2398 }
2399 }
2400 }
2401 }
2402 for (int qz = 0; qz < Q1D; ++qz)
2403 {
2404 for (int qy = 0; qy < Q1D; ++qy)
2405 {
2406 for (int qx = 0; qx < Q1D; ++qx)
2407 {
2408 y_(qx, qy, qz, e) *= m_(qx, qy, qz, e);
2409 }
2410 }
2411 }
2412 });
2413}
2414
2415void Restriction2D(const int NE, const int D1D, const int Q1D,
2416 const Vector& localH, Vector& localL,
2417 const Array<real_t>& Bt, const Vector& mask)
2418{
2419 auto x_ = Reshape(localH.Read(), Q1D, Q1D, NE);
2420 auto y_ = Reshape(localL.Write(), D1D, D1D, NE);
2421 auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
2422 auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
2423
2424 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
2425 {
2426 for (int dy = 0; dy < D1D; ++dy)
2427 {
2428 for (int dx = 0; dx < D1D; ++dx)
2429 {
2430 y_(dx, dy, e) = 0.0;
2431 }
2432 }
2433
2434 for (int qy = 0; qy < Q1D; ++qy)
2435 {
2436 real_t sol_x[DofQuadLimits::MAX_D1D];
2437 for (int dx = 0; dx < D1D; ++dx)
2438 {
2439 sol_x[dx] = 0.0;
2440 }
2441 for (int qx = 0; qx < Q1D; ++qx)
2442 {
2443 const real_t s = m_(qx, qy, e) * x_(qx, qy, e);
2444 for (int dx = 0; dx < D1D; ++dx)
2445 {
2446 sol_x[dx] += Bt_(dx, qx) * s;
2447 }
2448 }
2449 for (int dy = 0; dy < D1D; ++dy)
2450 {
2451 const real_t q2d = Bt_(dy, qy);
2452 for (int dx = 0; dx < D1D; ++dx)
2453 {
2454 y_(dx, dy, e) += q2d * sol_x[dx];
2455 }
2456 }
2457 }
2458 });
2459}
2460void Restriction3D(const int NE, const int D1D, const int Q1D,
2461 const Vector& localH, Vector& localL,
2462 const Array<real_t>& Bt, const Vector& mask)
2463{
2464 auto x_ = Reshape(localH.Read(), Q1D, Q1D, Q1D, NE);
2465 auto y_ = Reshape(localL.Write(), D1D, D1D, D1D, NE);
2466 auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
2467 auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
2468
2469 mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
2470 {
2471 for (int dz = 0; dz < D1D; ++dz)
2472 {
2473 for (int dy = 0; dy < D1D; ++dy)
2474 {
2475 for (int dx = 0; dx < D1D; ++dx)
2476 {
2477 y_(dx, dy, dz, e) = 0.0;
2478 }
2479 }
2480 }
2481
2482 for (int qz = 0; qz < Q1D; ++qz)
2483 {
2484 real_t sol_xy[DofQuadLimits::MAX_D1D][DofQuadLimits::MAX_D1D];
2485 for (int dy = 0; dy < D1D; ++dy)
2486 {
2487 for (int dx = 0; dx < D1D; ++dx)
2488 {
2489 sol_xy[dy][dx] = 0;
2490 }
2491 }
2492 for (int qy = 0; qy < Q1D; ++qy)
2493 {
2494 real_t sol_x[DofQuadLimits::MAX_D1D];
2495 for (int dx = 0; dx < D1D; ++dx)
2496 {
2497 sol_x[dx] = 0;
2498 }
2499 for (int qx = 0; qx < Q1D; ++qx)
2500 {
2501 const real_t s = m_(qx, qy, qz, e) * x_(qx, qy, qz, e);
2502 for (int dx = 0; dx < D1D; ++dx)
2503 {
2504 sol_x[dx] += Bt_(dx, qx) * s;
2505 }
2506 }
2507 for (int dy = 0; dy < D1D; ++dy)
2508 {
2509 const real_t wy = Bt_(dy, qy);
2510 for (int dx = 0; dx < D1D; ++dx)
2511 {
2512 sol_xy[dy][dx] += wy * sol_x[dx];
2513 }
2514 }
2515 }
2516 for (int dz = 0; dz < D1D; ++dz)
2517 {
2518 const real_t wz = Bt_(dz, qz);
2519 for (int dy = 0; dy < D1D; ++dy)
2520 {
2521 for (int dx = 0; dx < D1D; ++dx)
2522 {
2523 y_(dx, dy, dz, e) += wz * sol_xy[dy][dx];
2524 }
2525 }
2526 }
2527 }
2528 });
2529}
2530} // namespace TransferKernels
2531
2533 Vector& y) const
2534{
2535 if (lFESpace.GetMesh()->GetNE() == 0)
2536 {
2537 return;
2538 }
2539
2540 elem_restrict_lex_l->Mult(x, localL);
2541 if (dim == 2)
2542 {
2543 TransferKernels::Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
2544 }
2545 else if (dim == 3)
2546 {
2547 TransferKernels::Prolongation3D(NE, D1D, Q1D, localL, localH, B, mask);
2548 }
2549 else
2550 {
2551 MFEM_ABORT("TensorProductPRefinementTransferOperator::Mult not "
2552 "implemented for dim = "
2553 << dim);
2554 }
2555 elem_restrict_lex_h->MultTranspose(localH, y);
2556}
2557
2559 Vector& y) const
2560{
2561 if (lFESpace.GetMesh()->GetNE() == 0)
2562 {
2563 return;
2564 }
2565
2566 elem_restrict_lex_h->Mult(x, localH);
2567 if (dim == 2)
2568 {
2569 TransferKernels::Restriction2D(NE, D1D, Q1D, localH, localL, Bt, mask);
2570 }
2571 else if (dim == 3)
2572 {
2573 TransferKernels::Restriction3D(NE, D1D, Q1D, localH, localL, Bt, mask);
2574 }
2575 else
2576 {
2577 MFEM_ABORT("TensorProductPRefinementTransferOperator::MultTranspose not "
2578 "implemented for dim = "
2579 << dim);
2580 }
2581 elem_restrict_lex_l->MultTranspose(localL, y);
2582}
2583
2584
2586 const FiniteElementSpace& hFESpace_)
2587 : Operator(hFESpace_.GetTrueVSize(), lFESpace_.GetTrueVSize()),
2588 lFESpace(lFESpace_),
2589 hFESpace(hFESpace_)
2590{
2591 localTransferOperator = new TransferOperator(lFESpace_, hFESpace_);
2592
2593 P = lFESpace.GetProlongationMatrix();
2594 R = hFESpace.IsVariableOrder() ? hFESpace.GetHpRestrictionMatrix() :
2595 hFESpace.GetRestrictionMatrix();
2596
2597 // P and R can be both null
2598 // P can be null and R not null
2599 // If P is not null it is assumed that R is not null as well
2600 if (P) { MFEM_VERIFY(R, "Both P and R have to be not NULL") }
2601
2602 if (P)
2603 {
2604 tmpL.SetSize(lFESpace_.GetVSize());
2605 tmpH.SetSize(hFESpace_.GetVSize());
2606 }
2607 // P can be null and R not null
2608 else if (R)
2609 {
2610 tmpH.SetSize(hFESpace_.GetVSize());
2611 }
2612}
2613
2615{
2616 delete localTransferOperator;
2617}
2618
2620{
2621 if (P)
2622 {
2623 P->Mult(x, tmpL);
2624 localTransferOperator->Mult(tmpL, tmpH);
2625 R->Mult(tmpH, y);
2626 }
2627 else if (R)
2628 {
2629 localTransferOperator->Mult(x, tmpH);
2630 R->Mult(tmpH, y);
2631 }
2632 else
2633 {
2634 localTransferOperator->Mult(x, y);
2635 }
2636}
2637
2639{
2640 if (P)
2641 {
2642 R->MultTranspose(x, tmpH);
2643 localTransferOperator->MultTranspose(tmpH, tmpL);
2644 P->MultTranspose(tmpL, y);
2645 }
2646 else if (R)
2647 {
2648 R->MultTranspose(x, tmpH);
2649 localTransferOperator->MultTranspose(tmpH, y);
2650 }
2651 else
2652 {
2653 localTransferOperator->MultTranspose(x, y);
2654 }
2655}
2656
2657} // namespace mfem
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:758
int Size() const
Return the logical size of the array.
Definition array.hpp:147
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition array.hpp:935
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:337
static void Mult(const DenseTensor &A, const Vector &x, Vector &y)
Computes (e.g. by calling AddMult(A,x,y,1,0,Op::N)).
Definition batched.cpp:60
static void MultTranspose(const DenseTensor &A, const Vector &x, Vector &y)
Computes (e.g. by calling AddMult(A,x,y,1,0,Op::T)).
Definition batched.cpp:65
static void Invert(DenseTensor &A)
Replaces the block diagonal matrix with its inverse .
Definition batched.cpp:71
Abstract base class BilinearFormIntegrator.
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication: .
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Definition solvers.hpp:551
Data type for scaled Jacobi-type smoother of sparse matrix.
void Factor()
Factor the current DenseMatrix, *a.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication with the inverse of dense matrix.
void GetInverseMatrix(DenseMatrix &Ainv) const
Compute and return the inverse matrix in Ainv.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:125
void Transpose()
(*this) = (*this)^t
void SetRow(int r, const real_t *row)
real_t * GetData() const
Returns the matrix data array.
Definition densemat.hpp:118
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
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.
void GetRow(int r, Vector &row) const
Rank 3 tensor (array of matrices)
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Memory< real_t > & GetMemory()
void NewMemoryAndSize(const Memory< real_t > &mem, int i, int j, int k, bool own_mem)
Reset the DenseTensor to use the given external Memory mem and dimensions i, j, and k.
const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(GetMemory(), TotalSize(), on_dev).
int SizeI() const
int SizeK() const
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Definition device.hpp:278
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
Definition fe_base.hpp:141
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition fe_base.hpp:165
Array< real_t > B
Basis functions evaluated at quadrature points.
Definition fe_base.hpp:193
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:178
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition fe_base.hpp:182
Array< real_t > Bt
Transpose of B.
Definition fe_base.hpp:199
void TransformDual(real_t *v) const
Definition doftrans.cpp:81
void InvTransformPrimal(real_t *v) const
Definition doftrans.cpp:49
void InvTransformDual(real_t *v) const
Definition doftrans.cpp:113
void TransformPrimal(real_t *v) const
Definition doftrans.cpp:17
Operator that converts FiniteElementSpace L-vectors to E-vectors.
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:144
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:106
virtual int GetContType() const =0
@ DISCONTINUOUS
Field is discontinuous across element interfaces.
Definition fe_coll.hpp:48
@ CONTINUOUS
Field is continuous across element interfaces.
Definition fe_coll.hpp:45
Derefinement operator, used by the friend class InterpolationGridTransfer.
Definition fespace.hpp:558
GridFunction interpolation operator applicable after mesh refinement.
Definition fespace.hpp:532
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
void GetVDofs(int vd, Array< int > &dofs, int ndofs=-1) const
Returns the indices of all of the VDofs for the specified dimension 'vd'.
Definition fespace.cpp:238
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.cpp:1463
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:709
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:1287
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
Definition fespace.cpp:258
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:851
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
Definition fespace.cpp:3513
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition fespace.hpp:924
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.hpp:750
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:809
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition fespace.hpp:845
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
Definition fespace.cpp:1767
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:4045
virtual const Operator * GetProlongationMatrix() const
Definition fespace.hpp:727
virtual const Operator * GetRestrictionOperator() const
An abstract operator that performs the same action as GetRestrictionMatrix.
Definition fespace.hpp:746
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3835
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:891
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition fespace.cpp:1510
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const Table *coarse_elem_fos, const DenseTensor localP[]) const
Definition fespace.cpp:1642
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:790
const SparseMatrix * GetConformingProlongation() const
Definition fespace.cpp:1456
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:878
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:848
const Table * GetElementToFaceOrientationTable() const
Definition fespace.hpp:1283
int GetVDim() const
Returns the vector dimension of the finite element space.
Definition fespace.hpp:841
virtual const SparseMatrix * GetHpRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.hpp:754
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
Definition fespace.cpp:3871
static int DecodeDof(int dof)
Helper to return the DOF associated with a sign encoded DOF.
Definition fespace.hpp:1168
Abstract class for all finite elements.
Definition fe_base.hpp:244
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
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition fe_base.hpp:338
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
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition fe_base.hpp:360
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:400
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:331
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...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition fe_base.cpp:182
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Definition mesh.hpp:2937
Vector detJ
Determinants of the Jacobians at all quadrature points.
Definition mesh.hpp:2982
static const int NumGeom
Definition geom.hpp:46
bool Parallel() const
Definition transfer.hpp:50
FiniteElementSpace & dom_fes
Domain FE space.
Definition transfer.hpp:32
FiniteElementSpace & ran_fes
Range FE space.
Definition transfer.hpp:33
MemoryType d_mt
Definition transfer.hpp:45
Operator::Type oper_type
Desired Operator::Type for the construction of all operators defined by the underlying transfer algor...
Definition transfer.hpp:38
const Operator & MakeTrueOperator(FiniteElementSpace &fes_in, FiniteElementSpace &fes_out, const Operator &oper, OperatorHandle &t_oper)
Definition transfer.cpp:35
GridTransfer(FiniteElementSpace &dom_fes_, FiniteElementSpace &ran_fes_)
Definition transfer.cpp:20
The BoomerAMG solver in hypre.
Definition hypre.hpp:1717
void SetPrintLevel(int print_level)
Definition hypre.hpp:1797
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:408
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:2040
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition hypre.cpp:1730
IsoparametricTransformation Transf
Definition eltrans.hpp:733
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition eltrans.cpp:587
Class for integration point with weight.
Definition intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
const Array< real_t > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition intrules.cpp:86
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
BilinearFormIntegrator * mass_integ
Ownership depends on own_mass_integ.
Definition transfer.hpp:139
OperatorHandle F
Forward, coarse-to-fine, operator.
Definition transfer.hpp:142
bool own_mass_integ
Ownership flag for mass_integ.
Definition transfer.hpp:140
const Operator & BackwardOperator() override
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
Definition transfer.cpp:190
OperatorHandle B
Backward, fine-to-coarse, operator.
Definition transfer.hpp:143
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,...
Definition transfer.cpp:147
const Operator & ForwardOperator() override
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
Definition transfer.cpp:156
A standard isoparametric element transformation.
Definition eltrans.hpp:629
void SetPointMat(const DenseMatrix &pm)
Set the underlying point matrix describing the transformation.
Definition eltrans.hpp:668
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition eltrans.cpp:417
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:174
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 ...
H1SpaceLumpedMassOperator(const FiniteElementSpace *fes_ho_, const FiniteElementSpace *fes_lor_, Vector &ML_inv_)
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
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 ...
H1SpaceMixedMassOperator(const FiniteElementSpace *fes_ho_, const FiniteElementSpace *fes_lor_, Table *ho2lor_, Vector *M_LH_ea_)
void Mult(const Vector &x, Vector &y) const override
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...
std::unique_ptr< FiniteElementSpace > fes_ho_scalar
Definition transfer.hpp:475
std::unique_ptr< ParFiniteElementSpace > pfes_ho_scalar
Definition transfer.hpp:483
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.
void SetupPCG()
Sets up the PCG solver (sets parameters, operator, and preconditioner)
void Prolongate(const Vector &x, Vector &y) const override
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.
void SetFromTDofs(const FiniteElementSpace &fes, const Vector &X, Vector &x) const
Sets dof values given a vector of tdofs and a finite element space.
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...
L2ProjectionH1Space(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, const bool use_ea_, MemoryType d_mt_=Device::GetHostMemoryType())
std::unique_ptr< SparseMatrix > AllocR()
Computes sparsity pattern and initializes R matrix. Based on BilinearForm::AllocMat(),...
void MultTranspose(const Vector &x, Vector &y) const override
void LumpedMassInverse(Vector &ML_inv) const
Returns the inverse of an on-rank lumped mass matrix.
std::unique_ptr< FiniteElementSpace > fes_lor_scalar
Definition transfer.hpp:476
void SetAbsTol(real_t p_atol_) override
Sets absolute tolerance in preconditioned conjugate gradient solver.
void ProlongateTranspose(const Vector &x, Vector &y) const override
std::pair< std::unique_ptr< SparseMatrix >, std::unique_ptr< SparseMatrix > > ComputeSparseRAndM_LH()
Computes on-rank R and M_LH matrices. If true, computes mixed mass and/or inverse lumped mass matrix ...
void SetRelTol(real_t p_rtol_) override
Sets relative tolerance in preconditioned conjugate gradient solver.
std::unique_ptr< ParFiniteElementSpace > pfes_lor_scalar
Definition transfer.hpp:484
void ProlongateTranspose(const Vector &x, Vector &y) const override
Definition transfer.cpp:960
L2ProjectionL2Space(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, const bool use_ea_, MemoryType d_mt_=Device::GetHostMemoryType())
Definition transfer.cpp:477
void EAMultTranspose(const Vector &x, Vector &y) const
Definition transfer.cpp:887
void EAProlongate(const Vector &x, Vector &y) const
Definition transfer.cpp:946
void MultTranspose(const Vector &x, Vector &y) const override
Definition transfer.cpp:845
void EAMult(const Vector &x, Vector &y) const
Perform mult on the device (same as above)
Definition transfer.cpp:831
void EAProlongateTranspose(const Vector &x, Vector &y) const
void Prolongate(const Vector &x, Vector &y) const override
Definition transfer.cpp:901
void Mult(const Vector &x, Vector &y) const override
Definition transfer.cpp:792
void BuildHo2Lor(int nel_ho, int nel_lor, const CoarseFineTransformations &cf_tr)
Definition transfer.cpp:239
void MixedMassEA(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, Vector &M_LH, MemoryType d_mt_=Device::GetHostMemoryType())
Definition transfer.cpp:326
L2Projection(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, MemoryType d_mt_=Device::GetHostMemoryType())
Definition transfer.cpp:232
void ElemMixedMass(Geometry::Type geom, const FiniteElement &fe_ho, const FiniteElement &fe_lor, ElementTransformation *tr_ho, ElementTransformation *tr_lor, IntegrationPointTransformation &ip_tr, DenseMatrix &M_mixed_el) const
Definition transfer.cpp:259
L2Projection * F
Forward, coarse-to-fine, operator.
Definition transfer.hpp:511
const Operator & BackwardOperator() override
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
bool SupportsBackwardsOperator() const override
L2Prolongation * B
Backward, fine-to-coarse, operator.
Definition transfer.hpp:512
const Operator & ForwardOperator() override
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat) override
void AssembleEA(const FiniteElementSpace &fes, Vector &emat, const bool add) override
Method defining element assembly.
Class used by MFEM to store pointers to host and/or device memory.
void SyncAlias(const Memory &base, int alias_size) const
Update the alias Memory *this to match the memory location (all valid locations) of its base Memory,...
void CopyFrom(const Memory &src, int size)
Copy size entries from src to *this.
List of mesh geometries stored as Array<Geometry::Type>.
Definition mesh.hpp:1489
Mesh data type.
Definition mesh.hpp:64
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:7254
const CoarseFineTransformations & GetRefinementTransforms() const
Definition mesh.cpp:11407
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
const GeometricFactors * GetGeometricFactors(const IntegrationRule &ir, const int flags, MemoryType d_mt=MemoryType::DEFAULT)
Return the mesh geometric factors corresponding to the given integration rule.
Definition mesh.cpp:880
Geometry::Type GetElementBaseGeometry(int i) const
Definition mesh.hpp:1455
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition handle.hpp:120
Operator * Ptr() const
Access the underlying Operator pointer.
Definition handle.hpp:87
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition handle.hpp:145
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:333
Abstract operator.
Definition operator.hpp:25
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
@ ANY_TYPE
ID for the base class Operator, i.e. any type.
Definition operator.hpp:285
@ MFEM_SPARSEMAT
ID for class SparseMatrix.
Definition operator.hpp:286
@ Hypre_ParCSR
ID for class HypreParMatrix.
Definition operator.hpp:287
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
Matrix-free transfer operator between finite element spaces on the same mesh.
Definition transfer.hpp:567
PRefinementTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator from lFESpace to hFESpace which have different FE collections.
void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
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...
Class for parallel bilinear form.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:334
HYPRE_BigInt * GetTrueDofOffsets() const
Definition pfespace.hpp:343
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition pfespace.hpp:388
ParMesh * GetParMesh() const
Definition pfespace.hpp:338
General product operator: x -> (A*B)(x) = A(B(x)).
Definition operator.hpp:894
Data type sparse matrix.
Definition sparsemat.hpp:51
void MultTranspose(const Vector &x, Vector &y) const override
Multiply a vector with the transposed matrix. y = At * x.
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
void LoseData()
Call this if data has been stolen.
Definition table.hpp:188
int * GetJ()
Definition table.hpp:122
int RowSize(int i) const
Definition table.hpp:116
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition table.cpp:259
int * GetI()
Definition table.hpp:121
void SortRows()
Sort the column (TYPE II) indices in each row.
Definition table.cpp:271
void SetDims(int rows, int nnz)
Definition table.cpp:212
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:1273
Matrix-free transfer operator between finite element spaces on the same mesh exploiting the tensor pr...
Definition transfer.hpp:598
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...
TensorProductPRefinementTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator from lFESpace to hFESpace which have different FE collections.
void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
Matrix-free transfer operator between finite element spaces.
Definition transfer.hpp:536
virtual ~TransferOperator()
Destructor.
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...
void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
TransferOperator(const FiniteElementSpace &lFESpace, const FiniteElementSpace &hFESpace)
Constructs a transfer operator from lFESpace to hFESpace.
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition operator.hpp:847
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition operator.hpp:958
~TrueTransferOperator()
Destructor.
void MultTranspose(const Vector &x, Vector &y) const override
Restriction by applying the transpose of the Mult method.
TrueTransferOperator(const FiniteElementSpace &lFESpace_, const FiniteElementSpace &hFESpace_)
Constructs a transfer operator working on true degrees of freedom from lFESpace to hFESpace.
void Mult(const Vector &x, Vector &y) const override
Interpolation or prolongation of a true dof vector x to a true dof vector y.
Vector data type.
Definition vector.hpp:82
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:494
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:679
Memory< real_t > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition vector.hpp:257
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:745
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
Definition vector.cpp:337
void NewMemoryAndSize(const Memory< real_t > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
Definition vector.hpp:607
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:144
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
void Reciprocal()
(*this)(i) = 1.0 / (*this)(i)
Definition vector.cpp:383
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:653
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:502
void Prolongation2D(const int NE, const int D1D, const int Q1D, const Vector &localL, Vector &localH, const Array< real_t > &B, const Vector &mask)
void Restriction3D(const int NE, const int D1D, const int Q1D, const Vector &localH, Vector &localL, const Array< real_t > &Bt, const Vector &mask)
void Restriction2D(const int NE, const int D1D, const int Q1D, const Vector &localH, Vector &localL, const Array< real_t > &Bt, const Vector &mask)
void Prolongation3D(const int NE, const int D1D, const int Q1D, const Vector &localL, Vector &localH, const Array< real_t > &B, const Vector &mask)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:548
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:391
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition table.cpp:486
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 RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B, bool own_matrix)
Definition hypre.cpp:3007
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition fespace.hpp:83
@ NATIVE
Native ordering as defined by the FiniteElement.
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
void forall(int N, lambda &&body)
Definition forall.hpp:753
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:492
real_t p(const Vector &x, real_t t)
Defines the coarse-fine transformations of all fine elements.
Definition ncmesh.hpp:90
Array< Embedding > embeddings
Fine element positions in their parents.
Definition ncmesh.hpp:92
DenseTensor point_matrices[Geometry::NumGeom]
Definition ncmesh.hpp:96