Interpolate on boundary

General
1

Hello!
I’m using dolfinx 0.8.0 and I need to solve 3 problems.

Say I solve u on mesh1 and v on mesh2 and that both domains share a side \Gamma.

Now I want to solve a third problem on mesh2, using as Neumann data on \Gamma the difference between u and v.

How can I do that? because I have the feeling that by using create_non_matching_meshes_data I’m computing much more than needed (as I only need to interpolate the data on one boundary of the domain, which indeed is a box)

On top of that, I don’t understand how to correctly call create_non_matching_meshes_data in dolfinx 0.8.0.

Before I used to do:

create_non_matching_meshes_data(u._V, v._V)

passing the function spaces. Now I did

create_non_matching_meshes_data(mesh_to=mesh2, element=("Lagrange",1), mesh_from = mesh1)

But it didn’t work either. I seems to me that there should be a FiniteElement constructor class somwhere but I can’t find it.

Thank you very much!

2

Could you make a sketch of your domain, clarifying if:

  1. the boundaries are conforming, ie the facets on domain 1 align with the facets of domain 2?

Add your minimal attempt at solving this problem. One could for instance simply interpolate some know data into a function on mesh 1, and say that one want to use it in a Function on the other mesh.

Also note that you can restrict the number of cells in create_nonmatching_interpolation_data to be only those adjacent to the boundary Gamma on mesh2:
https://docs.fenicsproject.org/dolfinx/v0.8.0/python/generated/dolfinx.fem.html#dolfinx.fem.create_nonmatching_meshes_interpolation_data

There are unfortunately some bugs in 0.8.0 that has resulted in a complete refactoring of this code in 0.9.0/main
https://docs.fenicsproject.org/dolfinx/main/python/generated/dolfinx.fem.html#dolfinx.fem.Function.interpolate_nonmatching

3

The meshes are not conformanto on \Gamma since they are created separately by two calls to gmsh python addBox. One of them is simply a box so it could be entirely meshed in dolfinx, but the other one is not.

Taking that into account is it still possible to do it in 0.8.0 or may be I should export the boundary data to numpy and back?

I’ll add the m.nw.e later as I’m writting from home.

4

I would use non-matching interpolation with the restriction of cells close to the boundary, which you can get with dolfinx.mesh.compute_incident_entities(mesh, bndry_facets, mesh.topology.dim-1, mesh.topology.dim)

5

Hi, thank you for the answer but I’m not sure if I understood how to use the return of compute_incident_entities.

This would be the main non-working example:

from dolfinx.mesh import create_box
from dolfinx.mesh import locate_entities, meshtags
from mpi4py import MPI 
import numpy as np 

L = 1
N = 5

mesh_1 = create_box(comm=MPI.COMM_WORLD,points=[(0,0,0),(L,L,L)], n=[N,N,N])

Gamma= 1
Gamma_back = 2
Gamma_sides = 3

boundaries = [( Gammma, lambda x: np.isclose(x[0], 0)),
              ( Gamma_back, lambda x: np.isclose(x[0], L) ),
              ( Gamma_sides, lambda x: np.isclose( np.maximum(np.abs(x[1]), np.abs(x[2])), L) ) ] 

facet_indices, facet_markers = [], []
fdim = mesh_1.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh_1, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags_1 = meshtags( mesh_1 , fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

mesh_2 = create_box(comm=MPI.COMM_WORLD,points=[(-Lx,0,0),(0,Ly,Lz)], n=[Nx,Ny,Nz])


boundaries = [( Gammma, lambda x: np.isclose(x[0], 0)),
              ( Gamma_back, lambda x: np.isclose(x[0], -L) ),
              ( Gamma_sides, lambda x: np.isclose( np.maximum(np.abs(x[1]), np.abs(x[2])), L) ) ] 


facet_indices, facet_markers = [], []
fdim = mesh_2.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh_2, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags_2 = meshtags( mesh_2 , fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

# for creating ds1 and ds1 measures later

from dolfinx.fem import Function, FunctionSpace

H1 = FunctionSpace( mesh_1, ("Lagrange",1))
H2 = FunctionSpace( mesh_2, ("Lagrange",1))

u1 = Function(H1)

# ... solution of a variational formulation for u1

u2 = Function(H2)
u2.interpolate(u1,...) # here I would like to interpolate u1 restricted to Gamma, on u2 restricted to Gamma.
# I don't need to interpolate on the rest of the boundary.
6

I tried to strip the “example” from everything that is not relevant and now it looks as a domain decomposition problem. It is not.

7

Consider the minimal working example:

from mpi4py import MPI
import dolfinx
import numpy as np

L = 1
N = 5


def interface(x):
    return np.isclose(x[0], L)


mesh_1 = dolfinx.mesh.create_box(
    comm=MPI.COMM_WORLD, points=[(0, 0, 0), (L, L, L)], n=[N, N, N]
)
V = dolfinx.fem.functionspace(mesh_1, ("Lagrange", 2))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: x[0] + np.cos(x[1]) + np.sin(np.pi * x[2]))


Nx, Ny, Nz = 12, 13, 15
mesh_2 = dolfinx.mesh.create_box(
    comm=MPI.COMM_WORLD, points=[(L, -0.1, -0.1), (2 * L, 3 * L, 3 * L)], n=[Nx, Ny, Nz]
)

V_2 = dolfinx.fem.functionspace(mesh_2, ("Lagrange", 2))

boundary_2 = dolfinx.mesh.locate_entities_boundary(
    mesh_2, mesh_2.topology.dim - 1, interface
)
adjacent_cells = dolfinx.mesh.compute_incident_entities(
    mesh_2.topology, boundary_2, mesh_2.topology.dim - 1, mesh_2.topology.dim
).astype(np.int32)

interpolation_data = dolfinx.fem.create_nonmatching_meshes_interpolation_data(
    mesh_2.geometry, V_2.element, mesh_1, adjacent_cells, padding=1.0e-6
)


u2 = dolfinx.fem.Function(V_2)
u2.interpolate(u, cells=adjacent_cells, nmm_interpolation_data=interpolation_data)

with dolfinx.io.VTXWriter(mesh_1.comm, "u_1.bp", [u]) as bp:
    bp.write(0.0)


with dolfinx.io.VTXWriter(mesh_2.comm, "u_2.bp", [u2]) as bp:
    bp.write(0.0)

image
image1081×628 82.4 KB

8

Thank you very much for the example. I don’t have a VTX reader in paraview so I tried changing your code:

with dolfinx.io.VTXWriter(mesh_1.comm, "u_1.bp", [u]) as bp:
    bp.write(0.0)


with dolfinx.io.VTXWriter(mesh_2.comm, "u_2.bp", [u2]) as bp:
    bp.write(0.0)

with

from dolfinx.io import XDMFFile

with XDMFFile(mesh_1.comm, "u1.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_1)
    xdmf.write_function(u)

with XDMFFile(mesh_2.comm, "u2.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh_2)
    xdmf.write_function(u2)

However, when I try to run it ( dolfinx 0.8.0 ) I get the error (which I didn’t get before modifying it):

Traceback (most recent call last):
  File "/home/manuel/Documentos/MultiPhysics/src/MWE.py", line 48, in <module>
    xdmf.write_function(u)
  File "/home/manuel/anaconda3/envs/env3-10-complex/lib/python3.10/site-packages/dolfinx/io/utils.py", line 252, in write_function
    super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)
RuntimeError: Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?

This is weird, as I was not expecting the writing of function u to fail… maybe the u2 but not the u. Is this one of those bugs in 0.8.0 you mentioned earlier? (I hope not -.-’ )

Again, thank you very much for the code!

9

Forget it, I understood that it was a problem of XDMF with (“Lagrange”,2) elements. I changed it for (“Lagrange”,1) and it works fine.

Thank you very much!

10

You should also be able to use VTKFile with second order Lagrange;)