On 30 May 2016, at 15:24, Daniel Ruprecht <D.Ruprecht@leeds.ac.uk> wrote:
Hi Lawrence,
Thanks for the info - I do have two independent meshes (the first variant you describe), but the problem is not fully coupled: so yes, the sequence would be for every time step
1. Solve Problem A 2. Get Boundary Data of Problem A 3. Prescribe Neumann Data for Problem B based on Boundary Data from A (*) 4. Solve B
(*) We do not need to actually couple the meshes in terms of identifying DOFs or so - it would be enough to some how create an expression from the boundary data of A and feed this into the boundary condition for B.
So at the moment the idea is to derive from Expression and add an additional argument which can be a Function on mesh A. Then, in eval, I access values on points of mesh B through .at(X, dont_raise=True).
OK, so you effectively are evaluating (a function of) function_on_A at the nodes points of mesh B that live on the boundary? I would probably do this like so: meshA = ... meshB = ... V_a = FunctionSpace(meshA, element_a) V_b = FunctionSpace(meshB, "CG", 1) solution_A = Function(V_a) # Assume that meshB has boundary marker 1 on the appropriate boundary boundary_nodes = DirichletBC(V_b, 0, 1).nodes # Get the global coordinates of the boundary nodes by interpolating the # coordinate field into the vector-valued version of the space V_b c = Function(Vector_V_b) c.interpolate(SpatialCoordinate(meshB)) Xs = c.dat.data_ro[boundary_nodes, :] boundary_data_B = Function(V_b) # solve problem A ... # transfer to domain B. # You'll probably want to filter out the `None`s. boundary_data_B.dat.data[boundary_nodes] = solution_A.at(Xs, dont_raise=True) # solve problem B ... Note that .at is currently collective, so this won't work right in parallel. This is sort of what eval in the Expression is doing, but restricted to the boundary, and "vectorised" so it should be a little faster. Lawrence
I can't use a single mesh because the mesh of problem A will move in time - only in a simple way, so we can always identify where the two meshes would touch, but this still prevents us from doing everything on a single mesh.