How to define piecewise source term in FEniCSx?

General
1

For example, in heat PDEs, I need to define a piecewise source term in a 3D cuboid as shown in Fig 1.
and the python script is as follows:

z1 = 0
z2 = t_bulk

def Omega_left(x):
	return ((x[0] <= 0.5 * length) & ((x[2] >= z1) & (x[2] <= z2) ))

def Omega_right(x):
	return ((x[0] >= 0.5 * length) & ((x[2] >= z1) & (x[2] <= z2) ))

cells_left = locate_entities(mesh, mesh.topology.dim, Omega_left)
cells_right = locate_entities(mesh, mesh.topology.dim, Omega_right)

f.x.array[cells_left] = np.full_like(cells_left, 1e10, dtype=default_scalar_type)
f.x.array[cells_right] = np.full_like(cells_right, 1e11, dtype=default_scalar_type)

with XDMFFile(MPI.COMM_WORLD, "f.xdmf", "w") as xdmf:
	xdmf.write_mesh(mesh)
	xdmf.write_function(f)

However, the source term function visualized by paraview is shown in Fig. 2.
We can see that there is a gap with value 0 between the interface of the two piecewise areas.
And things get worse when more complex piecewise source term is used as shown in Fig. 3.
Therefore, how can we define piecewise source term without these undesired gaps? Thanks!!!

FEniCSx提问
FEniCSx提问1252×737 47.7 KB

2

Add a tolerance to deal with floating point here, ie.

def Omega_left(x, tol=1e-14):
	return ((x[0] <= 0.5 * length+tol) & ((x[2] >= z1-tol) & (x[2] <= z2+tol) ))

def Omega_right(x,tol=1e-14):
	return ((x[0] >= 0.5 * length-tol) & ((x[2] >= z1-tol) & (x[2] <= z2+tol) ))
3

problem solved. Thank you very much for your timely reply!!!