def METHOD_NAME(eval_dofs, eval_coors, n_el, ps, min_level=0, max_level=2,
                  eps=1e-4):
    """
    Create mesh with linear elements that approximates DOFs returned by
    `eval_dofs()` corresponding to a higher order approximation with a relative
    precision given by `eps`. The DOFs are evaluated in physical coordinates
    returned by `eval_coors()`.
    """
    def _get_msd(iels, rx, ree):
        rvals = eval_dofs(iels, rx)
        rng = rvals.max() - rvals.min()
        n_components = rvals.shape[-1]
        msd = 0.0
        for ic in range(n_components):
            rval = rvals[..., ic]
            sd = rval[:, ree]
            # ~ max. second derivative.
            msd += nm.abs(sd[..., 0] + sd[..., 2]
                          - 2.0 * sd[..., 1]).max(axis=-1)
        msd /= n_components
        return msd, rng
    rx0 = ps.geometry.coors
    rc0 = ps.geometry.conn[None, :]
    rx, rc, ree = refine_reference(ps.geometry, 1)
    factor = rc.shape[0] / rc0.shape[0]
    iels = nm.arange(n_el)
    msd, rng = _get_msd(iels, rx, ree)
    eps_r = rng * eps
    flag = msd > eps_r
    iels0 = flag0 = None
    coors = []
    conns = []
    vdofs = []
    inod = 0
    for level in range(max_level + 1):
        if level < min_level:
            flag.fill(True) # Force refinement everywhere.
        elif level == max_level:
            # Last level - take everything.
            flag.fill(False)
        # Deal with finished elements.
        if flag0 is not None:
            ii = nm.searchsorted(iels0, iels)
            expand_flag0 = flag0[ii].repeat(factor, axis=1)
        else:
            expand_flag0 = nm.ones_like(flag)
        ie, ir = nm.where((flag == False) & (expand_flag0 == True))
        if len(ie):
            uie, iies = nm.unique(ie, return_inverse=True)
            # Each (sub-)element has own coordinates - no shared vertices.
            xes = eval_coors(iels[uie], rx0)
            des = eval_dofs(iels[uie], rx0)
            # Vectorize (how??) or use cython?
            cc = []
            vd = []
            for ii, iie in enumerate(iies):
                ce = rc0[ir[ii]]
                xe = xes[iie]
                cc.append(xe[ce])
                de = des[iie]
                vd.append(de[ce])
            cc = nm.vstack(cc)
            vd = nm.vstack(vd)
            nc = cc.shape[0]
            np = rc0.shape[1]
            conn = nm.arange(nc, dtype=nm.int32).reshape((nc // np, np))
            coors.append(cc)
            conns.append(conn + inod)
            vdofs.append(vd)
            inod += nc
        if not flag.any():
            break
        iels0 = iels
        flag0 = flag
        # Deal with elements to refine.
        if level < max_level:
            eflag = flag.sum(axis=1, dtype=bool)
            iels = iels[eflag]
            rc0 = rc
            rx0 = rx
            rx, rc, ree = refine_reference(ps.geometry, level + 2)
            msd, rng = _get_msd(iels, rx, ree)
            eps_r = rng * eps
            flag = msd > eps_r
    all_coors = nm.concatenate(coors, axis=0)
    conn = nm.concatenate(conns, axis=0)
    all_vdofs = nm.concatenate(vdofs, axis=0)
    mat_ids = nm.zeros(conn.shape[0], dtype=nm.int32)
    return level, all_coors, conn, all_vdofs, mat_ids