1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "transfer.hpp"
13 #include "bilinearform.hpp"
14 #include "../general/forall.hpp"
15
16 namespace mfem
17 {
18
GridTransfer(FiniteElementSpace & dom_fes_,FiniteElementSpace & ran_fes_)19 GridTransfer::GridTransfer(FiniteElementSpace &dom_fes_,
20 FiniteElementSpace &ran_fes_)
21 : dom_fes(dom_fes_), ran_fes(ran_fes_),
22 oper_type(Operator::ANY_TYPE),
23 fw_t_oper(), bw_t_oper()
24 {
25 #ifdef MFEM_USE_MPI
26 const bool par_dom = dynamic_cast<ParFiniteElementSpace*>(&dom_fes);
27 const bool par_ran = dynamic_cast<ParFiniteElementSpace*>(&ran_fes);
28 MFEM_VERIFY(par_dom == par_ran, "the domain and range FE spaces must both"
29 " be either serial or parallel");
30 parallel = par_dom;
31 #endif
32 }
33
MakeTrueOperator(FiniteElementSpace & fes_in,FiniteElementSpace & fes_out,const Operator & oper,OperatorHandle & t_oper)34 const Operator &GridTransfer::MakeTrueOperator(
35 FiniteElementSpace &fes_in, FiniteElementSpace &fes_out,
36 const Operator &oper, OperatorHandle &t_oper)
37 {
38 if (t_oper.Ptr())
39 {
40 return *t_oper.Ptr();
41 }
42
43 if (!Parallel())
44 {
45 const SparseMatrix *in_cP = fes_in.GetConformingProlongation();
46 const SparseMatrix *out_cR = fes_out.GetConformingRestriction();
47 if (oper_type == Operator::MFEM_SPARSEMAT)
48 {
49 const SparseMatrix *mat = dynamic_cast<const SparseMatrix *>(&oper);
50 MFEM_VERIFY(mat != NULL, "Operator is not a SparseMatrix");
51 if (!out_cR)
52 {
53 t_oper.Reset(const_cast<SparseMatrix*>(mat), false);
54 }
55 else
56 {
57 t_oper.Reset(mfem::Mult(*out_cR, *mat));
58 }
59 if (in_cP)
60 {
61 t_oper.Reset(mfem::Mult(*t_oper.As<SparseMatrix>(), *in_cP));
62 }
63 }
64 else if (oper_type == Operator::ANY_TYPE)
65 {
66 const int RP_case = bool(out_cR) + 2*bool(in_cP);
67 switch (RP_case)
68 {
69 case 0:
70 t_oper.Reset(const_cast<Operator*>(&oper), false);
71 break;
72 case 1:
73 t_oper.Reset(
74 new ProductOperator(out_cR, &oper, false, false));
75 break;
76 case 2:
77 t_oper.Reset(
78 new ProductOperator(&oper, in_cP, false, false));
79 break;
80 case 3:
81 t_oper.Reset(
82 new TripleProductOperator(
83 out_cR, &oper, in_cP, false, false, false));
84 break;
85 }
86 }
87 else
88 {
89 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
90 }
91 }
92 else // Parallel() == true
93 {
94 #ifdef MFEM_USE_MPI
95 if (oper_type == Operator::Hypre_ParCSR)
96 {
97 const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
98 const ParFiniteElementSpace *pfes_in =
99 dynamic_cast<const ParFiniteElementSpace *>(&fes_in);
100 const ParFiniteElementSpace *pfes_out =
101 dynamic_cast<const ParFiniteElementSpace *>(&fes_out);
102 const SparseMatrix *sp_mat = dynamic_cast<const SparseMatrix *>(&oper);
103 const HypreParMatrix *hy_mat;
104 if (sp_mat)
105 {
106 SparseMatrix *RA = mfem::Mult(*out_R, *sp_mat);
107 t_oper.Reset(pfes_in->Dof_TrueDof_Matrix()->
108 LeftDiagMult(*RA, pfes_out->GetTrueDofOffsets()));
109 delete RA;
110 }
111 else if ((hy_mat = dynamic_cast<const HypreParMatrix *>(&oper)))
112 {
113 HypreParMatrix *RA =
114 hy_mat->LeftDiagMult(*out_R, pfes_out->GetTrueDofOffsets());
115 t_oper.Reset(mfem::ParMult(RA, pfes_in->Dof_TrueDof_Matrix()));
116 delete RA;
117 }
118 else
119 {
120 MFEM_ABORT("unknown Operator type");
121 }
122 }
123 else if (oper_type == Operator::ANY_TYPE)
124 {
125 const Operator *out_R = fes_out.GetRestrictionOperator();
126 t_oper.Reset(new TripleProductOperator(
127 out_R, &oper, fes_in.GetProlongationMatrix(),
128 false, false, false));
129 }
130 else
131 {
132 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
133 }
134 #endif
135 }
136
137 return *t_oper.Ptr();
138 }
139
140
~InterpolationGridTransfer()141 InterpolationGridTransfer::~InterpolationGridTransfer()
142 {
143 if (own_mass_integ) { delete mass_integ; }
144 }
145
SetMassIntegrator(BilinearFormIntegrator * mass_integ_,bool own_mass_integ_)146 void InterpolationGridTransfer::SetMassIntegrator(
147 BilinearFormIntegrator *mass_integ_, bool own_mass_integ_)
148 {
149 if (own_mass_integ) { delete mass_integ; }
150
151 mass_integ = mass_integ_;
152 own_mass_integ = own_mass_integ_;
153 }
154
ForwardOperator()155 const Operator &InterpolationGridTransfer::ForwardOperator()
156 {
157 if (F.Ptr())
158 {
159 return *F.Ptr();
160 }
161
162 // Construct F
163 if (oper_type == Operator::ANY_TYPE)
164 {
165 F.Reset(new FiniteElementSpace::RefinementOperator(&ran_fes, &dom_fes));
166 }
167 else if (oper_type == Operator::MFEM_SPARSEMAT)
168 {
169 Mesh::GeometryList elem_geoms(*ran_fes.GetMesh());
170
171 DenseTensor localP[Geometry::NumGeom];
172 for (int i = 0; i < elem_geoms.Size(); i++)
173 {
174 ran_fes.GetLocalRefinementMatrices(dom_fes, elem_geoms[i],
175 localP[elem_geoms[i]]);
176 }
177 F.Reset(ran_fes.RefinementMatrix_main(
178 dom_fes.GetNDofs(), dom_fes.GetElementToDofTable(), localP));
179 }
180 else
181 {
182 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
183 }
184
185 return *F.Ptr();
186 }
187
BackwardOperator()188 const Operator &InterpolationGridTransfer::BackwardOperator()
189 {
190 if (B.Ptr())
191 {
192 return *B.Ptr();
193 }
194
195 // Construct B, if not set, define a suitable mass_integ
196 if (!mass_integ && ran_fes.GetNE() > 0)
197 {
198 const FiniteElement *f_fe_0 = ran_fes.GetFE(0);
199 const int map_type = f_fe_0->GetMapType();
200 if (map_type == FiniteElement::VALUE ||
201 map_type == FiniteElement::INTEGRAL)
202 {
203 mass_integ = new MassIntegrator;
204 }
205 else if (map_type == FiniteElement::H_DIV ||
206 map_type == FiniteElement::H_CURL)
207 {
208 mass_integ = new VectorFEMassIntegrator;
209 }
210 else
211 {
212 MFEM_ABORT("unknown type of FE space");
213 }
214 own_mass_integ = true;
215 }
216 if (oper_type == Operator::ANY_TYPE)
217 {
218 B.Reset(new FiniteElementSpace::DerefinementOperator(
219 &ran_fes, &dom_fes, mass_integ));
220 }
221 else
222 {
223 MFEM_ABORT("Operator::Type is not supported: " << oper_type);
224 }
225
226 return *B.Ptr();
227 }
228
229
L2Projection(const FiniteElementSpace & fes_ho_,const FiniteElementSpace & fes_lor_)230 L2ProjectionGridTransfer::L2Projection::L2Projection(
231 const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
232 : Operator(fes_lor_.GetVSize(), fes_ho_.GetVSize()),
233 fes_ho(fes_ho_),
234 fes_lor(fes_lor_)
235 { }
236
BuildHo2Lor(int nel_ho,int nel_lor,const CoarseFineTransformations & cf_tr)237 void L2ProjectionGridTransfer::L2Projection::BuildHo2Lor(
238 int nel_ho, int nel_lor, const CoarseFineTransformations& cf_tr)
239 {
240 // Construct the mapping from HO to LOR
241 // ho2lor.GetRow(iho) will give all the LOR elements contained in iho
242 ho2lor.MakeI(nel_ho);
243 for (int ilor = 0; ilor < nel_lor; ++ilor)
244 {
245 int iho = cf_tr.embeddings[ilor].parent;
246 ho2lor.AddAColumnInRow(iho);
247 }
248 ho2lor.MakeJ();
249 for (int ilor = 0; ilor < nel_lor; ++ilor)
250 {
251 int iho = cf_tr.embeddings[ilor].parent;
252 ho2lor.AddConnection(iho, ilor);
253 }
254 ho2lor.ShiftUpI();
255 }
256
ElemMixedMass(Geometry::Type geom,const FiniteElement & fe_ho,const FiniteElement & fe_lor,ElementTransformation * el_tr,IntegrationPointTransformation & ip_tr,DenseMatrix & M_mixed_el) const257 void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
258 Geometry::Type geom, const FiniteElement& fe_ho,
259 const FiniteElement& fe_lor, ElementTransformation* el_tr,
260 IntegrationPointTransformation& ip_tr,
261 DenseMatrix& M_mixed_el) const
262 {
263 int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
264 const IntegrationRule* ir = &IntRules.Get(geom, order);
265 M_mixed_el = 0.0;
266 for (int i = 0; i < ir->GetNPoints(); i++)
267 {
268 const IntegrationPoint& ip_lor = ir->IntPoint(i);
269 IntegrationPoint ip_ho;
270 ip_tr.Transform(ip_lor, ip_ho);
271 Vector shape_lor(fe_lor.GetDof());
272 fe_lor.CalcShape(ip_lor, shape_lor);
273 Vector shape_ho(fe_ho.GetDof());
274 fe_ho.CalcShape(ip_ho, shape_ho);
275 el_tr->SetIntPoint(&ip_lor);
276 // For now we use the geometry information from the LOR space, which means
277 // we won't be mass conservative if the mesh is curved
278 double w = el_tr->Weight() * ip_lor.weight;
279 shape_lor *= w;
280 AddMultVWt(shape_lor, shape_ho, M_mixed_el);
281 }
282 }
283
L2ProjectionL2Space(const FiniteElementSpace & fes_ho_,const FiniteElementSpace & fes_lor_)284 L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space(
285 const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
286 : L2Projection(fes_ho_, fes_lor_)
287 {
288 Mesh *mesh_ho = fes_ho.GetMesh();
289 Mesh *mesh_lor = fes_lor.GetMesh();
290 int nel_ho = mesh_ho->GetNE();
291 int nel_lor = mesh_lor->GetNE();
292
293 // If the local mesh is empty, skip all computations
294 if (nel_ho == 0) { return; }
295
296 const CoarseFineTransformations &cf_tr = mesh_lor->GetRefinementTransforms();
297
298 int nref_max = 0;
299 Array<Geometry::Type> geoms;
300 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
301 for (int ig = 0; ig < geoms.Size(); ++ig)
302 {
303 Geometry::Type geom = geoms[ig];
304 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
305 }
306
307 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
308
309 offsets.SetSize(nel_ho+1);
310 offsets[0] = 0;
311 for (int iho = 0; iho < nel_ho; ++iho)
312 {
313 int nref = ho2lor.RowSize(iho);
314 const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
315 const FiniteElement &fe_lor = *fes_lor.GetFE(ho2lor.GetRow(iho)[0]);
316 offsets[iho+1] = offsets[iho] + fe_ho.GetDof()*fe_lor.GetDof()*nref;
317 }
318 // R will contain the restriction (L^2 projection operator) defined on each
319 // coarse HO element (and corresponding patch of LOR elements)
320 R.SetSize(offsets[nel_ho]);
321 // P will contain the corresponding prolongation operator
322 P.SetSize(offsets[nel_ho]);
323
324 IntegrationPointTransformation ip_tr;
325 IsoparametricTransformation &emb_tr = ip_tr.Transf;
326
327 for (int iho = 0; iho < nel_ho; ++iho)
328 {
329 Array<int> lor_els;
330 ho2lor.GetRow(iho, lor_els);
331 int nref = ho2lor.RowSize(iho);
332
333 Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
334 const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
335 const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
336 int ndof_ho = fe_ho.GetDof();
337 int ndof_lor = fe_lor.GetDof();
338
339 emb_tr.SetIdentityTransformation(geom);
340 const DenseTensor &pmats = cf_tr.point_matrices[geom];
341
342 DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
343 DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
344
345 DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
346 DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
347
348 MassIntegrator mi;
349 DenseMatrix M_lor_el(ndof_lor, ndof_lor);
350 DenseMatrixInverse Minv_lor_el(&M_lor_el);
351 DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
352 DenseMatrix M_mixed_el(ndof_lor, ndof_ho);
353
354 Minv_lor = 0.0;
355 M_lor = 0.0;
356
357 DenseMatrix RtMlor(ndof_ho, ndof_lor*nref);
358 DenseMatrix RtMlorR(ndof_ho, ndof_ho);
359 DenseMatrixInverse RtMlorR_inv(&RtMlorR);
360
361 for (int iref = 0; iref < nref; ++iref)
362 {
363 // Assemble the low-order refined mass matrix and invert locally
364 int ilor = lor_els[iref];
365 ElementTransformation *el_tr = fes_lor.GetElementTransformation(ilor);
366 mi.AssembleElementMatrix(fe_lor, *el_tr, M_lor_el);
367 M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
368 Minv_lor_el.Factor();
369 Minv_lor_el.GetInverseMatrix(M_lor_el);
370 // Insert into the diagonal of the patch LOR mass matrix
371 Minv_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
372
373 // Now assemble the block-row of the mixed mass matrix associated
374 // with integrating HO functions against LOR functions on the LOR
375 // sub-element.
376
377 // Create the transformation that embeds the fine low-order element
378 // within the coarse high-order element in reference space
379 emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
380
381 ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, M_mixed_el);
382
383 M_mixed.CopyMN(M_mixed_el, iref*ndof_lor, 0);
384 }
385 mfem::Mult(Minv_lor, M_mixed, R_iho);
386
387 mfem::MultAtB(R_iho, M_lor, RtMlor);
388 mfem::Mult(RtMlor, R_iho, RtMlorR);
389 RtMlorR_inv.Factor();
390 RtMlorR_inv.Mult(RtMlor, P_iho);
391 }
392 }
393
Mult(const Vector & x,Vector & y) const394 void L2ProjectionGridTransfer::L2ProjectionL2Space::Mult(
395 const Vector &x, Vector &y) const
396 {
397 int vdim = fes_ho.GetVDim();
398 Array<int> vdofs;
399 DenseMatrix xel_mat, yel_mat;
400 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
401 {
402 int nref = ho2lor.RowSize(iho);
403 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
404 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
405 xel_mat.SetSize(ndof_ho, vdim);
406 yel_mat.SetSize(ndof_lor*nref, vdim);
407 DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
408
409 fes_ho.GetElementVDofs(iho, vdofs);
410 x.GetSubVector(vdofs, xel_mat.GetData());
411 mfem::Mult(R_iho, xel_mat, yel_mat);
412 // Place result correctly into the low-order vector
413 for (int iref = 0; iref < nref; ++iref)
414 {
415 int ilor = ho2lor.GetRow(iho)[iref];
416 for (int vd=0; vd<vdim; ++vd)
417 {
418 fes_lor.GetElementDofs(ilor, vdofs);
419 fes_lor.DofsToVDofs(vd, vdofs);
420 y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
421 }
422 }
423 }
424 }
425
MultTranspose(const Vector & x,Vector & y) const426 void L2ProjectionGridTransfer::L2ProjectionL2Space::MultTranspose(
427 const Vector &x, Vector &y) const
428 {
429 int vdim = fes_ho.GetVDim();
430 Array<int> vdofs;
431 DenseMatrix xel_mat, yel_mat;
432 y = 0.0;
433 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
434 {
435 int nref = ho2lor.RowSize(iho);
436 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
437 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
438 xel_mat.SetSize(ndof_lor*nref, vdim);
439 yel_mat.SetSize(ndof_ho, vdim);
440 DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
441
442 // Extract the LOR DOFs
443 for (int iref=0; iref<nref; ++iref)
444 {
445 int ilor = ho2lor.GetRow(iho)[iref];
446 for (int vd=0; vd<vdim; ++vd)
447 {
448 fes_lor.GetElementDofs(ilor, vdofs);
449 fes_lor.DofsToVDofs(vd, vdofs);
450 x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
451 }
452 }
453 // Multiply locally by the transpose
454 mfem::MultAtB(R_iho, xel_mat, yel_mat);
455 // Place the result in the HO vector
456 fes_ho.GetElementVDofs(iho, vdofs);
457 y.AddElementVector(vdofs, yel_mat.GetData());
458 }
459 }
460
Prolongate(const Vector & x,Vector & y) const461 void L2ProjectionGridTransfer::L2ProjectionL2Space::Prolongate(
462 const Vector &x, Vector &y) const
463 {
464 int vdim = fes_ho.GetVDim();
465 Array<int> vdofs;
466 DenseMatrix xel_mat,yel_mat;
467 y = 0.0;
468 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
469 {
470 int nref = ho2lor.RowSize(iho);
471 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
472 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
473 xel_mat.SetSize(ndof_lor*nref, vdim);
474 yel_mat.SetSize(ndof_ho, vdim);
475 DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
476
477 // Extract the LOR DOFs
478 for (int iref = 0; iref < nref; ++iref)
479 {
480 int ilor = ho2lor.GetRow(iho)[iref];
481 for (int vd = 0; vd < vdim; ++vd)
482 {
483 fes_lor.GetElementDofs(ilor, vdofs);
484 fes_lor.DofsToVDofs(vd, vdofs);
485 x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
486 }
487 }
488 // Locally prolongate
489 mfem::Mult(P_iho, xel_mat, yel_mat);
490 // Place the result in the HO vector
491 fes_ho.GetElementVDofs(iho, vdofs);
492 y.AddElementVector(vdofs, yel_mat.GetData());
493 }
494 }
495
ProlongateTranspose(const Vector & x,Vector & y) const496 void L2ProjectionGridTransfer::L2ProjectionL2Space::ProlongateTranspose(
497 const Vector &x, Vector &y) const
498 {
499 int vdim = fes_ho.GetVDim();
500 Array<int> vdofs;
501 DenseMatrix xel_mat,yel_mat;
502 for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
503 {
504 int nref = ho2lor.RowSize(iho);
505 int ndof_ho = fes_ho.GetFE(iho)->GetDof();
506 int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
507 xel_mat.SetSize(ndof_ho, vdim);
508 yel_mat.SetSize(ndof_lor*nref, vdim);
509 DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
510
511 fes_ho.GetElementVDofs(iho, vdofs);
512 x.GetSubVector(vdofs, xel_mat.GetData());
513 mfem::MultAtB(P_iho, xel_mat, yel_mat);
514
515 // Place result correctly into the low-order vector
516 for (int iref = 0; iref < nref; ++iref)
517 {
518 int ilor = ho2lor.GetRow(iho)[iref];
519 for (int vd=0; vd<vdim; ++vd)
520 {
521 fes_lor.GetElementDofs(ilor, vdofs);
522 fes_lor.DofsToVDofs(vd, vdofs);
523 y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
524 }
525 }
526 }
527 }
528
L2ProjectionH1Space(const FiniteElementSpace & fes_ho_,const FiniteElementSpace & fes_lor_)529 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
530 const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_)
531 : L2Projection(fes_ho_, fes_lor_)
532 {
533 Mesh* mesh_ho = fes_ho.GetMesh();
534 Mesh* mesh_lor = fes_lor.GetMesh();
535 int nel_ho = mesh_ho->GetNE();
536 int nel_lor = mesh_lor->GetNE();
537 int ndof_lor = fes_lor.GetNDofs();
538
539 // If the local mesh is empty, skip all computations
540 if (nel_ho == 0) { return; }
541
542 const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
543
544 int nref_max = 0;
545 Array<Geometry::Type> geoms;
546 mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
547 for (int ig = 0; ig < geoms.Size(); ++ig)
548 {
549 Geometry::Type geom = geoms[ig];
550 nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
551 }
552
553 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
554
555 // ML_inv contains the inverse lumped (row sum) mass matrix. Note that the
556 // method will also work with a full (consistent) mass matrix, though this is
557 // not implemented here. L refers to the low-order refined mesh
558 Vector ML_inv(ndof_lor);
559 ML_inv = 0.0;
560
561 // Compute ML_inv
562 for (int iho = 0; iho < nel_ho; ++iho)
563 {
564 Array<int> lor_els;
565 ho2lor.GetRow(iho, lor_els);
566 int nref = ho2lor.RowSize(iho);
567
568 Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
569 const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
570 int nedof_lor = fe_lor.GetDof();
571
572 // Instead of using a MassIntegrator, manually loop over integration
573 // points so we can row sum and store the diagonal as a Vector.
574 Vector ML_el(nedof_lor);
575 Vector shape_lor(nedof_lor);
576 Array<int> dofs_lor(nedof_lor);
577
578 for (int iref = 0; iref < nref; ++iref)
579 {
580 int ilor = lor_els[iref];
581 ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
582
583 int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
584 const IntegrationRule* ir = &IntRules.Get(geom, order);
585 ML_el = 0.0;
586 for (int i = 0; i < ir->GetNPoints(); ++i)
587 {
588 const IntegrationPoint& ip_lor = ir->IntPoint(i);
589 fe_lor.CalcShape(ip_lor, shape_lor);
590 el_tr->SetIntPoint(&ip_lor);
591 ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
592 }
593 fes_lor.GetElementDofs(ilor, dofs_lor);
594 ML_inv.AddElementVector(dofs_lor, ML_el);
595 }
596 }
597 // DOF by DOF inverse of non-zero entries
598 for (int i = 0; i < ndof_lor; ++i)
599 {
600 ML_inv[i] = 1.0 / ML_inv[i];
601 }
602
603 // Compute sparsity pattern for R = M_L^(-1) M_LH and allocate
604 AllocR();
605 // Allocate M_LH (same sparsity pattern as R)
606 // L refers to the low-order refined mesh (DOFs correspond to rows)
607 // H refers to the higher-order mesh (DOFs correspond to columns)
608 M_LH = SparseMatrix(R.GetI(), R.GetJ(), NULL,
609 R.Height(), R.Width(), false, true, true);
610
611 IntegrationPointTransformation ip_tr;
612 IsoparametricTransformation& emb_tr = ip_tr.Transf;
613
614 // Compute M_LH and R
615 for (int iho = 0; iho < nel_ho; ++iho)
616 {
617 Array<int> lor_els;
618 ho2lor.GetRow(iho, lor_els);
619 int nref = ho2lor.RowSize(iho);
620
621 Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
622 const FiniteElement& fe_ho = *fes_ho.GetFE(iho);
623 const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
624
625 emb_tr.SetIdentityTransformation(geom);
626 const DenseTensor& pmats = cf_tr.point_matrices[geom];
627
628 int nedof_ho = fe_ho.GetDof();
629 int nedof_lor = fe_lor.GetDof();
630 DenseMatrix M_LH_el(nedof_lor, nedof_ho);
631 DenseMatrix R_el(nedof_lor, nedof_ho);
632
633 for (int iref = 0; iref < nref; ++iref)
634 {
635 int ilor = lor_els[iref];
636 ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
637
638 // Create the transformation that embeds the fine low-order element
639 // within the coarse high-order element in reference space
640 emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
641
642 ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, M_LH_el);
643
644 Array<int> dofs_lor(nedof_lor);
645 fes_lor.GetElementDofs(ilor, dofs_lor);
646 Vector R_row;
647 for (int i = 0; i < nedof_lor; ++i)
648 {
649 M_LH_el.GetRow(i, R_row);
650 R_el.SetRow(i, R_row.Set(ML_inv[dofs_lor[i]], R_row));
651 }
652 Array<int> dofs_ho(nedof_ho);
653 fes_ho.GetElementDofs(iho, dofs_ho);
654 M_LH.AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
655 R.AddSubMatrix(dofs_lor, dofs_ho, R_el);
656 }
657 }
658
659 // Create PCG solver
660 RTxM_LH = TransposeMult(R, M_LH);
661 pcg.SetPrintLevel(0);
662 pcg.SetMaxIter(1000);
663 // initial values for relative and absolute tolerance
664 SetRelTol(1e-13);
665 SetAbsTol(1e-13);
666 Ds = DSmoother(*RTxM_LH);
667 pcg.SetPreconditioner(Ds);
668 pcg.SetOperator(*RTxM_LH);
669 }
670
~L2ProjectionH1Space()671 L2ProjectionGridTransfer::L2ProjectionH1Space::~L2ProjectionH1Space()
672 {
673 delete RTxM_LH;
674 }
675
Mult(const Vector & x,Vector & y) const676 void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult(
677 const Vector& x, Vector& y) const
678 {
679 int vdim = fes_ho.GetVDim();
680 const int ndof_ho = fes_ho.GetNDofs();
681 const int ndof_lor = fes_lor.GetNDofs();
682 Array<int> dofs_ho(ndof_ho);
683 Array<int> dofs_lor(ndof_lor);
684 Vector x_dim(ndof_ho);
685 Vector y_dim(ndof_lor);
686
687 for (int d = 0; d < vdim; ++d)
688 {
689 fes_ho.GetVDofs(d, dofs_ho);
690 fes_lor.GetVDofs(d, dofs_lor);
691 x.GetSubVector(dofs_ho, x_dim);
692 R.Mult(x_dim, y_dim);
693 y.SetSubVector(dofs_lor, y_dim);
694 }
695 }
696
MultTranspose(const Vector & x,Vector & y) const697 void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose(
698 const Vector& x, Vector& y) const
699 {
700 int vdim = fes_ho.GetVDim();
701 const int ndof_ho = fes_ho.GetNDofs();
702 const int ndof_lor = fes_lor.GetNDofs();
703 Array<int> dofs_ho(ndof_ho);
704 Array<int> dofs_lor(ndof_lor);
705 Vector x_dim(ndof_lor);
706 Vector y_dim(ndof_ho);
707
708 for (int d = 0; d < vdim; ++d)
709 {
710 fes_ho.GetVDofs(d, dofs_ho);
711 fes_lor.GetVDofs(d, dofs_lor);
712 x.GetSubVector(dofs_lor, x_dim);
713 R.MultTranspose(x_dim, y_dim);
714 y.SetSubVector(dofs_ho, y_dim);
715 }
716 }
717
Prolongate(const Vector & x,Vector & y) const718 void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
719 const Vector& x, Vector& y) const
720 {
721 int vdim = fes_ho.GetVDim();
722 const int ndof_ho = fes_ho.GetNDofs();
723 const int ndof_lor = fes_lor.GetNDofs();
724 Array<int> dofs_ho(ndof_ho);
725 Array<int> dofs_lor(ndof_lor);
726 Vector x_dim(ndof_lor);
727 Vector y_dim(ndof_ho);
728 Vector xbar(ndof_ho);
729
730 for (int d = 0; d < vdim; ++d)
731 {
732 fes_lor.GetVDofs(d, dofs_lor);
733 x.GetSubVector(dofs_lor, x_dim);
734 // Compute y = P x = (R^T M_LH)^(-1) M_LH^T x = (R^T M_LH)^(-1) xbar
735 M_LH.MultTranspose(x_dim, xbar);
736 y_dim = 0.0;
737 pcg.Mult(xbar, y_dim);
738 fes_ho.GetVDofs(d, dofs_ho);
739 y.SetSubVector(dofs_ho, y_dim);
740 }
741 }
742
ProlongateTranspose(const Vector & x,Vector & y) const743 void L2ProjectionGridTransfer::L2ProjectionH1Space::ProlongateTranspose(
744 const Vector& x, Vector& y) const
745 {
746 int vdim = fes_ho.GetVDim();
747 const int ndof_ho = fes_ho.GetNDofs();
748 const int ndof_lor = fes_lor.GetNDofs();
749 Array<int> dofs_ho(ndof_ho);
750 Array<int> dofs_lor(ndof_lor);
751 Vector x_dim(ndof_ho);
752 Vector y_dim(ndof_lor);
753 Vector xbar(ndof_ho);
754
755 for (int d = 0; d < vdim; ++d)
756 {
757 fes_ho.GetVDofs(d, dofs_ho);
758 x.GetSubVector(dofs_ho, x_dim);
759 // Compute y = P^T x = M_LH (R^T M_LH)^(-1) x = M_LH xbar
760 xbar = 0.0;
761 pcg.Mult(x_dim, xbar);
762 M_LH.Mult(xbar, y_dim);
763 fes_lor.GetVDofs(d, dofs_lor);
764 y.SetSubVector(dofs_lor, y_dim);
765 }
766 }
767
SetRelTol(double p_rtol_)768 void L2ProjectionGridTransfer::L2ProjectionH1Space::SetRelTol(double p_rtol_)
769 {
770 pcg.SetRelTol(p_rtol_);
771 }
772
SetAbsTol(double p_atol_)773 void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(double p_atol_)
774 {
775 pcg.SetAbsTol(p_atol_);
776 }
777
AllocR()778 void L2ProjectionGridTransfer::L2ProjectionH1Space::AllocR()
779 {
780 const Table& elem_dof_ho = fes_ho.GetElementToDofTable();
781 const Table& elem_dof_lor = fes_lor.GetElementToDofTable();
782 const int ndof_ho = fes_ho.GetNDofs();
783 const int ndof_lor = fes_lor.GetNDofs();
784
785 Table dof_elem_lor;
786 Transpose(elem_dof_lor, dof_elem_lor, ndof_lor);
787
788 Mesh* mesh_lor = fes_lor.GetMesh();
789 const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
790
791 // mfem::Mult but uses ho2lor to map HO elements to LOR elements
792 const int* elem_dof_hoI = elem_dof_ho.GetI();
793 const int* elem_dof_hoJ = elem_dof_ho.GetJ();
794 const int* dof_elem_lorI = dof_elem_lor.GetI();
795 const int* dof_elem_lorJ = dof_elem_lor.GetJ();
796
797 Array<int> I(ndof_lor + 1);
798
799 // figure out the size of J
800 Array<int> dof_used_ho;
801 dof_used_ho.SetSize(ndof_ho, -1);
802
803 int sizeJ = 0;
804 for (int ilor = 0; ilor < ndof_lor; ++ilor)
805 {
806 for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
807 {
808 int el_lor = dof_elem_lorJ[jlor];
809 int iho = cf_tr.embeddings[el_lor].parent;
810 for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
811 {
812 int dof_ho = elem_dof_hoJ[jho];
813 if (dof_used_ho[dof_ho] != ilor)
814 {
815 dof_used_ho[dof_ho] = ilor;
816 ++sizeJ;
817 }
818 }
819 }
820 }
821
822 // initialize dof_ho_dof_lor
823 Table dof_lor_dof_ho;
824 dof_lor_dof_ho.SetDims(ndof_lor, sizeJ);
825
826 for (int i = 0; i < ndof_ho; ++i)
827 {
828 dof_used_ho[i] = -1;
829 }
830
831 // set values of J
832 int* dof_dofI = dof_lor_dof_ho.GetI();
833 int* dof_dofJ = dof_lor_dof_ho.GetJ();
834 sizeJ = 0;
835 for (int ilor = 0; ilor < ndof_lor; ++ilor)
836 {
837 dof_dofI[ilor] = sizeJ;
838 for (int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
839 {
840 int el_lor = dof_elem_lorJ[jlor];
841 int iho = cf_tr.embeddings[el_lor].parent;
842 for (int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
843 {
844 int dof_ho = elem_dof_hoJ[jho];
845 if (dof_used_ho[dof_ho] != ilor)
846 {
847 dof_used_ho[dof_ho] = ilor;
848 dof_dofJ[sizeJ] = dof_ho;
849 ++sizeJ;
850 }
851 }
852 }
853 }
854
855 dof_lor_dof_ho.SortRows();
856 double* data = Memory<double>(dof_dofI[ndof_lor]);
857
858 R = SparseMatrix(dof_dofI, dof_dofJ, data, ndof_lor, ndof_ho,
859 true, true, true);
860 R = 0.0;
861
862 dof_lor_dof_ho.LoseData();
863 }
864
~L2ProjectionGridTransfer()865 L2ProjectionGridTransfer::~L2ProjectionGridTransfer()
866 {
867 delete F;
868 delete B;
869 }
870
ForwardOperator()871 const Operator &L2ProjectionGridTransfer::ForwardOperator()
872 {
873 if (!F) { BuildF(); }
874 return *F;
875 }
876
BackwardOperator()877 const Operator &L2ProjectionGridTransfer::BackwardOperator()
878 {
879 if (!B)
880 {
881 if (!F) { BuildF(); }
882 B = new L2Prolongation(*F);
883 }
884 return *B;
885 }
886
BuildF()887 void L2ProjectionGridTransfer::BuildF()
888 {
889 if (!force_l2_space &&
890 dom_fes.FEColl()->GetContType() == FiniteElementCollection::CONTINUOUS)
891 {
892 F = new L2ProjectionH1Space(dom_fes, ran_fes);
893 }
894 else
895 {
896 F = new L2ProjectionL2Space(dom_fes, ran_fes);
897 }
898 }
899
900
TransferOperator(const FiniteElementSpace & lFESpace_,const FiniteElementSpace & hFESpace_)901 TransferOperator::TransferOperator(const FiniteElementSpace& lFESpace_,
902 const FiniteElementSpace& hFESpace_)
903 : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize())
904 {
905 if (lFESpace_.FEColl() == hFESpace_.FEColl())
906 {
907 OperatorPtr P(Operator::ANY_TYPE);
908 hFESpace_.GetTransferOperator(lFESpace_, P);
909 P.SetOperatorOwner(false);
910 opr = P.Ptr();
911 }
912 else if (lFESpace_.GetMesh()->GetNE() > 0
913 && hFESpace_.GetMesh()->GetNE() > 0
914 && dynamic_cast<const TensorBasisElement*>(lFESpace_.GetFE(0))
915 && dynamic_cast<const TensorBasisElement*>(hFESpace_.GetFE(0)))
916 {
917 opr = new TensorProductPRefinementTransferOperator(lFESpace_, hFESpace_);
918 }
919 else
920 {
921 opr = new PRefinementTransferOperator(lFESpace_, hFESpace_);
922 }
923 }
924
~TransferOperator()925 TransferOperator::~TransferOperator() { delete opr; }
926
Mult(const Vector & x,Vector & y) const927 void TransferOperator::Mult(const Vector& x, Vector& y) const
928 {
929 opr->Mult(x, y);
930 }
931
MultTranspose(const Vector & x,Vector & y) const932 void TransferOperator::MultTranspose(const Vector& x, Vector& y) const
933 {
934 opr->MultTranspose(x, y);
935 }
936
937
PRefinementTransferOperator(const FiniteElementSpace & lFESpace_,const FiniteElementSpace & hFESpace_)938 PRefinementTransferOperator::PRefinementTransferOperator(
939 const FiniteElementSpace& lFESpace_, const FiniteElementSpace& hFESpace_)
940 : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
941 hFESpace(hFESpace_)
942 {
943 }
944
~PRefinementTransferOperator()945 PRefinementTransferOperator::~PRefinementTransferOperator() {}
946
Mult(const Vector & x,Vector & y) const947 void PRefinementTransferOperator::Mult(const Vector& x, Vector& y) const
948 {
949 Mesh* mesh = hFESpace.GetMesh();
950 Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
951 DenseMatrix loc_prol;
952 Vector subY, subX;
953
954 Geometry::Type cached_geom = Geometry::INVALID;
955 const FiniteElement* h_fe = NULL;
956 const FiniteElement* l_fe = NULL;
957 IsoparametricTransformation T;
958
959 int vdim = lFESpace.GetVDim();
960
961 for (int i = 0; i < mesh->GetNE(); i++)
962 {
963 hFESpace.GetElementDofs(i, h_dofs);
964 lFESpace.GetElementDofs(i, l_dofs);
965
966 const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
967 if (geom != cached_geom)
968 {
969 h_fe = hFESpace.GetFE(i);
970 l_fe = lFESpace.GetFE(i);
971 T.SetIdentityTransformation(h_fe->GetGeomType());
972 h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
973 subY.SetSize(loc_prol.Height());
974 cached_geom = geom;
975 }
976
977 for (int vd = 0; vd < vdim; vd++)
978 {
979 l_dofs.Copy(l_vdofs);
980 lFESpace.DofsToVDofs(vd, l_vdofs);
981 h_dofs.Copy(h_vdofs);
982 hFESpace.DofsToVDofs(vd, h_vdofs);
983 x.GetSubVector(l_vdofs, subX);
984 loc_prol.Mult(subX, subY);
985 y.SetSubVector(h_vdofs, subY);
986 }
987 }
988 }
989
MultTranspose(const Vector & x,Vector & y) const990 void PRefinementTransferOperator::MultTranspose(const Vector& x,
991 Vector& y) const
992 {
993 y = 0.0;
994
995 Mesh* mesh = hFESpace.GetMesh();
996 Array<int> l_dofs, h_dofs, l_vdofs, h_vdofs;
997 DenseMatrix loc_prol;
998 Vector subY, subX;
999
1000 Array<char> processed(hFESpace.GetVSize());
1001 processed = 0;
1002
1003 Geometry::Type cached_geom = Geometry::INVALID;
1004 const FiniteElement* h_fe = NULL;
1005 const FiniteElement* l_fe = NULL;
1006 IsoparametricTransformation T;
1007
1008 int vdim = lFESpace.GetVDim();
1009
1010 for (int i = 0; i < mesh->GetNE(); i++)
1011 {
1012 hFESpace.GetElementDofs(i, h_dofs);
1013 lFESpace.GetElementDofs(i, l_dofs);
1014
1015 const Geometry::Type geom = mesh->GetElementBaseGeometry(i);
1016 if (geom != cached_geom)
1017 {
1018 h_fe = hFESpace.GetFE(i);
1019 l_fe = lFESpace.GetFE(i);
1020 T.SetIdentityTransformation(h_fe->GetGeomType());
1021 h_fe->GetTransferMatrix(*l_fe, T, loc_prol);
1022 loc_prol.Transpose();
1023 subY.SetSize(loc_prol.Height());
1024 cached_geom = geom;
1025 }
1026
1027 for (int vd = 0; vd < vdim; vd++)
1028 {
1029 l_dofs.Copy(l_vdofs);
1030 lFESpace.DofsToVDofs(vd, l_vdofs);
1031 h_dofs.Copy(h_vdofs);
1032 hFESpace.DofsToVDofs(vd, h_vdofs);
1033
1034 x.GetSubVector(h_vdofs, subX);
1035 for (int p = 0; p < h_dofs.Size(); ++p)
1036 {
1037 if (processed[lFESpace.DecodeDof(h_dofs[p])])
1038 {
1039 subX[p] = 0.0;
1040 }
1041 }
1042
1043 loc_prol.Mult(subX, subY);
1044 y.AddElementVector(l_vdofs, subY);
1045 }
1046
1047 for (int p = 0; p < h_dofs.Size(); ++p)
1048 {
1049 processed[lFESpace.DecodeDof(h_dofs[p])] = 1;
1050 }
1051 }
1052 }
1053
1054
1055 TensorProductPRefinementTransferOperator::
TensorProductPRefinementTransferOperator(const FiniteElementSpace & lFESpace_,const FiniteElementSpace & hFESpace_)1056 TensorProductPRefinementTransferOperator(
1057 const FiniteElementSpace& lFESpace_,
1058 const FiniteElementSpace& hFESpace_)
1059 : Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
1060 hFESpace(hFESpace_)
1061 {
1062 // Assuming the same element type
1063 Mesh* mesh = lFESpace.GetMesh();
1064 dim = mesh->Dimension();
1065 if (mesh->GetNE() == 0)
1066 {
1067 return;
1068 }
1069 const FiniteElement& el = *lFESpace.GetFE(0);
1070
1071 const TensorBasisElement* ltel =
1072 dynamic_cast<const TensorBasisElement*>(&el);
1073 MFEM_VERIFY(ltel, "Low order FE space must be tensor product space");
1074
1075 const TensorBasisElement* htel =
1076 dynamic_cast<const TensorBasisElement*>(hFESpace.GetFE(0));
1077 MFEM_VERIFY(htel, "High order FE space must be tensor product space");
1078 const Array<int>& hdofmap = htel->GetDofMap();
1079
1080 const IntegrationRule& ir = hFESpace.GetFE(0)->GetNodes();
1081 IntegrationRule irLex = ir;
1082
1083 // The quadrature points, or equivalently, the dofs of the high order space
1084 // must be sorted in lexicographical order
1085 for (int i = 0; i < ir.GetNPoints(); ++i)
1086 {
1087 irLex.IntPoint(i) = ir.IntPoint(hdofmap[i]);
1088 }
1089
1090 NE = lFESpace.GetNE();
1091 const DofToQuad& maps = el.GetDofToQuad(irLex, DofToQuad::TENSOR);
1092
1093 D1D = maps.ndof;
1094 Q1D = maps.nqpt;
1095 B = maps.B;
1096 Bt = maps.Bt;
1097
1098 elem_restrict_lex_l =
1099 lFESpace.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
1100
1101 MFEM_VERIFY(elem_restrict_lex_l,
1102 "Low order ElementRestriction not available");
1103
1104 elem_restrict_lex_h =
1105 hFESpace.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
1106
1107 MFEM_VERIFY(elem_restrict_lex_h,
1108 "High order ElementRestriction not available");
1109
1110 localL.SetSize(elem_restrict_lex_l->Height(), Device::GetMemoryType());
1111 localH.SetSize(elem_restrict_lex_h->Height(), Device::GetMemoryType());
1112 localL.UseDevice(true);
1113 localH.UseDevice(true);
1114
1115 MFEM_VERIFY(dynamic_cast<const ElementRestriction*>(elem_restrict_lex_h),
1116 "High order element restriction is of unsupported type");
1117
1118 mask.SetSize(localH.Size(), Device::GetMemoryType());
1119 static_cast<const ElementRestriction*>(elem_restrict_lex_h)
1120 ->BooleanMask(mask);
1121 mask.UseDevice(true);
1122 }
1123
1124 namespace TransferKernels
1125 {
Prolongation2D(const int NE,const int D1D,const int Q1D,const Vector & localL,Vector & localH,const Array<double> & B,const Vector & mask)1126 void Prolongation2D(const int NE, const int D1D, const int Q1D,
1127 const Vector& localL, Vector& localH,
1128 const Array<double>& B, const Vector& mask)
1129 {
1130 auto x_ = Reshape(localL.Read(), D1D, D1D, NE);
1131 auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, NE);
1132 auto B_ = Reshape(B.Read(), Q1D, D1D);
1133 auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
1134
1135 localH = 0.0;
1136
1137 MFEM_FORALL(e, NE,
1138 {
1139 for (int dy = 0; dy < D1D; ++dy)
1140 {
1141 double sol_x[MAX_Q1D];
1142 for (int qy = 0; qy < Q1D; ++qy)
1143 {
1144 sol_x[qy] = 0.0;
1145 }
1146 for (int dx = 0; dx < D1D; ++dx)
1147 {
1148 const double s = x_(dx, dy, e);
1149 for (int qx = 0; qx < Q1D; ++qx)
1150 {
1151 sol_x[qx] += B_(qx, dx) * s;
1152 }
1153 }
1154 for (int qy = 0; qy < Q1D; ++qy)
1155 {
1156 const double d2q = B_(qy, dy);
1157 for (int qx = 0; qx < Q1D; ++qx)
1158 {
1159 y_(qx, qy, e) += d2q * sol_x[qx];
1160 }
1161 }
1162 }
1163 for (int qy = 0; qy < Q1D; ++qy)
1164 {
1165 for (int qx = 0; qx < Q1D; ++qx)
1166 {
1167 y_(qx, qy, e) *= m_(qx, qy, e);
1168 }
1169 }
1170 });
1171 }
1172
Prolongation3D(const int NE,const int D1D,const int Q1D,const Vector & localL,Vector & localH,const Array<double> & B,const Vector & mask)1173 void Prolongation3D(const int NE, const int D1D, const int Q1D,
1174 const Vector& localL, Vector& localH,
1175 const Array<double>& B, const Vector& mask)
1176 {
1177 auto x_ = Reshape(localL.Read(), D1D, D1D, D1D, NE);
1178 auto y_ = Reshape(localH.ReadWrite(), Q1D, Q1D, Q1D, NE);
1179 auto B_ = Reshape(B.Read(), Q1D, D1D);
1180 auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
1181
1182 localH = 0.0;
1183
1184 MFEM_FORALL(e, NE,
1185 {
1186 for (int dz = 0; dz < D1D; ++dz)
1187 {
1188 double sol_xy[MAX_Q1D][MAX_Q1D];
1189 for (int qy = 0; qy < Q1D; ++qy)
1190 {
1191 for (int qx = 0; qx < Q1D; ++qx)
1192 {
1193 sol_xy[qy][qx] = 0.0;
1194 }
1195 }
1196 for (int dy = 0; dy < D1D; ++dy)
1197 {
1198 double sol_x[MAX_Q1D];
1199 for (int qx = 0; qx < Q1D; ++qx)
1200 {
1201 sol_x[qx] = 0;
1202 }
1203 for (int dx = 0; dx < D1D; ++dx)
1204 {
1205 const double s = x_(dx, dy, dz, e);
1206 for (int qx = 0; qx < Q1D; ++qx)
1207 {
1208 sol_x[qx] += B_(qx, dx) * s;
1209 }
1210 }
1211 for (int qy = 0; qy < Q1D; ++qy)
1212 {
1213 const double wy = B_(qy, dy);
1214 for (int qx = 0; qx < Q1D; ++qx)
1215 {
1216 sol_xy[qy][qx] += wy * sol_x[qx];
1217 }
1218 }
1219 }
1220 for (int qz = 0; qz < Q1D; ++qz)
1221 {
1222 const double wz = B_(qz, dz);
1223 for (int qy = 0; qy < Q1D; ++qy)
1224 {
1225 for (int qx = 0; qx < Q1D; ++qx)
1226 {
1227 y_(qx, qy, qz, e) += wz * sol_xy[qy][qx];
1228 }
1229 }
1230 }
1231 }
1232 for (int qz = 0; qz < Q1D; ++qz)
1233 {
1234 for (int qy = 0; qy < Q1D; ++qy)
1235 {
1236 for (int qx = 0; qx < Q1D; ++qx)
1237 {
1238 y_(qx, qy, qz, e) *= m_(qx, qy, qz, e);
1239 }
1240 }
1241 }
1242 });
1243 }
1244
Restriction2D(const int NE,const int D1D,const int Q1D,const Vector & localH,Vector & localL,const Array<double> & Bt,const Vector & mask)1245 void Restriction2D(const int NE, const int D1D, const int Q1D,
1246 const Vector& localH, Vector& localL,
1247 const Array<double>& Bt, const Vector& mask)
1248 {
1249 auto x_ = Reshape(localH.Read(), Q1D, Q1D, NE);
1250 auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, NE);
1251 auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
1252 auto m_ = Reshape(mask.Read(), Q1D, Q1D, NE);
1253
1254 localL = 0.0;
1255
1256 MFEM_FORALL(e, NE,
1257 {
1258 for (int qy = 0; qy < Q1D; ++qy)
1259 {
1260 double sol_x[MAX_D1D];
1261 for (int dx = 0; dx < D1D; ++dx)
1262 {
1263 sol_x[dx] = 0.0;
1264 }
1265 for (int qx = 0; qx < Q1D; ++qx)
1266 {
1267 const double s = m_(qx, qy, e) * x_(qx, qy, e);
1268 for (int dx = 0; dx < D1D; ++dx)
1269 {
1270 sol_x[dx] += Bt_(dx, qx) * s;
1271 }
1272 }
1273 for (int dy = 0; dy < D1D; ++dy)
1274 {
1275 const double q2d = Bt_(dy, qy);
1276 for (int dx = 0; dx < D1D; ++dx)
1277 {
1278 y_(dx, dy, e) += q2d * sol_x[dx];
1279 }
1280 }
1281 }
1282 });
1283 }
Restriction3D(const int NE,const int D1D,const int Q1D,const Vector & localH,Vector & localL,const Array<double> & Bt,const Vector & mask)1284 void Restriction3D(const int NE, const int D1D, const int Q1D,
1285 const Vector& localH, Vector& localL,
1286 const Array<double>& Bt, const Vector& mask)
1287 {
1288 auto x_ = Reshape(localH.Read(), Q1D, Q1D, Q1D, NE);
1289 auto y_ = Reshape(localL.ReadWrite(), D1D, D1D, D1D, NE);
1290 auto Bt_ = Reshape(Bt.Read(), D1D, Q1D);
1291 auto m_ = Reshape(mask.Read(), Q1D, Q1D, Q1D, NE);
1292
1293 localL = 0.0;
1294
1295 MFEM_FORALL(e, NE,
1296 {
1297 for (int qz = 0; qz < Q1D; ++qz)
1298 {
1299 double sol_xy[MAX_D1D][MAX_D1D];
1300 for (int dy = 0; dy < D1D; ++dy)
1301 {
1302 for (int dx = 0; dx < D1D; ++dx)
1303 {
1304 sol_xy[dy][dx] = 0;
1305 }
1306 }
1307 for (int qy = 0; qy < Q1D; ++qy)
1308 {
1309 double sol_x[MAX_D1D];
1310 for (int dx = 0; dx < D1D; ++dx)
1311 {
1312 sol_x[dx] = 0;
1313 }
1314 for (int qx = 0; qx < Q1D; ++qx)
1315 {
1316 const double s = m_(qx, qy, qz, e) * x_(qx, qy, qz, e);
1317 for (int dx = 0; dx < D1D; ++dx)
1318 {
1319 sol_x[dx] += Bt_(dx, qx) * s;
1320 }
1321 }
1322 for (int dy = 0; dy < D1D; ++dy)
1323 {
1324 const double wy = Bt_(dy, qy);
1325 for (int dx = 0; dx < D1D; ++dx)
1326 {
1327 sol_xy[dy][dx] += wy * sol_x[dx];
1328 }
1329 }
1330 }
1331 for (int dz = 0; dz < D1D; ++dz)
1332 {
1333 const double wz = Bt_(dz, qz);
1334 for (int dy = 0; dy < D1D; ++dy)
1335 {
1336 for (int dx = 0; dx < D1D; ++dx)
1337 {
1338 y_(dx, dy, dz, e) += wz * sol_xy[dy][dx];
1339 }
1340 }
1341 }
1342 }
1343 });
1344 }
1345 } // namespace TransferKernels
1346
1347
1348 TensorProductPRefinementTransferOperator::
~TensorProductPRefinementTransferOperator()1349 ~TensorProductPRefinementTransferOperator()
1350 {
1351 }
1352
Mult(const Vector & x,Vector & y) const1353 void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
1354 Vector& y) const
1355 {
1356 if (lFESpace.GetMesh()->GetNE() == 0)
1357 {
1358 return;
1359 }
1360
1361 elem_restrict_lex_l->Mult(x, localL);
1362 if (dim == 2)
1363 {
1364 TransferKernels::Prolongation2D(NE, D1D, Q1D, localL, localH, B, mask);
1365 }
1366 else if (dim == 3)
1367 {
1368 TransferKernels::Prolongation3D(NE, D1D, Q1D, localL, localH, B, mask);
1369 }
1370 else
1371 {
1372 MFEM_ABORT("TensorProductPRefinementTransferOperator::Mult not "
1373 "implemented for dim = "
1374 << dim);
1375 }
1376 elem_restrict_lex_h->MultTranspose(localH, y);
1377 }
1378
MultTranspose(const Vector & x,Vector & y) const1379 void TensorProductPRefinementTransferOperator::MultTranspose(const Vector& x,
1380 Vector& y) const
1381 {
1382 if (lFESpace.GetMesh()->GetNE() == 0)
1383 {
1384 return;
1385 }
1386
1387 elem_restrict_lex_h->Mult(x, localH);
1388 if (dim == 2)
1389 {
1390 TransferKernels::Restriction2D(NE, D1D, Q1D, localH, localL, Bt, mask);
1391 }
1392 else if (dim == 3)
1393 {
1394 TransferKernels::Restriction3D(NE, D1D, Q1D, localH, localL, Bt, mask);
1395 }
1396 else
1397 {
1398 MFEM_ABORT("TensorProductPRefinementTransferOperator::MultTranspose not "
1399 "implemented for dim = "
1400 << dim);
1401 }
1402 elem_restrict_lex_l->MultTranspose(localL, y);
1403 }
1404
1405 #ifdef MFEM_USE_MPI
TrueTransferOperator(const ParFiniteElementSpace & lFESpace_,const ParFiniteElementSpace & hFESpace_)1406 TrueTransferOperator::TrueTransferOperator(const
1407 ParFiniteElementSpace& lFESpace_,
1408 const ParFiniteElementSpace& hFESpace_)
1409 : Operator(hFESpace_.GetTrueVSize(), lFESpace_.GetTrueVSize()),
1410 lFESpace(lFESpace_),
1411 hFESpace(hFESpace_)
1412 {
1413 localTransferOperator = new TransferOperator(lFESpace_, hFESpace_);
1414
1415 tmpL.SetSize(lFESpace_.GetVSize());
1416 tmpH.SetSize(hFESpace_.GetVSize());
1417
1418 hFESpace.GetRestrictionMatrix()->BuildTranspose();
1419 }
1420
~TrueTransferOperator()1421 TrueTransferOperator::~TrueTransferOperator()
1422 {
1423 delete localTransferOperator;
1424 }
1425
Mult(const Vector & x,Vector & y) const1426 void TrueTransferOperator::Mult(const Vector& x, Vector& y) const
1427 {
1428 lFESpace.GetProlongationMatrix()->Mult(x, tmpL);
1429 localTransferOperator->Mult(tmpL, tmpH);
1430 hFESpace.GetRestrictionMatrix()->Mult(tmpH, y);
1431 }
1432
MultTranspose(const Vector & x,Vector & y) const1433 void TrueTransferOperator::MultTranspose(const Vector& x, Vector& y) const
1434 {
1435 hFESpace.GetRestrictionMatrix()->MultTranspose(x, tmpH);
1436 localTransferOperator->MultTranspose(tmpH, tmpL);
1437 lFESpace.GetProlongationMatrix()->MultTranspose(tmpL, y);
1438 }
1439 #endif
1440
1441 } // namespace mfem
1442