# Initial conditions for a mixed space

Hi everyone,

I’ve been working on implementing a code to solve a reaction-diffusion system and then move the mesh. However, I’ve been encountering some strange results, and I’m almost certain that the issue lies in the definition of the initial conditions for the mesh. I’ve set up a mixed function space for the coordinates of the mesh and the curvature. I suspect that the problem may be in this part of the code.

``````V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())

normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)

X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)

``````

Specifically, I suspect the issue lies in `X0.sub(0).interpolate(x_exp)` . I’ve already attempted to interpolate the mesh coordinate into `X0.sub(0)` , but the mesh always seems to move in a strange manner.

What is the correct way to define this function in the initial step?

Here is my MWE

``````# +
import os

from mpi4py import MPI
from petsc4py import PETSc

import numpy as np
import math

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, log, mesh
from dolfinx.fem import Function, functionspace, Expression, form
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, GhostMode
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner, div,dot

# Save all logging to file
log.set_output_file("log.txt")
# -

try:
import gmsh  # type: ignore
except ImportError:
import sys
print("This demo requires gmsh to be installed")
sys.exit(0)

def gmsh_circle(model: gmsh.model, name: str) -> gmsh.model:
"""Create a Gmsh model of a circle.

"""
model.setCurrent(name)

gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", 0.8)

circle = model.occ.addCircle(0, 0, 0, 1, tag=1)

# Synchronize OpenCascade representation with gmsh model
model.occ.synchronize()

# Add physical marker for cells. It is important to call this

# Generate the mesh
model.mesh.generate(dim=1)
model.mesh.setOrder(2)

return model

def create_mesh(comm: MPI.Comm, model: gmsh.model, name: str, filename: str, mode: str):
"""Create a DOLFINx from a Gmsh model and output to file.
"""
msh, ct, ft = gmshio.model_to_mesh(model, comm, rank=0,gdim=2)
msh.name = name
ct.name = f"{msh.name}_cells"
ft.name = f"{msh.name}_facets"
with XDMFFile(msh.comm, filename, mode) as file:
msh.topology.create_connectivity(0, 1)
file.write_mesh(msh)
file.write_meshtags(ct, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")
file.write_meshtags(ft, msh.geometry, geometry_xpath=f"/Xdmf/Domain/Grid[@Name='{msh.name}']/Geometry")

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)

# Create model
model = gmsh.model()
model = gmsh_circle(model, "circle")
model.setCurrent("circle")

create_mesh(MPI.COMM_SELF, model, "circle", f"out_gmsh/mesh_rankCircle_{MPI.COMM_WORLD.rank}.xdmf", "w")

with XDMFFile(MPI.COMM_WORLD, "out_gmsh/mesh_rankCircle_0.xdmf", "r") as xdmf:
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

gmsh.finalize

P1 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)
ME = functionspace(mesh, mixed_element([P1, P1], gdim=2))

n = ufl.CellNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
r = x/ufl.sqrt(ufl.dot(x, x))
sign = ufl.sign(ufl.dot(n, r))
n_oriented = sign*n

V1 = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,), gdim=2)
V2 = element("Lagrange", mesh.basix_cell(), 2, gdim=2)

ME2 = functionspace(mesh, mixed_element([V1, V2],gdim=2))

# Trial and test functions of the space `ME` are now defined:

q, v = ufl.TestFunctions(ME)

Xi, phi = ufl.TestFunctions(ME2)

# ```{index} split functions
# ```

# +
u = Function(ME)  # current solution
u0 = Function(ME)  # solution from previous converged step

# X = Function(ME2)  # current solution
X0 = Function(ME2)  # solution from previous converged step
X0.name = "X0"

(Xh, H) = ufl.TrialFunctions(ME2)  # current solution
#(Xh0, H0) = ufl.TrialFunctions(ME2)  # solution from previous converged step

# Split mixed functions
u1, u2 = ufl.split(u)
u10, u20 = ufl.split(u0)

# # Split mixed functions
# Xh, H = ufl.split(X)
Xh0, H0 = ufl.split(X0)

au = 0.9
bu = 1.1
av = 0.8
bv = 0.95

# Interpolate initial condition

u.sub(0).interpolate(lambda x: (bu-au)*np.random.rand(x.shape[1]) +au)
u.sub(1).interpolate(lambda x: (bv-av)*np.random.rand(x.shape[1]) +av)

u.x.scatter_forward()

V01, _ = ME2.sub(0).collapse()
V02, _ = ME2.sub(1).collapse()

x_exp = Expression(x, V01.element.interpolation_points())
H_exp = Expression(div(sign*n), V02.element.interpolation_points())

normal_expr = Expression(n_oriented, V01.element.interpolation_points())
normal_vec = Function(V01)

X0.sub(0).interpolate(x_exp)
X0.sub(1).interpolate(H_exp)

X0.x.scatter_forward()

print (f"X array  {X0.sub(0).x.array[:]}")

# -

#X0.x.array[:] = uh.x.array

u_vec = ufl.as_vector([u1, u2])
kp = ufl.as_vector([1.5 , 0])

# Next, various model parameters are defined:

dt = 5.0e-05  # time step

# Weak statement of the equations

k = dt
d = 10
gamma = 29
a = 0.1
b = 0.9
tol = 1e-4

d1 = 1.0
d2 = d

ks = 25
kb = 2

-(gamma*(u1*u1*u2-u1+a))*q*dx \
-(gamma*(-u1*u1*u2+b))*v*dx

ax = dot(((Xh) / k),n_oriented)*phi*dx + kb*inner(grad(H),grad(phi))*dx  - 1/2 *kb* H0**2 * H *phi*dx + ks*H*phi*dx\
Lx =  dot(((Xh0) / k),n_oriented)*phi*dx + dot(kp,u_vec)*phi*dx

# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver

ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()

X_h = Function(ME2)
X_h.name = "X_h"

problemX = LinearProblem(ax, Lx, u = X_h)

# Step in time
t = 0.0

# Output file
name = "mesh_at_t"+str(t)
mesh.name = name

#  Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 3 * dt
else:
T = 5 * dt

# Get the sub-space for c and the corresponding dofs

normal_vec.interpolate(normal_expr)

V0, dofs = ME.sub(0).collapse()
W02, dofs2 = ME2.sub(0).collapse()

u1 = u.sub(0)
u2 = u.sub(1)

Xh0 = X0.sub(0)
H0 = X0.sub(1)

file1 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/outputu1.xdmf", "w")
file3 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/normal.xdmf", "w")
file4 = XDMFFile(MPI.COMM_WORLD, "Schnakenberg/H.xdmf", "w")

file1.write_mesh(mesh)
file3.write_mesh(mesh)
file4.write_mesh(mesh)

file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file4.write_function(X0.sub(1), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")

norm_stop = 0.1e-6

fileL2 = open('L2.dat', 'w')

z = 0
inc = 0
impresion = 1

#print (f"MESH array {mesh.geometry.x[:,]}")

while (t < T):
t += dt

problemX.solve()
r = solver.solve(u)
inc += 1

print(f"Step {int(t/dt)}: num iterations: {r[0]}")

l2_norm = np.linalg.norm(u.x.array[dofs]-u0.x.array[dofs])/dt
#print (f"L2_norm = {l2_norm}")

fileL2.write(f"{t}\t{l2_norm}\n")
u0.x.array[:] = u.x.array
X0.x.array[:] = X_h.x.array

num=inc/impresion-1

# print (f"MESH array shape {mesh.geometry.x[:,].shape}")
# print (f"X array shape {X_h.x.array[dofs2].shape}")
# print (f"MESH array {mesh.geometry.x[:,]}")
# print (f"X array  {X_h.x.array[dofs2]}")
#print (f"X array  {X0.sub(0).x.array[:]}")

X_coord = X0.x.array[dofs2].reshape(len(mesh.geometry.x[:,1]),2)
DeltaX = mesh.geometry.x[:,[0,1]]-X_coord

#print (f"Delta X array  {DeltaX}")
mesh.geometry.x[:,0] = X_coord[:,0]
mesh.geometry.x[:,1] = X_coord[:,1]

normal_vec.interpolate(normal_expr)

if (int(num)==z) or (int(num-z)==0):
us, vs = ufl.split(u)

# Save solution to file (VTK)

name = "mesh_at_t"+str(t)
mesh.name = name

file1.write_mesh(mesh)
file1.write_function(u1, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file3.write_mesh(mesh)
file3.write_function(normal_vec, t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
file4.write_mesh(mesh)
file4.write_function(X0.sub(0), t, mesh_xpath=f"/Xdmf/Domain/Grid[@Name='{mesh.name}']")
z += 1

# if l2_norm < norm_stop:
#     print("l2_norm convergence")
#     break

file1.close()

fileL2.close()

print(l2_norm)
#print(H.x.array)

``````

Anyway, if you could find out if the error lies in another part, I would appreciate it,

Thank you

This is not a safe way to set displacements.
See for instance:

which shows how this should be done with mixed spaces.