【How to implement the RT projection locally】

[cglwdm]

I'm trying to implement the following RT interpolation locally, but there is not any examples that I can find to help to do this.

I also tried to implement it globally, but I can't find the suitable integrations， AddInteriorFaceIntegrator and AddInteriorFaceIntegrator are only for DG methds.

[tzanio]

It is global, but you can use IdentityInterpolator ector H1 into the RT space.

This is done for example in the construction of AMS which needs the similar Nedelec interpolant $\Pi$ from vector H1 into the NS space as a matrix, see https://github.com/mfem/mfem/blob/master/linalg/hypre.cpp#L5614

[cglwdm]

It is global, but you can use IdentityInterpolator ector H1 into the RT space.

This is done for example in the construction of AMS which needs the similar Nedelec interpolant Π from vector H1 into the NS space as a matrix, see https://github.com/mfem/mfem/blob/master/linalg/hypre.cpp#L5614

Thank you for your help! So IdentityInterpolator and ParDiscreteLinearOperator can define a Matrix from H1 space to RT space, is it correct?
The above Pi can be calculate locally, since (18) ensures that pi v\in H(div) spaces when v\in H1. I think the above pi can also be obtained in the help of MFEM locally. Since I'm investigating the postprocessing methods in the HDG branch of MFEM, where a postprocessing method is proposed and implemented element by element. We can solve (18) and (19) element by element for any given v\in [H^1(\Omega)]^d.

[cglwdm]

The below is the RTk projection written by my own code, I will try it in MFEM to implement the similar code.

int RTkprojection(
    MESH &Mesh,
    int degreeVecSigmah,
    Matrix<dataType, 2, 1>(*VecSigma)(const MESH &Mesh, int p_k, int id, int i, int i_face, int i_edge, dataType x, dataType y),
    dataType(*curl_VecSigma)(const MESH &Mesh, int p_k, int id, int i, int i_face, int i_edge, dataType x, dataType y),
    Matrix<dataType, Dynamic, 1> &coeffVecSigmah)
{
    int freedomVecSigmah = FREEDOM_RTK[degreeVecSigmah];
    coeffVecSigmah.resize(Mesh.getNumElement() * freedomVecSigmah);
    printf("Matrix ");
    fflush(stdout);

#pragma omp parallel for 
    for (int id = 0; id < Mesh.num_Element; id++) {
        int m_k = degreeVecSigmah;
        Matrix<dataType, Dynamic, 1> Fl(freedomVecSigmah);
        Matrix<dataType, Dynamic, Dynamic> Al(freedomVecSigmah, freedomVecSigmah);
        Fl.setZero();
        Al.setZero();

        for (int i = 0; i < 3 * (m_k + 1); i++) {
            int face_k = i / (m_k + 1);
            int face_i = i;
            for (int j = 0; j < FREEDOM_RTK[m_k]; j++) {

                Al(i, j) = qf5pE(Mesh, local_RTk, normalVecPkEdge, 1.0, m_k, m_k, id, j, face_i, face_k, face_k, 0, 0, face_k);
            }
            Fl(i) +=
                qf5pE(Mesh, VecSigma, normalVecPkEdge, 1.0, m_k, m_k, id, face_i, face_i, face_k, face_k, 0, 0, face_k);
        }
        if (m_k >= 1) {
            for (int i = 0; i < FREEDOM_VECTORPK[m_k - 1]; i++) {
                for (int j = 0; j < FREEDOM_RTK[m_k]; j++) {
                    Al(i + 3 * (m_k + 1), j) = qf25pT(Mesh, vectorPk, local_RTk, 1.0, m_k - 1, m_k, id, i, j);
                }
                Fl(i + 3 * (m_k + 1)) = qf25pT(Mesh, vectorPk, VecSigma, 1.0, m_k - 1, m_k, id, i, i);
            }
        }
        coeffVecSigmah.segment(id * freedomVecSigmah, freedomVecSigmah) = Al.inverse() * Fl;
    }

    return Mesh.getNumElement() * freedomVecSigmah;

}

[cglwdm]

I have implemented the above RT projection, may I will pull a request to add this new RT projection.

[tzanio]

Yes, please feel free to open a PR