Please I have a similar question but in my case, I intend to find the divergence of U since U is a vector.
Does it make sense to do this
U = as_vector((Ux, Uy))
where Ux and Uy are both scalar values depending on x and/or y.
Such that I can do this
Div(U)
Please make the problem reproducible, ie. add a code that can reproduce what you want to do.
I cannot guess what you want to do with what you have added in your post.
This is a minimal example
I have this equation to solve
\frac{\partial u_{_1}}{\partial t} = D_{_{1}}\Delta u_{_1} - \lambda u_{_1},\\ \frac{\partial u_{_2}}{\partial t} = \nabla \cdot (D_{_{2}}\nabla u_{_2} - \mu u_{_2}A) ,
Where
A(\mathbf{x}) = \frac{1}{R}\int\limits_{\mathbf{B}(0,R)}\displaystyle\mathbf{n}(\mathbf{y})(G(u(\mathbf{x+y}, t))\, \text{d}y_{_{1}} \text{d}y_{_{1}}, is the nonlocal interactions term,
G(u(\mathbf{x} + \mathbf{y}, t)) = u_{_1}(\mathbf{x}+\mathbf{y},t){(1-\rho (u))}^+,
\mathbf{B}(0,R):= \{ \zeta \in \mathbb{R}^{^{2}} : \|\zeta\|_{_{2}}\leq R\} is a closed ball of radius R (denoting the sensing region), centred at 0, and \mathbf{n(y)} denotes the unit radial vector originating from \mathbf{x} and moving towards \mathbf{x} + \mathbf{y}\in \mathbf{B}(0,R) for any \mathbf{y}\in \mathbf{B}(0,R) and is defined as follows:
\mathbf{n(y)} := \begin{cases}
\frac{\mathbf{y}}{\|\mathbf{y}\|_{_{2}}}, & \text{if } \mathbf{y}\in \mathbf{B}(0,R)\setminus \{(0, 0)\}\\
(0, 0), & \text{otherwise, }\end{cases}
and
To define this nonlocal interaction, I defined a subdomain region as follows
# 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)
# define DG order 0 element, one value per cell since we want to assign a constant value to each cell
DG_Function_f = Function(V0)
DG_Function_m = 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 = []
brc_e=[]
# iterate over the triangles in the fine mesh
for cell in cells(mesh):
#if cell.index()==0:
#define the barycenter of 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)
#print(f"elements marked {materials.array}")
# 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
R_1 = 0.05
#iterate on the elements of the mesh
for cell in cells(mesh):
# 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 elemens in the subdomain
ele_list = np.array(ele_list).tolist()
R_1 = 0.05
sum_w = []
ele_list[0].remove(cell.index())
for j in ele_list[0]:
sem = MeshEntity(mesh, 2, j).midpoint() # get the barycentre of each elemet in the subdomain
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)
local_cells = list(np.flatnonzero(materials.array()==mark_value))
I got my normal vector \mathbf{n(y)} this way
# get the normal vector
nd_x, nd_y, nd_z = (sem.array()-cell_bc.array())/np.linalg.norm(sem.array()-cell_bc.array())
Calculated my A component-wise this way though
Au_x = (1/R)*assemble(nd_x*Gu*dx(mark_value))
Au_y = (1/R)*assemble(nd_y*Gu*dx(mark_value))
# Since A should be a function of $x$ since we integrated with respect to $y$
x = SpatialCoordinate(mesh)
# store the x and y values of the non-local term as a vector
Au = as_vector((Au_x*x[0], Au_y*x[1]))
materials.set_all(0)
materials.array()[global_cells[cell.index()]] = mark_value
# calculate the integral on each cell as a constant value
DG_Function.vector()[cell.index()] = assemble(mu*(div(u2*Au))*dx(mark_value))
I then project this result to a V1 space then joined it with the rest of the weak form. So ultimately,
should I have done this?
Au = as_vector((Au_x*x[0], Au_y*x[1]))
I really hope you respond. I keep getting this error,
dijitso failed to load existing file:
/home/fenics/.cache/fenics/dijitso/lib/libdijitso-ffc_form_c716e16b7bf831ef1d725a817e5d4873838e3825.so
error is:
/home/fenics/.cache/fenics/dijitso/lib/libdijitso-ffc_form_c716e16b7bf831ef1d725a817e5d4873838e3825.so: cannot apply additional memory protection after relocation: Cannot allocate memory
For the projection you need to do something like:
StressSpace = TensorFunctionSpace(mesh, "DG", 0) S_func = project(S, StressSpace)
Could you please explain why we should use “DG” and 0 here?
Can we use StressSpace = TensorFunctionSpace(mesh, "CG", 1) here?
Thank you very much in advance.
If your displacement is in a space of piecewise continuous functions, its derivative lives in the space of cell-wise constants. Thus projecting into a CG-1 space will give rise to Gibb’s phenomenon, see for instance: Projection and interpolation — FEniCS Tutorial @ Sorbonne and references therein for more information
Your reply is very clear and helpful. I understand now. Thank you very much.