Change integration points for projection of DG solution

Question

Hi,

In an attempt to visualize a P>1 DG solution I wanted to project my (coarse) DG solution onto a finer DG mesh of P=1 I encountered the problem that the fine DG elements that touch the coarse DG element-interfaces do not project well. I use the following code.

Say we have a DG solution 'u', on a unit-square mesh 'mesh' of size (nrElsW x nrElsH) corresponding to functionspace 'V'. Then:

class evalCoarseSol(Expression):
    def eval(self, values, x):
        values[0] = self.coarseSol(x)
    def value_shape(self):
        return ()

meshVisualize = UnitSquareMesh(nrElsW*10,nrElsH*10)
VVisualize = FunctionSpace(meshVisualize, "DG", 1)
p = TrialFunction(VVisualize)
q = TestFunction(VVisualize)
uEval = evalCoarseSol(element=V.ufl_element()) #The problem lies here
uEval.coarseSol = u

#L2 projection:
aVisualize = p*q*dx
LVisualize = uEval*q*dx
uVisualize = Function(VVisualize)
solve(aVisualize == LVisualize, uVisualize)

If i look at which points are evaluated for solving the integration of the L2 weak form, then these turn out to lie on the triangle boundaries. For fine triangles touching the coarse-element interfaces this gives the wrong value and is a problem. Why is this happening, these are not the correct Gauss-quadrature points right? If it would use the 'normal' Gauss points that lie in the triangle surface/volume my problem should disappear.

How to solve this? I feel like i should pass a different keyword argument to the Expression constructor to make it use different integration points.

Thanks,
Stein

Answer 1

Did you try to use built-in project in python? I am not sure whether project works for DG functions in different meshes with DG functions.

Let's simplify the problem. Suppose we have u1 in mesh1 (coarse) which is DG (degree k), and u2 in mesh2 (e.g., uniform refinement of mesh1) which is also DG (degree k). Then how can we obtain the correct interpolation/projection in FEniCS? i.e.,

interpolate u1 to u1_2 in mesh2
project u2 to u2_1 in mesh1

As far as I know, direct application u1_2 = interpolate(u1, V2) where V2 = FunctionSpace(mesh2, 'DG', k) does not work. Basically, the problem is similar to your situation, the function evaluation of DG function at interface causes a problem (it enforce the same value at that point). By the way, if we check u1(p) where p lies on the interface, then it gives us a value (although this should not be true since u1 is not well-defined at this point).

In paper, interpolate are straightforward (just go cell by cell in mesh1 and compute the values in children cells). But is there already a method to do this job in FEniCS?

Comments

I did originally try the built-in Project function, which gave me the same issue. That's why I wrote out the L2 projection weak formulation instead. In the end it did indeed seem simplest to write a custom interpolation loop that runs over all elements, which is what I did (code added below if someone is interested).

However the essence of my initial question still remains:
Why is FEniCS using integration points on the boundary of the triangle at all?
That doesn't make sense to me and worries me a little. Standard Gauss quadrature points are located in the triangle interior. What is FEniCS doing? I am also surprised by how difficult it is to change the integration rules. Maybe I am missing something.

Here is the code:

from dolfin import *
import mpl_toolkits.mplot3d as a3
import matplotlib
import matplotlib.colors as colors
import pylab as pl
import numpy as np


def plotDG(u,refinefactor=3,exactSol=None,showedges = False,plot=True):
    V = u.function_space()
    mesh = V.mesh()
    PolDeg = V.ufl_element().degree()
    if PolDeg > 1:
        for i in range(refinefactor):
            mesh = refine(mesh)
        V = FunctionSpace(mesh, "DG", 1)


    dofmap = V.dofmap()
    gdim = mesh.geometry().dim()
    dofs = dofmap.dofs()
    dofs_x = V.tabulate_dof_coordinates().reshape((-1, gdim))

    ax = a3.Axes3D(pl.figure())
    cmap = matplotlib.cm.get_cmap('jet')
    minz = min(u.vector().array())
    maxz = max(u.vector().array())
    for c in cells(V.mesh()):
        c_dofs = dofmap.cell_dofs(c.index())

        XY_mp = c.midpoint()
        z_mp = u(XY_mp)
        vec_mp = np.array([XY_mp.x(),XY_mp.y()])
        zcolor = cmap(  (z_mp-minz)/(maxz-minz)  )

        vtx = np.zeros((3,3))
        for i,d in enumerate(c_dofs):
            x,y = dofs_x[d][:]
            vec = np.array([x,y])
            z = u(Point(vec))
            vec2 = vec+0.001*(vec_mp-vec)
            z2 = u(  Point(vec2)  )
            if abs(z-z2) > 0.001:
                z = z2
                vec = vec2
            if not exactSol == None:
                z = exactSol(vec)-z

            vtx[i][:] = [x,y,z]

        tri = a3.art3d.Poly3DCollection([vtx])
        tri.set_color(colors.rgb2hex(zcolor) )
        if showedges:
            tri.set_edgecolor('k')
        ax.add_collection3d(tri)

    if plot:
        pl.show()

That makes a matplotlib render of the individual facets. For higher order bases it subdivides the triangles to capture the smoothness. This function is meant for quick inspection, as he plot is very slow and due to the limited 3d-capabilties of Matplotlib the depth of the facets is often a little buggy. I also wrote a little script that outputs a vtk file, which can be read in paraview. That looks very nice:

def exportDG(u,name,refinefactor=5,exactSol=None):
    print "Exporting solution"
    V = u.function_space()
    mesh = V.mesh()
    PolDeg = V.ufl_element().degree()
    if PolDeg > 1:
        for i in range(refinefactor):
            mesh = refine(mesh)
        V = FunctionSpace(mesh, "DG", 1)


    dofmap = V.dofmap()
    gdim = mesh.geometry().dim()
    dofs = dofmap.dofs()
    dofs_x = V.tabulate_dof_coordinates().reshape((-1, gdim))


    ptxyz = ""
    connec = ""
    offsets = ""
    types = ""
    dataval = ""
    ct = 0
    for c in cells(V.mesh()):
        c_dofs = dofmap.cell_dofs(c.index())

        XY_mp = c.midpoint()
        z_mp = u(XY_mp)
        vec_mp = np.array([XY_mp.x(),XY_mp.y()])
        for i,d in enumerate(c_dofs):
            x,y = dofs_x[d][:]
            vec = np.array([x,y])
            z = u(Point(vec))
            vec2 = vec+0.00001*(vec_mp-vec)
            z2 = u(  Point(vec2)  )
            if abs(z-z2) > 0.0000001:
                z = z2
                vec = vec2
            if not exactSol == None:
                z = exactSol(vec)-z

            ptxyz+=" "+str(x)+" "+str(y)+" "+str(z)
            dataval+= " "+str(z)
            connec+= " "+str(ct)
            ct+=1
        offsets+=" "+str(ct)
        types+=" 5"

    vtktxt = """<?xml version="1.0"?>
<VTKFile type="UnstructuredGrid"  version="0.1"  >
<UnstructuredGrid>
<Piece  NumberOfPoints=" """+str(ct)+""" " NumberOfCells=" """+str(ct/3)+""" ">
<Points>
<DataArray  type="Float64"  NumberOfComponents="3"  format="ascii">"""+ptxyz+"""</DataArray>
</Points>
<Cells>
<DataArray  type="UInt32"  Name="connectivity"  format="ascii">"""+connec+"""</DataArray>
<DataArray  type="UInt32"  Name="offsets"  format="ascii">"""+offsets+"""</DataArray>
<DataArray  type="UInt8"  Name="types"  format="ascii">"""+types+"""</DataArray>
</Cells>
<PointData  Scalars="f_31"> 
<DataArray  type="Float64"  Name=" """+u.name()+""" "  format="ascii">"""+dataval+"""</DataArray> 
</PointData> 
</Piece>
</UnstructuredGrid>
</VTKFile>"""
    f = file(name+'.vtu','w')
    f.write(vtktxt)
    f.close()