You don’t need to instantiate all these objects every iteration…
import dolfinx
import dolfinx.io
from mpi4py import MPI
import ufl
from ufl import inner, grad, dx
import numpy as np
def set_bc(V):
u_bc = dolfinx.Function(V)
bc_facets = np.where(
np.array(dolfinx.cpp.mesh.compute_boundary_facets(mesh.topology)) == 1)[0]
bc_dofs = dolfinx.fem.locate_dofs_topological(V,
mesh.topology.dim-1,
bc_facets)
with u_bc.vector.localForm() as loc:
loc.setValues(bc_dofs, np.full(len(bc_dofs), 0))
return dolfinx.DirichletBC(u_bc, bc_dofs)
mesh = dolfinx.UnitSquareMesh(MPI.COMM_WORLD, 4, 4)
mesh.topology.create_connectivity(mesh.topology.dim-1,
mesh.topology.dim)
Hcurl = ufl.FiniteElement("Nedelec 1st kind H(curl)", mesh.ufl_cell(), 1)
H1 = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
N = 1000
V = dolfinx.FunctionSpace(mesh, H1)
f = dolfinx.Function(V)
V = dolfinx.FunctionSpace(mesh, Hcurl)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
i_plus_one = dolfinx.Constant(mesh, 1.0)
a_p = inner(i_plus_one*u, v) * dx
L_p = inner(v, grad(f)) * dx
solution = dolfinx.fem.LinearProblem(a_p, L_p,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
for i in range(N):
print(f"{i+1}/{N}")
i_plus_one.value = float(i + 1)
f.interpolate(lambda x: (i+1) * x[0]**2 + 0*x[1])
f.x.scatter_forward()
u = solution.solve()