Project matrix from one function space to another function space

Linear Algebra
1

I have to function spaces, like

mesh1 = UnitSquareMesh(8, 8)
mesh2 = refine(mesh1)

V1 = FunctionSpace(mesh1, "CG", 1)
V2 = FunctionSpace(mesh2, "CG", 1)

project_matrix = something_like_project(V2, V1) # I want to obtain the matrix

So how to implement?

2

I’m not quite sure what you mean, but it sounds like maybe you’re looking for something like the PETScDMCollection functionality. Take a look at the following example:

from dolfin import *

# Two different meshes and function spaces:
mesh1 = UnitIntervalMesh(7)
mesh2 = UnitIntervalMesh(100)
V1 = FunctionSpace(mesh1,"Lagrange",1)
V2 = FunctionSpace(mesh2,"Lagrange",1)

# Create some function in V1:
f1 = Function(V1)
f1.interpolate(Expression("x[0]*x[0]",degree=2))

# Create a transfer matrix from V1 to V2, and use it to interpolate f1 in V2:
A = PETScDMCollection.create_transfer_matrix(V1,V2)
f2 = Function(V2)
f2.vector()[:] = A*f1.vector()

# Plot results:
import matplotlib.pyplot as plt
plot(f2)
plt.show()
3 Likes
3

Thank you very much! That’s exactly what I want.

4

Sorry to bother you again. I find something is wrong if I obtain mesh2 from refining mesh1, like below.

from dolfin import *

# Two different meshes and function spaces:
mesh1 = UnitIntervalMesh(7)
# mesh2 = UnitIntervalMesh(100)
mesh2 = refine(refine(mesh1)) # ???
V1 = FunctionSpace(mesh1,"Lagrange",1)
V2 = FunctionSpace(mesh2,"Lagrange",1)

# Create some function in V1:
f1 = Function(V1)
f1.interpolate(Expression("x[0]*x[0]",degree=2))

# Create a transfer matrix from V1 to V2, and use it to interpolate f1 in V2:
A = PETScDMCollection.create_transfer_matrix(V1,V2)
f2 = Function(V2)
f2.vector()[:] = A*f1.vector()

# Plot results:
import matplotlib.pyplot as plt
plot(f2)
plt.show()

And the result is like

figure_1
figure_1.png800×600 26.9 KB

And I believe it’s about DOFs numbering, but how to fix it?

5

You can try using the solution from this old thread

https://fenicsproject.org/qa/1327/plotting-high-order-function-in-1d/

of sorting the mesh coordinates:

# ...

# mesh2 = UnitIntervalMesh(100)
mesh2 = refine(refine(mesh1)) # ???

####### Added #######
mesh2.coordinates().sort(axis=0)

V1 = FunctionSpace(mesh1,"Lagrange",1)
V2 = FunctionSpace(mesh2,"Lagrange",1)

# ...
6

Thank you for your time and patience. This one solves my problem, but I find another problem.
When I add
mesh2.coordinates().sort(axis=0),
my stiffness matrix become singular. I think this function may have a bug.
The whole code is below

from fenics import *
import matplotlib.pyplot as plt
import sympy


def AllExpressions(sigma, cell): # the PDE's exact solution and source term
    x, y = sympy.symbols("x[0] x[1]")
    u = sympy.exp(x*y) * sympy.sin(pi*x) * sympy.sin(pi*y) # exact solution, satisfies boundary condition

    u_x = u.diff(x, 1)
    u_y = u.diff(y, 1)
    u_xx = u_x.diff(x, 1)
    u_yy = u_y.diff(y, 1)
    f = -(u_xx + u_yy) \
        + (sigma.values()[0]*u_x + sigma.values()[1]*u_y)

    u = Expression(sympy.printing.ccode(u), degree=6, cell=cell)
    f = Expression(sympy.printing.ccode(f), degree=4, cell=cell)
    return u, f


def solve_coarse_problem():
    # ------------------------------------------------- coarse mesh problem -------------------------------------------
    p, q = TrialFunction(VH), TestFunction(VH)
    AH = assemble(inner(grad(p), grad(q)) * dx + inner(sigma, grad(p)) * q * dx)
    bH = assemble(f * q * dx)
    bcH.apply(AH)
    bcH.apply(bH)

    solver = KrylovSolver()
    solver.set_operator(AH)
    uH = Function(VH)
    solver.solve(uH.vector(), bH)

    plt.figure()
    plot(uH, mode="warp", title="solution on coarse space")


def solve_fine_problem():
    # ----------------------------------------------- fine mesh problem -----------------------------------------------
    u, v = TrialFunction(Vh), TestFunction(Vh)
    Ah = assemble(inner(grad(u), grad(v)) * dx + inner(sigma, grad(u)) * v * dx)
    bh = assemble(f * v * dx)
    bch.apply(Ah)
    bch.apply(bh)

    solver = KrylovSolver()
    solver.set_operator(Ah)
    uh = Function(Vh)
    solver.solve(uh.vector(), bh)

    plt.figure()
    plot(uh, mode="warp", title="solution on fine space")


if __name__ == '__main__':
    sigma = Constant((8, 8))
    cell = "triangle"
    u_exact, f = AllExpressions(sigma, cell)
    element = FiniteElement("CG", cell, 1)

    coarse = UnitSquareMesh(8, 8)
    fine = refine(refine(refine(coarse)))
    fine.coordinates().sort(axis=0) # comment this line, I can get the right result!!!

    VH = FunctionSpace(coarse, element)
    bcH = DirichletBC(VH, Constant(0.0), DomainBoundary())

    Vh = FunctionSpace(fine, element)
    bch = DirichletBC(Vh, Constant(0.0), DomainBoundary())


    solve_coarse_problem()
    solve_fine_problem()
    plt.show()
7

It’s really just a hack to fix interactive plotting on refined meshes in 1D. In 2D or 3D, or if plotting 1D results with an external software like Paraview, it’s not necessary.

8

I am not interested in interactive plotting. I want to get a transfer matrix from coarse space to fine space by coarse_to_fine = PETScDMCollection.create_transfer_matrix(VH, Vh). After refining coarse mesh several times, I get the fine mesh and I define fine function space with it.
Now if I add fine.coordinates().sort(axis=0), I get a singular stiffness matrix. But without fine.coordinates().sort(axis=0), the stiffness matrix is not singular and I get the right numerical solution.
So where is wrong? Thank you very much in advance!

9

The issue is that renumbering of vertices due to coordinates sort is not reflected elsewhere, for example in encoding of cells. In the following the sort introduces cells with 0 volume

from dolfin import *

mesh = UnitSquareMesh(32, 32)
print(min(c.volume() for c in cells(mesh)))

mesh.coordinates().sort(axis=0)
print(min(c.volume() for c in cells(mesh)))