Skip to content

Commit

Permalink
Merge branch 'master' into mesh-reflect
Browse files Browse the repository at this point in the history
  • Loading branch information
@dylan-copeland
dylan-copeland authored Feb 15, 2023
2 parents 4210208 + 61b3944 commit 92e11e4
Show file tree
Hide file tree
Showing 3 changed files with 67 additions and 35 deletions.
34 changes: 27 additions & 7 deletions fem/transfer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -291,6 +291,10 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space(
int nel_ho = mesh_ho->GetNE();
int nel_lor = mesh_lor->GetNE();

// The prolongation operation is only well-defined when the LOR space has at
// least as many DOFs as the high-order space.
const bool build_P = fes_lor.GetTrueVSize() >= fes_ho.GetTrueVSize();

// If the local mesh is empty, skip all computations
if (nel_ho == 0) { return; }

Expand Down Expand Up @@ -319,8 +323,11 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space(
// R will contain the restriction (L^2 projection operator) defined on each
// coarse HO element (and corresponding patch of LOR elements)
R.SetSize(offsets[nel_ho]);
// P will contain the corresponding prolongation operator
P.SetSize(offsets[nel_ho]);
if (build_P)
{
// P will contain the corresponding prolongation operator
P.SetSize(offsets[nel_ho]);
}

IntegrationPointTransformation ip_tr;
IsoparametricTransformation &emb_tr = ip_tr.Transf;
Expand All @@ -341,7 +348,6 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space(
const DenseTensor &pmats = cf_tr.point_matrices[geom];

DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);

DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
Expand Down Expand Up @@ -385,10 +391,15 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space(
}
mfem::Mult(Minv_lor, M_mixed, R_iho);

mfem::MultAtB(R_iho, M_lor, RtMlor);
mfem::Mult(RtMlor, R_iho, RtMlorR);
RtMlorR_inv.Factor();
RtMlorR_inv.Mult(RtMlor, P_iho);
if (build_P)
{
DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);

mfem::MultAtB(R_iho, M_lor, RtMlor);
mfem::Mult(RtMlor, R_iho, RtMlorR);
RtMlorR_inv.Factor();
RtMlorR_inv.Mult(RtMlor, P_iho);
}
}
}

Expand Down Expand Up @@ -462,6 +473,8 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::MultTranspose(
void L2ProjectionGridTransfer::L2ProjectionL2Space::Prolongate(
const Vector &x, Vector &y) const
{
if (fes_ho.GetNE() == 0) { return; }
MFEM_VERIFY(P.Size() > 0, "Prolongation not supported for these spaces.")
int vdim = fes_ho.GetVDim();
Array<int> vdofs;
DenseMatrix xel_mat,yel_mat;
Expand Down Expand Up @@ -497,6 +510,8 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::Prolongate(
void L2ProjectionGridTransfer::L2ProjectionL2Space::ProlongateTranspose(
const Vector &x, Vector &y) const
{
if (fes_ho.GetNE() == 0) { return; }
MFEM_VERIFY(P.Size() > 0, "Prolongation not supported for these spaces.")
int vdim = fes_ho.GetVDim();
Array<int> vdofs;
DenseMatrix xel_mat,yel_mat;
Expand Down Expand Up @@ -898,6 +913,11 @@ void L2ProjectionGridTransfer::BuildF()
}
}

bool L2ProjectionGridTransfer::SupportsBackwardsOperator() const
{
return ran_fes.GetTrueVSize() >= dom_fes.GetTrueVSize();
}


TransferOperator::TransferOperator(const FiniteElementSpace& lFESpace_,
const FiniteElementSpace& hFESpace_)
Expand Down
4 changes: 4 additions & 0 deletions fem/transfer.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,8 @@ class GridTransfer
{
return MakeTrueOperator(ran_fes, dom_fes, BackwardOperator(), bw_t_oper);
}

virtual bool SupportsBackwardsOperator() const { return true; }
};


Expand Down Expand Up @@ -346,6 +348,8 @@ class L2ProjectionGridTransfer : public GridTransfer
virtual const Operator &ForwardOperator();

virtual const Operator &BackwardOperator();

virtual bool SupportsBackwardsOperator() const;
private:
void BuildF();
};
Expand Down
64 changes: 36 additions & 28 deletions miniapps/tools/lor-transfer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,32 +170,36 @@ int main(int argc, char *argv[])
gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
}
const Operator &R = gt->ForwardOperator();
const Operator &P = gt->BackwardOperator();

// HO->LOR restriction
direction = "HO -> LOR @ LOR";
R.Mult(rho, rho_lor);
compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO) ");
if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy); Wx += offx; }

// LOR->HO prolongation
direction = "HO -> LOR @ HO";
GridFunction rho_prev = rho;
P.Mult(rho_lor, rho);
compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy); Wx = 0; Wy += offy; }

rho_prev -= rho;
cout.precision(12);
cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl;
if (gt->SupportsBackwardsOperator())
{
const Operator &P = gt->BackwardOperator();
// LOR->HO prolongation
direction = "HO -> LOR @ HO";
GridFunction rho_prev = rho;
P.Mult(rho_lor, rho);
compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy); Wx = 0; Wy += offy; }

rho_prev -= rho;
cout.precision(12);
cout << "|HO - P(R(HO))|_∞ = " << rho_prev.Normlinf() << endl;
}

// HO* to LOR* dual fields
GridFunction ones(&fespace), ones_lor(&fespace_lor);
ones = 1.0;
ones_lor = 1.0;
LinearForm M_rho(&fespace), M_rho_lor(&fespace_lor);
if (!use_pointwise_transfer)
if (!use_pointwise_transfer && gt->SupportsBackwardsOperator())
{
const Operator &P = gt->BackwardOperator();
M_ho.Mult(rho, M_rho);
P.MultTranspose(M_rho, M_rho_lor);
cout << "HO -> LOR dual field: " << fabs(M_rho(ones)-M_rho_lor(ones_lor))
Expand All @@ -209,22 +213,26 @@ int main(int argc, char *argv[])
double lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR ");
if (vis) { visualize(LOR_dc, "LOR", Wx, Wy); Wx += offx; }

// Prolongate to HO space
direction = "LOR -> HO @ HO";
P.Mult(rho_lor, rho);
compute_mass(&fespace, lor_mass, HO_dc, "P(LOR) ");
if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy); Wx += offx; }

// Restrict back to LOR space. This won't give the original function because
// the rho_lor doesn't necessarily live in the range of R.
direction = "LOR -> HO @ LOR";
R.Mult(rho, rho_lor);
compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy); }

rho_lor_prev -= rho_lor;
cout.precision(12);
cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl;
if (gt->SupportsBackwardsOperator())
{
const Operator &P = gt->BackwardOperator();
// Prolongate to HO space
direction = "LOR -> HO @ HO";
P.Mult(rho_lor, rho);
compute_mass(&fespace, lor_mass, HO_dc, "P(LOR) ");
if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy); Wx += offx; }

// Restrict back to LOR space. This won't give the original function because
// the rho_lor doesn't necessarily live in the range of R.
direction = "LOR -> HO @ LOR";
R.Mult(rho, rho_lor);
compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy); }

rho_lor_prev -= rho_lor;
cout.precision(12);
cout << "|LOR - R(P(LOR))|_∞ = " << rho_lor_prev.Normlinf() << endl;
}

// LOR* to HO* dual fields
if (!use_pointwise_transfer)
Expand Down

0 comments on commit 92e11e4

Please sign in to comment.