As f and g are finite element functions this will return
\int_{\Omega_1} 0.5\nabla f \cdot \nabla g~\mathrm{d}x, a scalar value, where \Omega_1\subset \Omega such that \Omega_1 = \{k\in cells(\Omega)\vert marker(k)=1 \}
As non_local_interaction is a scalar value, you can add it to any variational form.
Dear @dokken, if I need to do:
materials = MeshFunction("size_t", mesh, mesh.topology().dim()) # over cells
dx = Measure('dx', domain=mesh, subdomain_data = materials)
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))
In C++ how can I go about it?
I actually need to carry out this subdomain integration either at every time step or 10-time steps depending on my choice but it’s expensive to loop through all the elements at every time step for nx=ny=256
Or is there a faster way to iterate through all the elements in materials without a loop?
Do you mean how to code this in C++?
There are many C++ demos illustrating how dolfin works. I do not have the time to translate your code to another language
Yes, to code it in C++ either still within python (where we do it in quotes) or otherwise, perhaps you can give me a link? or is there a faster way to do it without involving C++?
Are you moving your mesh between time steps? If not the interaction range is constant, and does not need to recompute the range at each time step.
You would however have to store the range for each cell, but that can be done with numpy arrays, and not full meshfunction objects.
The mesh is constant at each time step. How do I store the interaction range for each cell using NumPy arrays?
I tried
import numpy as np
arr = np.array([])
for cell_no in range(len(materials.array())): # type: ignore
#define the cell
cell = MeshEntity(mesh, 2, cell_no)
#define the barycenterof the cell
cell_bc = cell.midpoint()
# create subdomain class of radius r for each cell on the mesh
class Interaction_range(SubDomain):
def inside(self, x, on_boundary):
return (x[0]-cell_bc[0])*(x[0]-cell_bc[0]) + (x[1]-cell_bc[1])*(x[1]-cell_bc[1]) <= r*r + tol
'''Interaction_range = CompiledSubDomain('(x[0]-cell_bc_0)*(x[0]-cell_bc_0) + (x[1]-cell_bc_1)*(x[1]-cell_bc_1) <= r*r + tol', cell_bc_0=cell_bc[0],cell_bc_1=cell_bc[1], tol=tol, r=r)'''
# mark the interaction range
subdomain_1 = Interaction_range()
subdomain_1.mark(materials, 1)
arr = np.append(arr, subdomain_1)
but it returns as
<__main__.Interaction_range object at 0x7fc70c274d00>
How do I specify dx(1) for each cell after collecting these interaction ranges for each cell when integrating for each time step? I apologize if my questions are a bit dumb.
Here is an improved approach.
from dolfin import *
import numpy as np
import time
r = 0.1
tol = 10*DOLFIN_EPS
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
N = 35
mesh = UnitSquareMesh(N, N)
start = time.perf_counter()
value = 3
data = np.array([], dtype=np.int32)
offsets = [0]
materials = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for cell in cells(mesh):
#define the barycenterof the cell
cell_bc = cell.midpoint()
# create subdomain class of radius r for each cell on the mesh
subdomain = Interaction_range(cell_bc, r)
materials.set_all(0)
subdomain.mark(materials, value)
local_cells = np.flatnonzero(materials.array()==value)
data = np.hstack([data, local_cells])
offsets.append(len(data))
end = time.perf_counter()
print(f"Generate radius {end-start:.2e}")
start = time.perf_counter()
for cell in cells(mesh):
materials.set_all(0)
materials.array()[data[offsets[cell.index()]:offsets[cell.index()+1]]] = value
# Pick a random cell to visualize
if cell.index() == 3*N:
print(cell.midpoint()[0], cell.midpoint()[1])
File("mat.pvd") << materials
end = time.perf_counter()
print(f"Assign marker for each domain {end-start:.2e}")
however it does scale with the number of elements, which is not good.
Generate radius 2.43e+01
0.49523809523809526 0.047619047619047616
Assign marker for each domain 3.10e-02