Subdomain integration with kernel (i.e., where each cells has a weight function)

General
1

@dokken @nate
Suppose I want to integrate on a subdomain while specifying the integration kernel (i.e., assigning a weight to every single cell in the marked subdomain) Where my subdomain is created in this way.

materials = MeshFunction("size_t", mesh, mesh.topology().dim()) # over cells

dx = Measure('dx', domain=mesh, subdomain_data = materials)
# interaction radius
r = 0.02

for cell_no in range(len(materials.array())):

	#define the cell
	cell = MeshEntity(mesh, 2, cell_no)

	#define the barycenter
	cell_bc = cell.midpoint()

	# create a subdomain class of radius r for each cell on the mesh
	class Interaction_range(SubDomain):
		def inside(self, x, on_boundary):
			return on_boundary and (x[0]-cell_bc[0])*(x[0]-cell_bc[0]) + (x[1]-cell_bc[1])*(x[1]-cell_bc[1]) <= r*r + tol

	# mark the interaction range
	subdomain_1 = Interaction_range()
	subdomain_1.mark(materials, 1)

	# integrate over the subdomain (removing the text function to get float results instead of a vector)
	nonlocal_interaction = assemble(0.5*dot(grad(f), grad(g))*dx(1))

I believe the line

nonlocal_interaction = assemble(0.5*dot(grad(f), grad(g))*dx(1))

Integrates the subdomain assigning an equal weight of value 1 to each cell in the subdomain by default. How can I specify a weight for each cell in the subdomain? i.e., if given a kernel K defined as


K =||midpoint of current cell||/r

Thank you

2

You can use x = SpatialCoordinate(mesh) to compute

midpoint = Constant(cell.midpoint())
r = sqrt((x[0]-midpoint[0])**2 + (x[1]-midpoint[1])**2 )

Cellmidpoint is not implemented in ufl, as far as I know, but you could use (x[0], x[1]) which is the coordinate of the integration point.

3

I know how to calculate the midpoints, the issue is how to specify a weight for each cell in the subdomain when using dx(tag) to integrate over the entire subdomain.

4

As I said, you can add it to

By defining a w Which could be defined with spatial coordinates and the midpoints as illustrated above, and then call
assemble(w*0.5*dot(grad(f), grad(g))*dx(1))

5

Okay. Thank you. What I mean is that if I have a different weight for say only cell index 0 in the subdomain while I want others (i.e., other cells in the subdomain) to maintain their initial weights. How can I make sure the w only acts on a single cell in the subdomain when applying assemble(w*0.5*dot(grad(f), grad(g))*dx(1))?

6

Choose w to be something that is unique in each cell (for instance a DG0 function, Which is a piecewise constant per cell). You can project any variable into DG 0, and it will then have a unique value per cell.

You could Also think of w as a function of some physical coordinates (x,y,z), here represented with SpatialCoordinate(mesh). When using numerical integration, x[i] is the ith component physical coordinate representing a quadrature point.

7

Thank you very much.

8

I tried this, but my answers were similar for the two kernels

from __future__ import print_function
from fenics import *
import numpy as np
import matplotlib.pyplot as plt


T = 100.0 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size

# Create mesh and define function space
nx = ny = 32
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)

# diffusivity contants
D_g = 0.0035 # m^2 per seconds
D_f = 0.0008 #m^2 per seconds

# fraction of volume space ocuupied by each specie 
w_g = 1.0
w_f = 1.0

# taxis
chi_f = 0.08



# Define function space for system of concentrations
P1 = FiniteElement('P', triangle, 1)
P2 = FiniteElement('P', triangle, 2)
P0 = FiniteElement('DG', triangle, 0)

element = MixedElement([P1, P1])
element2 = MixedElement([P2, P2])
V = FunctionSpace(mesh, element2)
V1 = FunctionSpace(mesh, P1)
V0 = FunctionSpace(mesh, P0)
V2 = FunctionSpace(mesh, P2)

# Define test functions
v_g, v_f = TestFunctions(V)

# Define functions
u = Function(V)

u_0 = Expression(('0.1',
				  '0.4*((0.5+0.5*tanh(20*x[0]-3))+(0.5+0.5*tanh(-20*x[0]-3)))'), degree=2)


u_n = interpolate(u_0, V)

# Split system functions to access components
g, f = split(u)
g_n, f_n = split(u_n)

# define the Neumann boundary conditions for g and f
Nbc_g = Constant(0) 
Nbc_f = Constant(0) 

# starting time
t = 0

# The LHS of the trapezoidal rule
a_b = (g - g_n)*v_g*dx + (f - f_n)*v_f*dx 

# interaction radius
R = 0.1

# defining the weak form
func_g = dot(D_g*grad(g), grad(v_g))*dx + D_g*Nbc_g*v_g*ds
func_f = dot(D_f*grad(f), grad(v_f))*dx - D_f*Nbc_f*v_f*ds 

# tolerance
tol = 10*DOLFIN_EPS

# to play with a mesh e.g, generate subdomain, we creat a MeshFunction
materials = MeshFunction("size_t", mesh, mesh.topology().dim()) # over cells f
bd = MeshFunction("size_t", mesh, mesh.topology().dim()-1) # over facets for boundary

# define the measures for the intending subdomain and its boundary
ds = Measure('ds', domain=mesh, subdomain_data=bd)
dx = Measure('dx', domain=mesh, subdomain_data = materials)

# no
Af_n = (0.5*w_f*f_n*chi_f)*grad(f_n) 

# define DG order 0 element, one value per cell since we want to assign a constant value to each cell
DG_Function_f1 = Function(V0) # for the two kernels
DG_Function_f2 = Function(V0)

# Functions to keep the weight of every cell in the subdomains
DG_S_f1 = Function(V0)
DG_S_f2 = Function(V0)

# create a subdomain class
class Interaction_range(SubDomain):
	def __init__(self, cell_bc, r, **kwargs):
		super().__init__(**kwargs)

		self.r = r
		self.cell_bc = cell_bc

	def inside(self, x, on_boundary):

		return (x[0]-self.cell_bc[0])*(x[0]-self.cell_bc[0]) +\
			 (x[1]-self.cell_bc[1])*(x[1]-self.cell_bc[1]) <= self.r*self.r + tol

# mark value for the subdomain
mark_value = 5

# a list of lists to store all the triangles in each sudomain generated by each triangles in the fine mesh globally
global_cells = []

# iterate over the triangles in the fine mesh
for cell in cells(mesh):
	
	#if cell.index()==0:

	#define the barycenterof the triangle
	cell_bc = cell.midpoint()
	
	# create subdomain class of radius r for each cell on the mesh
	subdomain = Interaction_range(cell_bc, R)

	# mark the triangles as zero
	materials.set_all(0)
	
	
	# mark all the triangles in the interaction region of the chosen triangle
	subdomain.mark(materials, mark_value)
	
	# keep the marked triangles in the interaction region of the chosen triangle in a list
	local_cells = list(np.flatnonzero(materials.array()==mark_value)) # get the triangles in the subdomain as a list
	global_cells.append(local_cells) # store the list of the triangles for each subdomain in a global list
		

# step into to time
for n in range(num_steps):
	
	# Update current time
	t += dt

	for cell in cells(mesh):# type: ignore

		# mark all as zero
		materials.set_all(0)

		cell_bc = cell.midpoint()

		# pick the appropriate subdomain (conisting of triangles) for each triangle by specifying the corresponding index and marking it
		materials.array()[global_cells[cell.index()]] = mark_value
		
		ele_list = np.where(materials.array()==mark_value) # get the index of the elements in the subdomain
		ele_list = np.array(ele_list).tolist() # to convert a turple to a list


		# a radius to make sure only the present cell remain in the sudomaim (checked and verified)
		R_1 = 0.05
		
		# loop throught the cells in each the subdoamin 
		for j in ele_list[0]: # picked the first element becuase I got a list of lists [[a, b]] and I need [a, b]
			
			sem = MeshEntity(mesh, 2, j).midpoint() # get the barycentre of each element in the subdomain
			
			# make each element a subdomain consisting of only the element
			subdomain_s = Interaction_range(sem, R_1)

			# mark the triangles as zero
			materials.set_all(0)
			
			# mark a single triangle
			subdomain_s.mark(materials, mark_value)

			# kernels of interaction 
			k_1 = 3/(pi*R*R)*(1 - np.linalg.norm(sem.array()-cell_bc.array())/R)
			
			s = 0.04 # std of the Gaussian kernel

			# Gaussian kernel
			k_2 = 1/(2*pi*s*s)*exp(-(1/(2*s*s))*(np.linalg.norm(sem.array()-cell_bc.array()))**2)

			# Carry out this expressions for each cell (subdomain) for each of the kernels
			DG_S_f1.vector()[j] = assemble(chi_f*0.5*(f_n*div(f_n*Af_n)*k_1*dx(mark_value)))
			DG_S_f2.vector()[j] = assemble(chi_f*0.5*(f_n*div(f_n*Af_n)*k_2*dx(mark_value)))

		# Step out to the main subdomain
		materials.set_all(0)	
		materials.array()[global_cells[cell.index()]] = mark_value
		

		# integrate over the marked subdomain 		
		nonlocal_interaction_f1 = DG_S_f1*dx(mark_value)
		nonlocal_interaction_f2 = DG_S_f2*dx(mark_value)

		DG_Function_f1.vector()[cell.index()]=assemble(nonlocal_interaction_f1)
		DG_Function_f2.vector()[cell.index()]=assemble(nonlocal_interaction_f2)

	# project it to a V1 space
	Nonlocal_interaction_f1=interpolate(DG_Function_f1, V1)
	Nonlocal_interaction_f2=interpolate(DG_Function_f2, V1)

	# get the weak form and add it to the rest of the weak form
	func_f+=Nonlocal_interaction_f1*v_f*dx
	
	# uncomment to use the second kerne;
	#func_f+=Nonlocal_interaction_f2*v_f*dx

	# Adding the functions (weak forms) together
	f_u = func_g + func_f
	
	# Euler 
	F = a_b + dt*f_u

	# solve 
	du = TrialFunction(V)
	J = derivative(F, u, du)
	Problem = NonlinearVariationalProblem(F, u, J = J)
	Solver  = NonlinearVariationalSolver(Problem)
	Solver.parameters['newton_solver']['convergence_criterion'] = 'incremental'
	Solver.parameters['newton_solver']['linear_solver'] = 'mumps'
	Solver.parameters['newton_solver']['relative_tolerance'] = 1.e-10
	Solver.parameters['newton_solver']['absolute_tolerance'] = 1.e-11
	Solver.parameters['newton_solver']['maximum_iterations'] = 100
	
	# Compute solution
	Solver.solve()

	# Extract the solutions from u
	gn_sol, fn_sol = u.split()
	
	# save the solutions for external view
	File('non_local_trial/trial_growth_factor.pvd') << (gn_sol,t)
	File('non_local_trial/trial_fibroblast.pvd') << (fn_sol,t)
	
	# Update previous solution
	u_n.assign(u)

Does it make sense to assign value to each cell in the subdomain by doing this

# loop throught the cells in each the subdoamin 
		for j in ele_list[0]: # picked the first element becuase I got a list of lists [[a, b]] and I need [a, b]
			
			sem = MeshEntity(mesh, 2, j).midpoint() # get the barycentre of each element in the subdomain
			
			# make each element a subdomain consisting of only the element
			subdomain_s = Interaction_range(sem, R_1)

			# mark the triangles as zero
			materials.set_all(0)
			
			# mark a single triangle
			subdomain_s.mark(materials, mark_value)

			# kernels of interaction 
			k_1 = 3/(pi*R*R)*(1 - np.linalg.norm(sem.array()-cell_bc.array())/R)
			
			s = 0.04 # std of the Gaussian kernel

			# Gaussian kernel
			k_2 = 1/(2*pi*s*s)*exp(-(1/(2*s*s))*(np.linalg.norm(sem.array()-cell_bc.array()))**2)

			# Carry out this expressions for each cell (subdomain) for each of the kernels
			DG_S_f1.vector()[j] = assemble(chi_f*0.5*(f_n*div(f_n*Af_n)*k_1*dx(mark_value)))
			DG_S_f2.vector()[j] = assemble(chi_f*0.5*(f_n*div(f_n*Af_n)*k_2*dx(mark_value)))