Interpolating/Projecting From One Mesh to Another
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-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 <= 1+1e-8 or x >= 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', 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.