How to interpolate the gradient of numerical solution into another space？

General
1

Hello all,

I am trying to solve a multi-domain problem. For example, there are two cube domains(domain 1 and domain 2), the interface of two domains is z=0 plane.

I solve one standard Poisson equation on domain 1, the weak formulation is (grad(u), grad(v))=(f, v). The numerical solution on domain 1 is u.

then, I want to solve another Poisson equation on domain 2 by weak formulation: (grad(u_s), grad(v_s))=(f_s, v_s)+(grad(u), grad(v_s)){interface}, where, (grad(u), grad(v_s)){interface} is boundary integral on the interface, u is the numerical solution on domain 1.

my question is how to coding this boundary integral term: (grad(u), grad(v_s))_{interface} ?

If this boundary integral term is (u, v_s)_{interface}, I can solve it by interpolation, see MWE. but if this term have grad, I do not know how to solve it?

Thanks!

Here is MWE

import matplotlib.pyplot as plt
from dolfin import *
from mshr import *
import numpy as np

''' Poisson-1 (this part is standard Poisson equation code, demo_poisson.py)'''

mesh = UnitCubeMesh(16, 16, 16)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)',degree=2)
g = Expression('sin(5*x[0])',degree=2)
a = inner(grad(u), grad(v))*dx
L = f*v*dx + g*v*ds


# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Save solution in VTK format
file = File("poisson.pvd")
file << u






'''Poisson-2'''

domain_s = Box(Point(0,0,-1), Point(1,1,0))

mesh_s = generate_mesh(domain_s, 32)


# define interface space

class BoundaryX1(SubDomain):
    def inside(self, x, on_boundary):
        return  near(x[2], 0, DOLFIN_EPS)
    
interface_s = BoundaryX1()

boundary_markers = MeshFunction('size_t', mesh_s, 2, mesh_s.domains())
interface_s.mark(boundary_markers,0)
ds_s = Measure('ds', domain=mesh_s, subdomain_data=boundary_markers)


exteriormesh_s = BoundaryMesh(mesh_s, "exterior")
interfacemesh_s=SubMesh(exteriormesh_s,interface_s)
V_interface_s = FunctionSpace(interfacemesh_s, "Lagrange", 1)

# interpolate u obtained by Poisson-1 to Iu_s

Iu_s=Function(V_interface_s)
Iu_s.set_allow_extrapolation(True)
Iu_s.interpolate(u)



# solve Poisson-2

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary_s(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS 

V_s = FunctionSpace(mesh_s, "Lagrange", 1)

# Define boundary condition
u0_s = Constant(1.0)
bc_s = DirichletBC(V_s, u0_s, boundary_s)

u_s = TrialFunction(V_s)
v_s = TestFunction(V_s)
f_s = Constant(1.0)
a_s = inner(grad(u_s), grad(v_s))*dx
#L_s= Iu_s*v_s*ds_s(0)+f_s*v_s*dx                # this works
L_s= dot(grad(Iu_s),grad(v_s))*ds_s(0)+f_s*v_s*dx    # this fails

u_s = Function(V_s)
solve(a_s == L_s, u_s, bc_s)

file = File("poisson_s.pvd")
file << u_s
1 Like
2

Interpolate u from one mesh to the other, and then compute the gradient in the ufl form

3

Hi dokken,

Thanks for the reply.
In MWE, I Interpolate u from one mesh to interface space, obtained Iu_s, and compute the gradient of Iu_s in ufl form ,
L_s= dot(grad(Iu_s),grad(v_s))ds_s(0)+f_sv_s*dx

The error produced is

------------------- Start compiler output ------------------------
/var/folders/zy/2nk6sy516r99j63dk36ym15r0000gn/T/tmp7xnvek0r/ffc_form_40a208a52c346d0bbc3f6296b3fef794f578a15d.cpp:182:18: error: redefinition of 'J_c4'
    const double J_c4 = coordinate_dofs[1] * FE15_C0_D01_F_Q1[0][0][0] + coordinate_dofs[7] * FE15_C0_D01_F_Q1[0][0][1];
                 ^
/var/folders/zy/2nk6sy516r99j63dk36ym15r0000gn/T/tmp7xnvek0r/ffc_form_40a208a52c346d0bbc3f6296b3fef794f578a15d.cpp:176:18: note: previous definition is here
    const double J_c4 = coordinate_dofs[2] * FE15_C0_D01_F_Q1[0][0][0] + coordinate_dofs[5] * FE15_C0_D01_F_Q1[0][0][1];
                 ^
/var/folders/zy/2nk6sy516r99j63dk36ym15r0000gn/T/tmp7xnvek0r/ffc_form_40a208a52c346d0bbc3f6296b3fef794f578a15d.cpp:184:18: error: redefinition of 'J_c5'
    const double J_c5 = coordinate_dofs[1] * FE15_C0_D01_F_Q1[0][0][0] + coordinate_dofs[10] * FE15_C0_D01_F_Q1[0][0][1];
4

I would not use submeshes, as they dont work very well. Have you tried using MeshView to define the mesh? Or simply interpolate on to the full mesh (with extrapolation)?

6

Sorry, I didn’t use MeshView，simply interpolate on to the full mesh with extrapolation.

In MWE, if
L_s= dot(grad(Iu_s),grad(v_s))ds_s(0)+f_s v_sdx
is replaced by
L_s= Iu_sv_sds_s(0)+f_sv_s*dx

It works.

7

This is what you have in your code, which does not match what you say in:

8

Sorry, my reply is not clear. I define interface mesh and build interface space on this mesh. Then, I interpolate u into interface space, obtained Iu_s, and set Iu_s allow extrapolation.

9

And that is what Im suggesting you shouldn’t do.
Either

  1. Make a MeshView of the boundary, and used the Mixed assembly functionality described in
    https://dl.acm.org/doi/abs/10.1145/3471138 (And various forum posts
  2. use interpolation with extrapolation on the whole other mesh.
10

Thank you for your suggestions.

I’m confused about the second suggestion. How to interpolate with extrapolation on the whole other mesh? I try to interpolate u to the whole other mesh directly, but failed, so I have to interpolate u into interface mesh.
Could you present any minicode to learn this method?

11

I am not at a computer (Christmas holiday), so its hard for me to give a functional code.
Pseudocode

mesh_a = …
V_a = FunctionSpace(mesh_a, …)
u_a = Function(V_a)
u_a.set_allow_extrapolation(True)

mesh_b = …
V_b = FunctionSpace(mesh_b, …)
u_b = Function(V_b)
u_b.interpolate(u_a)
12

Thank you. I got it！ Have a nice Christmas holiday.