### Interpolating/Projecting From One Mesh to Another

99
views
0
7 weeks ago by
I have a system of PDE's where one of the 5 equations, c_e, is defined in a 2D space (doesn't have to be this one exactly, but takes the same shape):

and the remaining 4 equations are coupled by the solution at the top boundary (c_e,s). How can I interpolate the 1D solution c_e,s onto the main solution mesh which is also 1D? My main mesh has a varied point distribution, but the interval mesh with 280 points from 0 to 3 is pretty close. My strategy was to try and define c_e,s as a  subdomain of the c_e mesh first to get my 1D variable, then re-interpolate the solution onto a subset of my main mesh. The subset of the main mesh is [0,1],[2,3] to match the x axis of the 2D mesh. But in the step going from c_es defined on the x axis of the 2D mesh to the main 1D submesh everything falls apart. This is my example 2D function:

This is after finding the solution at the top of the mesh (expected result):

And this is after trying to project onto the subset of my main mesh:

I expect my final result to *appear* to be identical to the solution for the top of the 2D mesh, but not even the magnitude is right.

Here's my MWE:

from mshr import *
import numpy as np
import dolfin as fin
import matplotlib.pyplot as plt

# Create 2D mesh
r1 = Rectangle(fin.Point(0,0), fin.Point(1,1))
r2 = Rectangle(fin.Point(2,0), fin.Point(3,1))
domain = r1 + r2
domain.set_subdomain(1, r1)
domain.set_subdomain(3, r2)
mesh = generate_mesh(domain, 18)

# Find values along the top
bmesh = fin.BoundaryMesh(mesh, 'exterior')
cc = fin.MeshFunction('size_t', bmesh, bmesh.topology().dim())
top = fin.AutoSubDomain(lambda x: x[1] > 1-1e-8)
top.mark(cc, 1)
submesh = fin.SubMesh(bmesh, cc, 1)

# Main mesh, 1D
dmesh = fin.IntervalMesh(280, 0, 3)

# Sub mesh to project onto, derived from the main mesh
dcc = fin.MeshFunction('size_t', dmesh, dmesh.topology().dim())
elect = fin.AutoSubDomain(lambda x: x[0] <= 1+1e-8 or x[0] >= 2-1e-8)
elect.mark(dcc, 5)
dsubmesh = fin.SubMesh(dmesh, dcc, 5)
#fin.plot(dsubmesh)
#plt.show()

V = fin.FunctionSpace(dsubmesh, 'Lagrange', 1)
W = fin.FunctionSpace(mesh, 'Lagrange', 1)
W_1 = fin.FunctionSpace(submesh, 'Lagrange', 1)

u_1 = fin.interpolate(fin.Expression('x[0]', degree=1), W)
#fin.plot(u_1)
#plt.show()

us = fin.interpolate(u_1, W_1)
us.set_allow_extrapolation(True)
#plt.plot(submesh.coordinates()[:][:,0], us.vector().get_local()[fin.vertex_to_dof_map(W_1)], 'o')
#plt.show()

us_v = fin.interpolate(us, V)
#plt.plot(dsubmesh.coordinates()[:][:,0], us_v.vector().get_local()[fin.vertex_to_dof_map(V)], 'o')
#plt.show()

Additionally, this sometimes returns an error (probably has something to do with the randomness of the mshr mesh):

*** Error:   Unable to create mesh entity.
*** Reason:  Mesh entity index -1 out of range [0, 22] for entity of dimension 1.
*** Where:   This error was encountered inside MeshEntity.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset:  unknown

Thanks in advance for any help! I'm running FEniCS 2017.2.0 on Ubuntu 16.04.

Community: FEniCS Project
I think I figured out how to do what I was looking for, I just don't know how valid this overall approach is. So I'd love to hear any feedback. I took the submesh and normalized it to zero before interpolating:
submesh.coordinates()[:][:, 1] -= submesh.coordinates()[:][:, 1]​
and removed:
us.set_allow_extrapolation(True)​
written 7 weeks ago by Chris Macklen
After a lot more searching, I think what I'm really trying to do is map the dof's along the top surface of the 2D mesh to a function in my 1D mesh. That way I can run the solve command an keep everything updated. Is this correct? And if so, how do I implement it? I can't seem to find the right verbiage to find my answer.
written 7 weeks ago by Chris Macklen