22 : dom_fes(dom_fes_), ran_fes(ran_fes_),
24 fw_t_oper(), bw_t_oper(), use_ea(false), d_mt(
Device::GetHostMemoryType())
29 MFEM_VERIFY(par_dom == par_ran,
"the domain and range FE spaces must both"
30 " be either serial or parallel");
51 MFEM_VERIFY(mat != NULL,
"Operator is not a SparseMatrix");
67 const int RP_case = bool(out_cR) + 2*bool(in_cP);
84 out_cR, &oper, in_cP,
false,
false,
false));
90 MFEM_ABORT(
"Operator::Type is not supported: " <<
oper_type);
121 MFEM_ABORT(
"unknown Operator type");
129 false,
false,
false));
133 MFEM_ABORT(
"Operator::Type is not supported: " <<
oper_type);
138 return *t_oper.
Ptr();
173 for (
int i = 0; i < elem_geoms.
Size(); i++)
176 localP[elem_geoms[i]]);
184 MFEM_ABORT(
"Operator::Type is not supported: " <<
oper_type);
214 MFEM_ABORT(
"unknown type of FE space");
225 MFEM_ABORT(
"Operator::Type is not supported: " <<
oper_type);
235 :
Operator(fes_lor_.GetVSize(), fes_ho_.GetVSize()),
236 fes_ho(fes_ho_), fes_lor(fes_lor_), d_mt(d_mt_)
244 ho2lor.MakeI(nel_ho);
245 for (
int ilor = 0; ilor < nel_lor; ++ilor)
248 ho2lor.AddAColumnInRow(iho);
251 for (
int ilor = 0; ilor < nel_lor; ++ilor)
254 ho2lor.AddConnection(iho, ilor);
313 for (
int j=0; j<shape_lor.
Size(); ++j)
315 B_L(i, j) = shape_lor(j);
318 for (
int j=0; j<shape_ho.
Size(); ++j)
320 B_H(i, j) = shape_ho(j);
333 int nel_ho = mesh_ho->
GetNE();
334 int nel_lor = mesh_lor->
GetNE();
341 for (
int ig = 0; ig < geoms.
Size(); ++ig)
347 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
359 ho2lor.GetRow(iho, lor_els);
360 int nref = ho2lor.RowSize(iho);
380 MFEM_ASSERT(nel_ho*nref == nel_lor,
"we expect nel_ho*nref == nel_lor");
386 const auto d_D =
Reshape(D.
Write(), qPts, nref, nel_ho);
388 mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (
int tid)
390 const int q = tid % qPts;
391 const int iref = (tid / qPts) % nref;
392 const int iho = (tid / (qPts * nref)) % nel_ho;
394 const int lo_el_id = iref + nref*iho;
395 const real_t detJ = J(q, lo_el_id);
397 d_D(q, iref, iho) = W(q) * detJ;
405 for (
int iref = 0; iref < nref; ++iref)
407 int ilor = lor_els[iref];
419 ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, b_lo, b_ho);
431 ho2lor.GetRow(iho, lor_els);
432 int nref = ho2lor.RowSize(iho);
436 const int ndof_ho = fe_ho.
GetDof();
437 const int ndof_lor = fe_lor.
GetDof();
439 const int qPts = D.
SizeI();
446 auto v_M_LH =
Reshape(M_LH.
Write(), ndof_lor, ndof_ho, nref,
449 const int fe_ho_ndof = fe_ho.
GetDof();
450 const int fe_lor_ndof = fe_lor.
GetDof();
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);
456 mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (
int idx)
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;
462 for (
int bl = 0; bl < fe_lor_ndof; ++bl)
465 for (
int qi=0; qi<qPts; ++qi)
467 dot += d_B_L(qi, bl, iref) * d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
470 v_M_LH(bl, bh, iref, iho) = dot;
490 int nel_ho = mesh_ho->
GetNE();
491 int nel_lor = mesh_lor->
GetNE();
498 if (nel_ho == 0) {
return; }
505 for (
int ig = 0; ig < geoms.
Size(); ++ig)
515 for (
int iho = 0; iho < nel_ho; ++iho)
534 for (
int iho = 0; iho < nel_ho; ++iho)
543 int ndof_ho = fe_ho.
GetDof();
544 int ndof_lor = fe_lor.
GetDof();
553 DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
569 for (
int iref = 0; iref < nref; ++iref)
572 int ilor = lor_els[iref];
575 M_lor.
CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
579 Minv_lor.
CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
589 ElemMixedMass(geom, fe_ho, fe_lor, tr_ho, tr_lor, ip_tr, M_mixed_el);
591 M_mixed.
CopyMN(M_mixed_el, iref*ndof_lor, 0);
602 RtMlorR_inv.
Mult(RtMlor, P_iho);
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();
618 const bool build_P = fes_lor.GetTrueVSize() >= fes_ho.GetTrueVSize();
621 if (nel_ho == 0) {
return; }
628 for (
int ig = 0; ig < geoms.
Size(); ++ig)
634 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
636 offsets.SetSize(nel_ho+1);
638 for (
int iho = 0; iho < nel_ho; ++iho)
640 int nref = ho2lor.RowSize(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;
648 R.SetSize(offsets[nel_ho]);
653 P.SetSize(offsets[nel_ho]);
658 MixedMassEA(fes_ho, fes_lor, M_mixed_all,
d_mt);
674 ho2lor.GetRow(iho, lor_els);
675 nref = ho2lor.RowSize(iho);
680 ndof_lor = fe_lor.
GetDof();
685 const bool add =
false;
689 Minv_ear_lor.
SetSize(ndof_lor, ndof_lor, nel_lor,
d_mt);
696 auto v_M_mixed_all =
Reshape(M_mixed_all.
Read(), ndof_lor, ndof_ho, nref,
700 auto v_Minv_ear_lor =
Reshape(Minv_ear_lor.
Read(), ndof_lor, ndof_lor,
704 auto v_R =
Reshape(R.Write(), ndof_lor, nref, ndof_ho, nel_ho);
706 MFEM_VERIFY(nel_lor==nel_ho*nref,
"nel_lor != nel_ho*nref");
709 mfem::forall(ndof_lor * nref * ndof_ho * nel_ho, [=] MFEM_HOST_DEVICE (
int tid)
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;
717 const int lor_idx = iref + iho * nref;
721 for (
int k=0; k<ndof_lor; ++k)
723 dot += v_Minv_ear_lor(i, k, lor_idx) * v_M_mixed_all(k, j, iref, iho);
725 v_R(i, iref, j, iho) = dot;
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);
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);
744 mfem::forall(ndof_lor * nref * ndof_ho * nel_ho, [=] MFEM_HOST_DEVICE (
int tid)
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;
752 const int lor_idx = iref + e * nref;
755 for (
int t=0; t<ndof_lor; ++t)
757 dot += v_R(t, iref, iho, e) * v_M_ea_lor(t, jlo, lor_idx);
760 v_RtM_L(iho, jlo, iref, e) = dot;
771 Vector RtM_LR(ndof_ho * ndof_ho * nel_ho,
d_mt);
802 int vdim = fes_ho.GetVDim();
805 for (
int iho = 0; iho < fes_ho.GetNE(); ++iho)
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);
814 fes_ho.GetElementVDofs(iho, vdofs);
818 for (
int iref = 0; iref < nref; ++iref)
820 int ilor = ho2lor.GetRow(iho)[iref];
821 for (
int vd=0; vd<vdim; ++vd)
823 fes_lor.GetElementDofs(ilor, vdofs);
824 fes_lor.DofsToVDofs(vd, vdofs);
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();
841 R_dt.
NewMemoryAndSize(R.GetMemory(), ndof_lor*nref, ndof_ho, nel_ho,
false);
851 return EAMultTranspose(x,y);
854 int vdim = fes_ho.GetVDim();
858 for (
int iho = 0; iho < fes_ho.GetNE(); ++iho)
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);
868 for (
int iref=0; iref<nref; ++iref)
870 int ilor = ho2lor.GetRow(iho)[iref];
871 for (
int vd=0; vd<vdim; ++vd)
873 fes_lor.GetElementDofs(ilor, vdofs);
874 fes_lor.DofsToVDofs(vd, vdofs);
881 fes_ho.GetElementVDofs(iho, vdofs);
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();
897 R_dt.
NewMemoryAndSize(R.GetMemory(), ndof_lor*nref, ndof_ho, nel_ho,
false);
905 if (fes_ho.GetNE() == 0) {
return; }
909 return EAProlongate(x,y);
912 MFEM_VERIFY(P.Size() > 0,
"Prolongation not supported for these spaces.")
913 int vdim = fes_ho.GetVDim();
917 for (
int iho = 0; iho < fes_ho.GetNE(); ++iho)
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);
927 for (
int iref = 0; iref < nref; ++iref)
929 int ilor = ho2lor.GetRow(iho)[iref];
930 for (
int vd = 0; vd < vdim; ++vd)
932 fes_lor.GetElementDofs(ilor, vdofs);
933 fes_lor.DofsToVDofs(vd, vdofs);
940 fes_ho.GetElementVDofs(iho, vdofs);
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();
956 P_dt.
NewMemoryAndSize(P.GetMemory(), ndof_ho, ndof_lor * nref, nel_ho,
false);
966 return EAProlongateTranspose(x,y);
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();
975 for (
int iho = 0; iho < fes_ho.GetNE(); ++iho)
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);
984 fes_ho.GetElementVDofs(iho, vdofs);
989 for (
int iref = 0; iref < nref; ++iref)
991 int ilor = ho2lor.GetRow(iho)[iref];
992 for (
int vd=0; vd<vdim; ++vd)
994 fes_lor.GetElementDofs(ilor, vdofs);
995 fes_lor.DofsToVDofs(vd, vdofs);
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();
1013 P_dt.
NewMemoryAndSize(P.GetMemory(), ndof_ho, ndof_lor * nref, nel_ho,
false);
1038 std::unique_ptr<SparseMatrix> R_mat, M_LH_mat;
1049 R_mat.reset(
RAP(*P_lor, *R_mat, *P_ho));
1050 M_LH_mat.reset(
RAP(*P_lor, *M_LH_mat, *P_ho));
1055 M_LH_mat.reset(
mfem::Mult(*M_LH_mat, *P_ho));
1060 M_LH_mat.reset(
mfem::Mult(*P_lor, *M_LH_mat));
1069 R = std::move(R_mat);
1070 M_LH = std::move(M_LH_mat);
1081 use_ea(use_ea_), pcg(pfes_ho.GetComm())
1119 std::unique_ptr<HypreParMatrix> R_T(R_mat->
Transpose());
1126 M_LH.reset(M_LH_mat);
1140 pcg.SetPrintLevel(0);
1142 pcg.SetMaxIter(1000);
1144 pcg.SetRelTol(1e-13);
1145 pcg.SetAbsTol(1e-13);
1146 pcg.SetPreconditioner(*precon);
1147 pcg.SetOperator(*RTxM_LH);
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();
1170 for (
int ig = 0; ig < geoms.
Size(); ++ig)
1176 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
1198 Mho.
Mult(ones_ho, M_H);
1206 Mlor.
Mult(ones_lor, ML_inv_ea);
1209 LumpedMassInverse(ML_inv_ea);
1212 MixedMassEA(fes_ho, fes_lor, M_LH_ea,
d_mt);
1216 fes_lor_scalar.get(),
1221 fes_lor_scalar.get(),
1223 M_LH.reset(M_LH_local_op);
1242 int nel_ho = mesh_ho->
GetNE();
1243 int nel_lor = mesh_lor->
GetNE();
1245 int ndof_lor = pfes_lor.
GetNDofs();
1259 for (
int ig = 0; ig < geoms.
Size(); ++ig)
1265 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
1287 pMho.
Mult(ones_ho, M_H);
1295 pMlor.
Mult(ones_lor, ML_inv_ea);
1299 LumpedMassInverse(ML_inv_ea);
1302 MixedMassEA(*pfes_ho_scalar.get(), *pfes_lor_scalar.get(), M_LH_ea,
d_mt);
1306 pfes_lor_scalar.get(),
1309 const Operator *P_ho = pfes_ho_scalar->GetProlongationMatrix();
1310 const Operator *P_lor = pfes_lor_scalar->GetProlongationMatrix();
1319 RML_inv.SetSize(pfes_lor_scalar->GetTrueVSize());
1320 GetTDofs(*pfes_lor_scalar, ML_inv_ea, RML_inv);
1322 pfes_lor_scalar.get(),
1327 Vector RM_H(pfes_ho_scalar->GetTrueVSize());
1328 GetTDofsTranspose(*pfes_ho_scalar, M_H, RM_H);
1334 pfes_lor_scalar.get(),
1338 Vector RM_H(pfes_ho_scalar->GetTrueVSize());
1339 GetTDofsTranspose(*pfes_ho_scalar.get(), M_H, RM_H);
1345 RML_inv.SetSize(pfes_lor_scalar->GetTrueVSize());
1346 GetTDofsTranspose(*pfes_lor_scalar, ML_inv_ea, RML_inv);
1348 pfes_lor_scalar.get(),
1359 pfes_lor_scalar.get(),
1361 M_LH.reset(M_LH_local_op);
1380 Vector X(fes_ho.GetTrueVSize());
1381 Vector X_dim(R->Width());
1383 Vector Y_dim(R->Height());
1384 Vector Y(fes_lor.GetTrueVSize());
1388 GetTDofs(fes_ho, x, X);
1390 for (
int d = 0; d < fes_ho.GetVDim(); ++d)
1392 TDofsListByVDim(fes_ho, d, vdofs_list);
1394 R->Mult(X_dim, Y_dim);
1395 TDofsListByVDim(fes_lor, d, vdofs_list);
1399 SetFromTDofs(fes_lor, Y, y);
1406 Vector X(fes_lor.GetTrueVSize());
1407 Vector X_dim(R->Height());
1409 Vector Y_dim(R->Width());
1410 Vector Y(fes_ho.GetTrueVSize());
1414 GetTDofsTranspose(fes_lor, x, X);
1416 for (
int d = 0; d < fes_ho.GetVDim(); ++d)
1418 TDofsListByVDim(fes_lor, d, vdofs_list);
1420 R->MultTranspose(X_dim, Y_dim);
1421 TDofsListByVDim(fes_ho, d, vdofs_list);
1425 SetFromTDofsTranspose(fes_ho, Y, y);
1433 Vector X(fes_lor.GetTrueVSize());
1434 Vector X_dim(M_LH->Height());
1435 Vector Xbar(pcg.Width());
1437 Vector Y_dim(pcg.Height());
1439 Vector Y(fes_ho.GetTrueVSize());
1443 GetTDofs(fes_lor, x, X);
1445 for (
int d = 0; d < fes_ho.GetVDim(); ++d)
1447 TDofsListByVDim(fes_lor, d, vdofs_list);
1450 M_LH->MultTranspose(X_dim, Xbar);
1451 pcg.Mult(Xbar, Y_dim);
1452 TDofsListByVDim(fes_ho, d, vdofs_list);
1456 SetFromTDofs(fes_ho, Y, y);
1463 Vector X(fes_ho.GetTrueVSize());
1464 Vector X_dim(pcg.Width());
1465 Vector Xbar(pcg.Height());
1467 Vector Y_dim(M_LH->Height());
1468 Vector Y(fes_lor.GetTrueVSize());
1472 GetTDofsTranspose(fes_ho, x, X);
1474 for (
int d = 0; d < fes_ho.GetVDim(); ++d)
1476 TDofsListByVDim(fes_ho, d, vdofs_list);
1480 pcg.Mult(X_dim, Xbar);
1481 M_LH->Mult(Xbar, Y_dim);
1482 TDofsListByVDim(fes_lor, d, vdofs_list);
1486 SetFromTDofsTranspose(fes_lor, Y, y);
1492 pcg.SetRelTol(p_rtol_);
1497 pcg.SetAbsTol(p_atol_);
1501std::unique_ptr<SparseMatrix>,
1502std::unique_ptr<SparseMatrix>>
1505 std::pair<std::unique_ptr<SparseMatrix>,
1506 std::unique_ptr<SparseMatrix>> r_and_mlh;
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();
1517 return std::make_pair(
1528 for (
int ig = 0; ig < geoms.
Size(); ++ig)
1534 BuildHo2Lor(nel_ho, nel_lor, cf_tr);
1543 for (
int iho = 0; iho < nel_ho; ++iho)
1546 ho2lor.GetRow(iho, lor_els);
1547 int nref = ho2lor.RowSize(iho);
1551 int nedof_lor = fe_lor.
GetDof();
1556 Vector shape_lor(nedof_lor);
1559 for (
int iref = 0; iref < nref; ++iref)
1561 int ilor = lor_els[iref];
1572 ML_el += (shape_lor *= (el_tr->
Weight() * ip_lor.
weight));
1574 fes_lor.GetElementDofs(ilor, dofs_lor);
1579 LumpedMassInverse(ML_inv);
1582 r_and_mlh.first = AllocR();
1587 for (
int icol = 0; icol < r_and_mlh.first->Height() + 1; ++icol)
1589 I[icol] = r_and_mlh.first->GetI()[icol];
1591 Memory<int> J(r_and_mlh.first->NumNonZeroElems());
1592 for (
int jcol = 0; jcol < r_and_mlh.first->NumNonZeroElems(); ++jcol)
1594 J[jcol] = r_and_mlh.first->GetJ()[jcol];
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));
1604 offsets.SetSize(nel_ho+1);
1606 for (
int iho = 0; iho < nel_ho; ++iho)
1609 ho2lor.GetRow(iho, lor_els);
1610 int nref = ho2lor.RowSize(iho);
1615 offsets[iho+1] = offsets[iho] + fe_ho.
GetDof()*fe_lor.
GetDof()*nref;
1622 int nedof_ho = fe_ho.
GetDof();
1623 int nedof_lor = fe_lor.
GetDof();
1627 for (
int iref = 0; iref < nref; ++iref)
1629 int ilor = lor_els[iref];
1636 ElemMixedMass(geom, fe_ho, fe_lor, tr_ho, tr_lor, ip_tr, M_LH_el);
1639 fes_lor.GetElementDofs(ilor, dofs_lor);
1641 for (
int i = 0; i < nedof_lor; ++i)
1643 M_LH_el.
GetRow(i, R_row);
1644 R_el.
SetRow(i, R_row.
Set(ML_inv[dofs_lor[i]], R_row));
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);
1724 R_mat->
BooleanMult(x_vdofs_marker, X_vdofs_marker);
1741 auto * fes = pfes_lor_scalar ==
nullptr ? fes_lor_scalar.get() :
1742 pfes_lor_scalar.get();
1744 auto * fes = fes_lor_scalar.get();
1746 MFEM_ASSERT(fes !=
nullptr,
"[p]fes_lor_scalar is nullptr");
1748 Vector ML_inv_true(fes->GetTrueVSize());
1749 const Operator *P = fes->GetProlongationMatrix();
1751 else { ML_inv_true = ML_inv; }
1755 if (P) { P->
Mult(ML_inv_true, ML_inv); }
1756 else { ML_inv = ML_inv_true; }
1760std::unique_ptr<SparseMatrix>
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();
1769 Transpose(elem_dof_lor, dof_elem_lor, ndof_lor);
1771 Mesh* mesh_lor = fes_lor.GetMesh();
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();
1784 dof_used_ho.
SetSize(ndof_ho, -1);
1787 for (
int ilor = 0; ilor < ndof_lor; ++ilor)
1789 for (
int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
1791 int el_lor = dof_elem_lorJ[jlor];
1793 for (
int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
1795 int dof_ho = elem_dof_hoJ[jho];
1796 if (dof_used_ho[dof_ho] != ilor)
1798 dof_used_ho[dof_ho] = ilor;
1806 Table dof_lor_dof_ho;
1807 dof_lor_dof_ho.
SetDims(ndof_lor, sizeJ);
1809 for (
int i = 0; i < ndof_ho; ++i)
1811 dof_used_ho[i] = -1;
1815 int* dof_dofI = dof_lor_dof_ho.
GetI();
1816 int* dof_dofJ = dof_lor_dof_ho.
GetJ();
1818 for (
int ilor = 0; ilor < ndof_lor; ++ilor)
1820 dof_dofI[ilor] = sizeJ;
1821 for (
int jlor = dof_elem_lorI[ilor]; jlor < dof_elem_lorI[ilor + 1]; ++jlor)
1823 int el_lor = dof_elem_lorJ[jlor];
1825 for (
int jho = elem_dof_hoI[iho]; jho < elem_dof_hoI[iho + 1]; ++jho)
1827 int dof_ho = elem_dof_hoJ[jho];
1828 if (dof_used_ho[dof_ho] != ilor)
1830 dof_used_ho[dof_ho] = ilor;
1831 dof_dofJ[sizeJ] = dof_ho;
1841 std::unique_ptr<SparseMatrix> R_local(
new SparseMatrix(
1842 dof_dofI, dof_dofJ, data, ndof_lor,
1843 ndof_ho,
true,
true,
true));
1856 fes_ho(fes_ho_), fes_lor(fes_lor_), ho2lor(ho2lor_),
1863 const Operator* elem_restrict_ho = fes_ho->GetElementRestriction(
1865 const Operator* elem_restrict_lor = fes_lor->GetElementRestriction(
1868 const int vdim = fes_ho->GetVDim();
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();
1877 elem_restrict_ho->
Mult(x, tempx);
1879 Vector tempy(ndof_lor*nref*vdim*nel_ho);
1881 auto v_M_mixed_ea =
Reshape(M_LH_ea->Read(), ndof_lor, ndof_ho, nref,
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);
1887 mfem::forall(ndof_lor * nref * vdim * nel_ho, [=] MFEM_HOST_DEVICE (
int tid)
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;
1895 for (
int k=0; k<ndof_ho; ++k)
1897 dot += v_M_mixed_ea(j, k, i, iho) * v_tempx(k, v, iho);
1900 v_tempy(j, i, v, iho) = dot;
1909 const Operator* elem_restrict_ho = fes_ho->GetElementRestriction(
1911 const Operator* elem_restrict_lor = fes_lor->GetElementRestriction(
1914 const int vdim = fes_ho->GetVDim();
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();
1923 elem_restrict_lor->
Mult(x, tempx);
1925 Vector tempy(ndof_ho*vdim*nel_ho);
1927 auto v_M_mixed_ea =
Reshape(M_LH_ea->Read(), ndof_lor, ndof_ho, nref,
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);
1932 mfem::forall(ndof_ho * vdim * nel_ho, [=] MFEM_HOST_DEVICE (
int tid)
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;
1939 for (
int i=0; i<nref; ++i)
1941 for (
int j=0; j<ndof_lor; ++j)
1943 dot += v_M_mixed_ea(j, k, i, iho) * v_tempx(j, i, v, iho);
1945 v_tempy(k, v, iho) = dot;
1956 Operator(ML_inv_.Size(), ML_inv_.Size()),
1957 fes_ho(fes_ho_), fes_lor(fes_lor_),
1964 MFEM_ASSERT(ML_inv->Size() == x.
Size(),
"sizes not the same");
1965 auto v_ML_inv =
Reshape(ML_inv->Read(), ML_inv->Size());
1969 mfem::forall(ML_inv->Size(), [=] MFEM_HOST_DEVICE(
int i)
1970 { v_y(i) = v_ML_inv(i) * v_x(i); });
1987 if (!
F) { BuildF(); }
1995 if (!
F) { BuildF(); }
2001void L2ProjectionGridTransfer::BuildF()
2018 F =
new L2ProjectionH1Space(dom_pfes, ran_pfes,
2038 :
Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize())
2041 if (lFESpace_.
FEColl() == hFESpace_.
FEColl() && !isvar_order)
2048 else if (lFESpace_.
GetVDim() == 1
2081 :
Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
2099 int vdim = lFESpace.
GetVDim();
2103 for (
int i = 0; i < mesh->
GetNE(); i++)
2109 if (geom != cached_geom || isvar_order)
2111 h_fe = hFESpace.
GetFE(i);
2112 l_fe = lFESpace.
GetFE(i);
2119 for (
int vd = 0; vd < vdim; vd++)
2121 l_dofs.
Copy(l_vdofs);
2123 h_dofs.
Copy(h_vdofs);
2130 loc_prol.
Mult(subX, subY);
2158 int vdim = lFESpace.
GetVDim();
2160 for (
int i = 0; i < mesh->
GetNE(); i++)
2166 if (geom != cached_geom || isvar_order)
2168 h_fe = hFESpace.
GetFE(i);
2169 l_fe = lFESpace.
GetFE(i);
2177 for (
int vd = 0; vd < vdim; vd++)
2179 l_dofs.
Copy(l_vdofs);
2181 h_dofs.
Copy(h_vdofs);
2189 for (
int p = 0;
p < h_dofs.
Size(); ++
p)
2191 if (processed[lFESpace.
DecodeDof(h_dofs[
p])])
2197 loc_prol.
Mult(subX, subY);
2205 for (
int p = 0;
p < h_dofs.
Size(); ++
p)
2207 processed[lFESpace.
DecodeDof(h_dofs[
p])] = 1;
2217 :
Operator(hFESpace_.GetVSize(), lFESpace_.GetVSize()), lFESpace(lFESpace_),
2227 MFEM_VERIFY(ltel,
"Low order FE space must be tensor product space");
2231 MFEM_VERIFY(htel,
"High order FE space must be tensor product space");
2241 int j = hdofmap[i] >=0 ? hdofmap[i] : -1 - hdofmap[i];
2245 NE = lFESpace.
GetNE();
2253 elem_restrict_lex_l =
2256 MFEM_VERIFY(elem_restrict_lex_l,
2257 "Low order ElementRestriction not available");
2259 elem_restrict_lex_h =
2262 MFEM_VERIFY(elem_restrict_lex_h,
2263 "High order ElementRestriction not available");
2271 "High order element restriction is of unsupported type");
2275 ->BooleanMask(mask);
2279namespace TransferKernels
2292 for (
int qy = 0; qy < Q1D; ++qy)
2294 for (
int qx = 0; qx < Q1D; ++qx)
2296 y_(qx, qy, e) = 0.0;
2300 for (
int dy = 0; dy < D1D; ++dy)
2302 real_t sol_x[DofQuadLimits::MAX_Q1D];
2303 for (
int qy = 0; qy < Q1D; ++qy)
2307 for (
int dx = 0; dx < D1D; ++dx)
2309 const real_t s = x_(dx, dy, e);
2310 for (
int qx = 0; qx < Q1D; ++qx)
2312 sol_x[qx] += B_(qx, dx) * s;
2315 for (
int qy = 0; qy < Q1D; ++qy)
2317 const real_t d2q = B_(qy, dy);
2318 for (
int qx = 0; qx < Q1D; ++qx)
2320 y_(qx, qy, e) += d2q * sol_x[qx];
2324 for (
int qy = 0; qy < Q1D; ++qy)
2326 for (
int qx = 0; qx < Q1D; ++qx)
2328 y_(qx, qy, e) *= m_(qx, qy, e);
2338 auto x_ =
Reshape(localL.
Read(), D1D, D1D, D1D, NE);
2341 auto m_ =
Reshape(mask.
Read(), Q1D, Q1D, Q1D, NE);
2345 for (
int qz = 0; qz < Q1D; ++qz)
2347 for (
int qy = 0; qy < Q1D; ++qy)
2349 for (
int qx = 0; qx < Q1D; ++qx)
2351 y_(qx, qy, qz, e) = 0.0;
2356 for (
int dz = 0; dz < D1D; ++dz)
2358 real_t sol_xy[DofQuadLimits::MAX_Q1D][DofQuadLimits::MAX_Q1D];
2359 for (
int qy = 0; qy < Q1D; ++qy)
2361 for (
int qx = 0; qx < Q1D; ++qx)
2363 sol_xy[qy][qx] = 0.0;
2366 for (
int dy = 0; dy < D1D; ++dy)
2368 real_t sol_x[DofQuadLimits::MAX_Q1D];
2369 for (
int qx = 0; qx < Q1D; ++qx)
2373 for (
int dx = 0; dx < D1D; ++dx)
2375 const real_t s = x_(dx, dy, dz, e);
2376 for (
int qx = 0; qx < Q1D; ++qx)
2378 sol_x[qx] += B_(qx, dx) * s;
2381 for (
int qy = 0; qy < Q1D; ++qy)
2383 const real_t wy = B_(qy, dy);
2384 for (
int qx = 0; qx < Q1D; ++qx)
2386 sol_xy[qy][qx] += wy * sol_x[qx];
2390 for (
int qz = 0; qz < Q1D; ++qz)
2392 const real_t wz = B_(qz, dz);
2393 for (
int qy = 0; qy < Q1D; ++qy)
2395 for (
int qx = 0; qx < Q1D; ++qx)
2397 y_(qx, qy, qz, e) += wz * sol_xy[qy][qx];
2402 for (
int qz = 0; qz < Q1D; ++qz)
2404 for (
int qy = 0; qy < Q1D; ++qy)
2406 for (
int qx = 0; qx < Q1D; ++qx)
2408 y_(qx, qy, qz, e) *= m_(qx, qy, qz, e);
2426 for (
int dy = 0; dy < D1D; ++dy)
2428 for (
int dx = 0; dx < D1D; ++dx)
2430 y_(dx, dy, e) = 0.0;
2434 for (
int qy = 0; qy < Q1D; ++qy)
2436 real_t sol_x[DofQuadLimits::MAX_D1D];
2437 for (
int dx = 0; dx < D1D; ++dx)
2441 for (
int qx = 0; qx < Q1D; ++qx)
2443 const real_t s = m_(qx, qy, e) * x_(qx, qy, e);
2444 for (
int dx = 0; dx < D1D; ++dx)
2446 sol_x[dx] += Bt_(dx, qx) * s;
2449 for (
int dy = 0; dy < D1D; ++dy)
2451 const real_t q2d = Bt_(dy, qy);
2452 for (
int dx = 0; dx < D1D; ++dx)
2454 y_(dx, dy, e) += q2d * sol_x[dx];
2464 auto x_ =
Reshape(localH.
Read(), Q1D, Q1D, Q1D, NE);
2467 auto m_ =
Reshape(mask.
Read(), Q1D, Q1D, Q1D, NE);
2471 for (
int dz = 0; dz < D1D; ++dz)
2473 for (
int dy = 0; dy < D1D; ++dy)
2475 for (
int dx = 0; dx < D1D; ++dx)
2477 y_(dx, dy, dz, e) = 0.0;
2482 for (
int qz = 0; qz < Q1D; ++qz)
2484 real_t sol_xy[DofQuadLimits::MAX_D1D][DofQuadLimits::MAX_D1D];
2485 for (
int dy = 0; dy < D1D; ++dy)
2487 for (
int dx = 0; dx < D1D; ++dx)
2492 for (
int qy = 0; qy < Q1D; ++qy)
2494 real_t sol_x[DofQuadLimits::MAX_D1D];
2495 for (
int dx = 0; dx < D1D; ++dx)
2499 for (
int qx = 0; qx < Q1D; ++qx)
2501 const real_t s = m_(qx, qy, qz, e) * x_(qx, qy, qz, e);
2502 for (
int dx = 0; dx < D1D; ++dx)
2504 sol_x[dx] += Bt_(dx, qx) * s;
2507 for (
int dy = 0; dy < D1D; ++dy)
2509 const real_t wy = Bt_(dy, qy);
2510 for (
int dx = 0; dx < D1D; ++dx)
2512 sol_xy[dy][dx] += wy * sol_x[dx];
2516 for (
int dz = 0; dz < D1D; ++dz)
2518 const real_t wz = Bt_(dz, qz);
2519 for (
int dy = 0; dy < D1D; ++dy)
2521 for (
int dx = 0; dx < D1D; ++dx)
2523 y_(dx, dy, dz, e) += wz * sol_xy[dy][dx];
2540 elem_restrict_lex_l->
Mult(x, localL);
2551 MFEM_ABORT(
"TensorProductPRefinementTransferOperator::Mult not "
2552 "implemented for dim = "
2566 elem_restrict_lex_h->
Mult(x, localH);
2577 MFEM_ABORT(
"TensorProductPRefinementTransferOperator::MultTranspose not "
2578 "implemented for dim = "
2587 :
Operator(hFESpace_.GetTrueVSize(), lFESpace_.GetTrueVSize()),
2588 lFESpace(lFESpace_),
2600 if (P) { MFEM_VERIFY(R,
"Both P and R have to be not NULL") }
2616 delete localTransferOperator;
2624 localTransferOperator->
Mult(tmpL, tmpH);
2629 localTransferOperator->
Mult(x, tmpH);
2634 localTransferOperator->
Mult(x, y);
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
void Copy(Array ©) const
Create a copy of the internal array to the provided copy.
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
static void Mult(const DenseTensor &A, const Vector &x, Vector &y)
Computes (e.g. by calling AddMult(A,x,y,1,0,Op::N)).
static void MultTranspose(const DenseTensor &A, const Vector &x, Vector &y)
Computes (e.g. by calling AddMult(A,x,y,1,0,Op::T)).
static void Invert(DenseTensor &A)
Replaces the block diagonal matrix with its inverse .
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
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.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void Transpose()
(*this) = (*this)^t
void SetRow(int r, const real_t *row)
real_t * GetData() const
Returns the matrix data array.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
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).
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
static MemoryType GetMemoryType()
(DEPRECATED) Equivalent to GetDeviceMemoryType().
Structure representing the matrices/tensors needed to evaluate (in reference space) the values,...
@ TENSOR
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Array< real_t > B
Basis functions evaluated at quadrature points.
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number.
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Array< real_t > Bt
Transpose of B.
Operator that converts FiniteElementSpace L-vectors to E-vectors.
virtual int GetContType() const =0
@ DISCONTINUOUS
Field is discontinuous across element interfaces.
@ CONTINUOUS
Field is continuous across element interfaces.
Derefinement operator, used by the friend class InterpolationGridTransfer.
GridFunction interpolation operator applicable after mesh refinement.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
void GetVDofs(int vd, Array< int > &dofs, int ndofs=-1) const
Returns the indices of all of the VDofs for the specified dimension 'vd'.
const SparseMatrix * GetConformingRestriction() const
The returned SparseMatrix is owned by the FiniteElementSpace.
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element,...
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
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...
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
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...
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
void GetLocalRefinementMatrices(Geometry::Type geom, DenseTensor &localP) const
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,...
virtual const Operator * GetProlongationMatrix() const
virtual const Operator * GetRestrictionOperator() const
An abstract operator that performs the same action as GetRestrictionMatrix.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
int GetNE() const
Returns number of elements in the mesh.
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
SparseMatrix * RefinementMatrix_main(const int coarse_ndofs, const Table &coarse_elem_dof, const Table *coarse_elem_fos, const DenseTensor localP[]) const
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
const SparseMatrix * GetConformingProlongation() const
const FiniteElementCollection * FEColl() const
Mesh * GetMesh() const
Returns the mesh.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
const Table * GetElementToFaceOrientationTable() const
int GetVDim() const
Returns the vector dimension of the finite element space.
virtual const SparseMatrix * GetHpRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
const FiniteElement * GetTypicalFE() const
Return GetFE(0) if the local mesh is not empty; otherwise return a typical FE based on the Geometry t...
static int DecodeDof(int dof)
Helper to return the DOF associated with a sign encoded DOF.
Abstract class for all finite elements.
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...
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
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...
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
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.
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 ...
Structure for storing mesh geometric factors: coordinates, Jacobians, and determinants of the Jacobia...
Vector detJ
Determinants of the Jacobians at all quadrature points.
FiniteElementSpace & dom_fes
Domain FE space.
FiniteElementSpace & ran_fes
Range FE space.
Operator::Type oper_type
Desired Operator::Type for the construction of all operators defined by the underlying transfer algor...
const Operator & MakeTrueOperator(FiniteElementSpace &fes_in, FiniteElementSpace &fes_out, const Operator &oper, OperatorHandle &t_oper)
GridTransfer(FiniteElementSpace &dom_fes_, FiniteElementSpace &ran_fes_)
The BoomerAMG solver in hypre.
void SetPrintLevel(int print_level)
Wrapper for hypre's ParCSR matrix class.
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...
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Class for integration point with weight.
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
const Array< real_t > & GetWeights() const
Return the quadrature weights in a contiguous array.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
virtual ~InterpolationGridTransfer()
BilinearFormIntegrator * mass_integ
Ownership depends on own_mass_integ.
OperatorHandle F
Forward, coarse-to-fine, operator.
bool own_mass_integ
Ownership flag for mass_integ.
const Operator & BackwardOperator() override
Return an Operator that transfers GridFunctions from the range FE space back to GridFunctions in the ...
OperatorHandle B
Backward, fine-to-coarse, operator.
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,...
const Operator & ForwardOperator() override
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
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 EAL2ProjectionH1Space()
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
std::unique_ptr< ParFiniteElementSpace > pfes_ho_scalar
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.
std::unique_ptr< Solver > precon
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.
std::unique_ptr< Operator > RTxM_LH
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(),...
std::unique_ptr< Operator > M_LH
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
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 ...
std::unique_ptr< Operator > R
void SetRelTol(real_t p_rtol_) override
Sets relative tolerance in preconditioned conjugate gradient solver.
std::unique_ptr< ParFiniteElementSpace > pfes_lor_scalar
void ProlongateTranspose(const Vector &x, Vector &y) const override
L2ProjectionL2Space(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, const bool use_ea_, MemoryType d_mt_=Device::GetHostMemoryType())
void EAMultTranspose(const Vector &x, Vector &y) const
void EAProlongate(const Vector &x, Vector &y) const
void MultTranspose(const Vector &x, Vector &y) const override
void EAMult(const Vector &x, Vector &y) const
Perform mult on the device (same as above)
void EAL2ProjectionL2Space()
void EAProlongateTranspose(const Vector &x, Vector &y) const
void Prolongate(const Vector &x, Vector &y) const override
void Mult(const Vector &x, Vector &y) const override
void BuildHo2Lor(int nel_ho, int nel_lor, const CoarseFineTransformations &cf_tr)
const FiniteElementSpace & fes_ho
void MixedMassEA(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, Vector &M_LH, MemoryType d_mt_=Device::GetHostMemoryType())
const FiniteElementSpace & fes_lor
L2Projection(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, MemoryType d_mt_=Device::GetHostMemoryType())
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
L2Projection * F
Forward, coarse-to-fine, operator.
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.
const Operator & ForwardOperator() override
Return an Operator that transfers GridFunctions from the domain FE space to GridFunctions in the rang...
virtual ~L2ProjectionGridTransfer()
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>.
void GetGeometries(int dim, Array< Geometry::Type > &el_geoms) const
Return all element geometries of the given dimension present in the mesh.
const CoarseFineTransformations & GetRefinementTransforms() const
int GetNE() const
Returns number of elements.
int Dimension() const
Dimension of the reference space used within the elements.
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.
Geometry::Type GetElementBaseGeometry(int i) const
Pointer to an Operator of a specified type.
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Operator * Ptr() const
Access the underlying Operator pointer.
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Jacobi smoothing for a given bilinear form (no matrix necessary).
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
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().
@ ANY_TYPE
ID for the base class Operator, i.e. any type.
@ MFEM_SPARSEMAT
ID for class SparseMatrix.
@ Hypre_ParCSR
ID for class HypreParMatrix.
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 ...
Matrix-free transfer operator between finite element spaces on the same mesh.
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...
Abstract parallel finite element space.
HYPRE_BigInt * GetTrueDofOffsets() const
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
ParMesh * GetParMesh() const
General product operator: x -> (A*B)(x) = A(B(x)).
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.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
void SortRows()
Sort the column (TYPE II) indices in each row.
void SetDims(int rows, int nnz)
const Array< int > & GetDofMap() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Matrix-free transfer operator between finite element spaces on the same mesh exploiting the tensor pr...
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.
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().
General triple product operator x -> A*B*C*x, with ownership of the factors.
~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.
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Memory< real_t > & GetMemory()
Return a reference to the Memory object used by the Vector.
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...
Vector & Set(const real_t a, const Vector &x)
(*this) = a * x
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.
int Size() const
Returns the size of the vector.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
void SetSize(int s)
Resize the vector to size s.
void Reciprocal()
(*this)(i) = 1.0 / (*this)(i)
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
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)
void add(const Vector &v1, const Vector &v2, Vector &v)
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
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)
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
MemoryType
Memory types supported by MFEM.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
@ 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)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
real_t p(const Vector &x, real_t t)