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). 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. Cheers, Daniel On 30/05/16 15:13, Lawrence Mitchell wrote:
Dear Daniel,
On 30 May 2016, at 13:29, Daniel Ruprecht <D.Ruprecht@leeds.ac.uk> wrote:
Dear Firedrakers,
I am trying to solve two coupled problems on different domains (essentially a small box next to a larger box), where the solution at the boundary from one problem determines the Neumann boundary condition for the second problem at the interface.
I am not sure how to efficiently use the solution of the first problem to prescribe the boundary condition for the second. My idea was to derive from the Expression class and overwrite the eval function including a couple of if statements to make sure the coordinates are in the region where both meshes actually match. This will work, but as you note, is very slow. I presume the identification of the degrees of freedom on the boundary that you need to move from one mesh to the other is static (i.e. you only need to do so once). If that is the case you could do this up front and then (assuming you solve in A, transfer to B, solve in B) you could do something like:
solve in A
boundary_function_in_B.dat.data[boundary_nodes_B] = solution_in_A.dat.data_ro[boundary_nodes_in_A]
solve in B
Is this the correct way to do it (it seems quite slow, probably because the derived expression can't be compiled in C) or is there a cleverer way to do this? How are you describing the coupled problem. As two completely separate meshes, or as a single mesh in which you mask out the equations in one half or the other?
Moreover, is the coupling of the form:
solve in domain A
transfer solution to domain B
solve in domain B using the solution from A as a boundary condition
Or do you solve the full coupled A-B problem in one go.
The current approach in firedrake for solving in subdomains is currently a bit hacky, and what I describe now only works if you have a single mesh.
What you can do is as follows. Mark the two regions of your domain (firedrake as of earlier last week correct reads cell subdomains from Gmsh files).
Now you can define an indicator function in a DG0 functionspace
I_A = Function(DG0) I_B = Function(DG0)
which you can initialise using a par_loop to be 1 in the appropriate part of the domain.
You can use these indicator functions in your forms to only select contributions from the appropriate domain in the appropriate parts.
e.g. you can use:
dx_A = I_A*dx dx_B = I_B*dx ds_A = I_A*ds ds_B = I_B*ds
To apply the neumann boundary term you'll need to use an interior facet integral. This requires a little more thought as to how to restrict the indicator functions, but is not too bad.
You'll then want to mask out the appropriate parts of the operators in the two domains with strong boundary conditions in the "other" half. You can do this by subclassing DirichletBC and overriding the "nodes" property to select all function space nodes that are in the wrong domain.
A better approach to this latter problem is on our roadmap (although we do not currently have a timescale for it). See https://github.com/firedrakeproject/firedrake/issues/723 for details. Please comment if you have ideas on what should be available.
Note that coupling two *unrelated* meshes is a harder problem https://github.com/firedrakeproject/firedrake/issues/714 effectively because the operator mapping between meshes (even if they are topologically conforming) is almost always not the identity.
Hope this is helpful!
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-- Dr. Daniel Ruprecht University Academic Fellow School of Mechanical Engineering, Room 4.50 University of Leeds Leeds LS2 9JT, UK Email: d.ruprecht@leeds.ac.uk Phone: +44 (0)113 343 22 01