Zero extrapolation of a function to a greater domain

mesh
1

Hi,
I wonder how I can get a scalar product of functions defined on partially overlapping submeshes in the way that only common finite elements contribute into the integral (i.e. we assume that a function is equal to zero at points outside the submesh it was originally defined on)

from dolfin import *
# parameters['allow_extrapolation'] = True

mesh = UnitSquareMesh(32, 32)
submesh1 = UnitSquareMesh(16, 16)
submesh1.scale(.5)

submesh2 = UnitSquareMesh(16, 16)
submesh2.scale(.5)
submesh2.translate(Point(.25, 0))

W = FunctionSpace(mesh, 'P', 1)
W1 = FunctionSpace(submesh1, 'P', 1)
W2 = FunctionSpace(submesh2, 'P', 1)

w1 = Function(W1)
w2 = Function(W2)

w1 = project(w1, W)
w2 = project(w2, W)

dx = Measure('dx', domain=mesh)
w1 * w2 * dx

If I don’t specify parameters['allow_extrapolation'] = True, I get an error (which is reasonable).
If I specify it, I won’t get what I expect (it won’t be a zero extrapolation, kind of ‘nearest’ instead)

I have seen this work-around. Didn’t test it, still wonder if there is a “more fenics-like” way to do this.

I am doing GMsFEM implementation, by the way. So, maybe I have chosen the wrong approach for this problem. I will write about it in more detail later, in the comments.

EDITS:
Change 8 elements in submeshes (code part) to 16, because it is more logical in the problem articulation

2

Did you consider simply just integrating over the overlapping region?

from dolfin import *
import numpy as np
# parameters['allow_extrapolation'] = True

def overlap(mesh, mesh1, mesh2):
    '''Cell function of mesh marking mesh1 \cap mesh2'''
    meshes = (mesh, mesh1, mesh2)
    V, V1, V2 = (FunctionSpace(m, 'DG', 0) for m in meshes)

    chi1, chi2 = Function(V), Function(V)
    LagrangeInterpolator.interpolate(chi1, interpolate(Constant(1), V1))
    LagrangeInterpolator.interpolate(chi2, interpolate(Constant(1), V2))
    # Intersect
    chi = chi1.vector().get_local() * chi2.vector().get_local()
    chi = np.fromiter(chi, dtype='uintp')
    # As cell function
    f = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
    f.array()[:] = chi

    return f

mesh = UnitSquareMesh(32, 32)
submesh1 = UnitSquareMesh(8, 8)
submesh1.scale(.5)

submesh2 = UnitSquareMesh(8, 8)
submesh2.scale(.5)
submesh2.translate(Point(.25, 0))

W = FunctionSpace(mesh, 'P', 1)
W1 = FunctionSpace(submesh1, 'P', 1)
W2 = FunctionSpace(submesh2, 'P', 1)

w1 = interpolate(Constant(1), W1)
w2 = interpolate(Constant(2), W2)

w1W, w2W = Function(W), Function(W)
LagrangeInterpolator.interpolate(w1W, w1)
LagrangeInterpolator.interpolate(w2W, w2)

chi = overlap(mesh, submesh1, submesh2)

dx = Measure('dx', domain=mesh, subdomain_data=chi)
print(assemble(w1 * w2 * dx(1)))
1 Like