Hi Lawrence,

We set up the problem according to your info and it works initially. However, as we move meshA, at the some point all the results from 
solution_A.at(Xs, dont_raise=True)
become None's and we have absolutely no idea why this is happening - it works for about 5 time steps and then it doesn't anymore. I attached a small example reproducing the problem and a short video showing what is happening.

To illustrate, we print out

solution_A.at((0.0, 2.2), dont_raise=True)
and the point is marked in the video. For the first three steps, the result is the correct 5.0 (an arbitrary value we assign to the solution on meshA at the moment) but after that it becomes None, even though the point is still clearly ``in contact'' with meshA.

Any ideas what is happening here?

Cheers, Daniel

On 30/05/16 16:04, Lawrence Mitchell wrote:

      

        
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.


      

      

      

      

      

      
_______________________________________________
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