Fast interpolation between two meshes

Asked by Artur Palha on 2013-03-04

I have a mesh and I solve a problem with FEniCS using that mesh.

Every time step I need to evaluate the solution obtained with FEniCS in a specific set of points (large). I can use the function interpolate but I was thinking that since my points are always the same if I could use something similar to assemble. This assemble would correspond to generating a matrix with the evaluation of each basis functions at each point where I wish to interpolate. Then, interpolation would be obtained just by multiplying my solution vector by this matrix.

Is this possible?

Thank you.

-artur palha

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Mikael Mortensen
Solved:
2013-03-06
Last query:
2013-03-06
Last reply:
2013-03-05

Hi

Den Mar 4, 2013 kl. 2:41 PM skrev Artur Palha:

> New question #223355 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/223355
>
> I have a mesh and I solve a problem with FEniCS using that mesh.
>
> Every time step I need to evaluate the solution obtained with FEniCS in a specific set of points (large). I can use the function interpolate but I was thinking that since my points are always the same if I could use something similar to assemble. This assemble would correspond to generating a matrix with the evaluation of each basis functions at each point where I wish to interpolate. Then, interpolation would be obtained just by multiplying my solution vector by this matrix.
>
> Is this possible?

If I understand you correctly it sort of is, but it's not part of the basic FEniCS package. I do turbulence simulations and typically need to probe one or many locations many times during a simulation. I use the pdesys (https://launchpad.net/cbcpdesys) package that has functionality for efficiently probing the same point many times. The code basically does what a Function eval does, but precomputes as much as possible. Also, it works in parallel since the probes are only created when the point is found on the local mesh. There is a fast C++ version of the code in branch lp:~mikael-mortensen/dolfin/probe, and a slower Python version in cbc.cfd.tools.Probe.py that you obtain by installing pdesys.

If you install CBC.PDESys the following demo will interpolate a solution onto a structured 2D mesh each timestep (you could alternatively use another FEniCS mesh to achieve the same)

from cbc.cfd.tools.Probe import *

# Create a solution Function that we want to probe many times. Any mesh will do
mesh = UnitSquareMesh(20, 20)
V = VectorFunctionSpace(mesh, 'CG', 1)
f0 = Expression(("sin(pi*x[0]*t)", "cos(pi*x[1]*t)"),
                element=V.ufl_element(), t=0)

# Create a structured mesh on a rectangle [0.25, 0.75] x [0.25, 0.75]
N = 20
xx = linspace(0.25, 0.75, N)
xx = xx.repeat(N).reshape((N, N))
x = zeros((N*N, 3))
for i in range(N):
    for j in range(N):
        x[i*N + j, 0 ] = xx[j, i] # x-value
        x[i*N + j, 1 ] = xx[i, j] # y-value
# Note that we could also use another FEniCS mesh here instead, but that only works in serial

# Create probes for all points in x
num_evals = 10
#probes = CppProbes(x, V) # Fast C++ version
probes = Probes(x, V, num_evals) # Slower Python version

# Evaluate the solution Function f inside a loop and probe it
for j in range(num_evals):
    # Modify the solution Function
    f0.t = j*0.5
    f = interpolate(f0, V)
    # Evaluate all probes
    probes(f)

# Check the first probe on each processor
probe_0 = probes[0]
print "Global index and position of probe = ", probe_0[0], probe_0[1].coordinates()
print "Values f[0] = ", probe_0[1].get_probe(0)
print "Values f[1] = ", probe_0[1].get_probe(1)

if MPI.num_processes() == 1:
    # plot the solution if running with just one proc
    from scitools.std import surfc
    z = zeros((N, N, num_evals))
    for i in range(N):
        for j in range(N):
            val = probes[i*N + j][1].get_probe(0)
            z[i, j, :] = val[:]

    plot(f[0])
    surfc(xx.transpose(), xx, z[:, :, num_evals-1], indexing='xy')

Best regards

Mikael

>
> Thank you.
>
> -artur palha
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Artur Palha (artur-palha) said : #2

Thank you Mikael, I will give a try on that.

-artur palha

Artur Palha (artur-palha) said : #3

Thanks Mikael Mortensen, that solved my question.