I fear my question was not clear, so I have tried to reformulate it in a simpler setting, independent from any previous notebook or application.
I attach a MWE below, at least showing concretely what I am trying to do.
In essence I define a problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) and necessary arguments a, L, bcs.
Now I need to solve it in a loop where one argument form, say a is changing, but cannot call fem.petsc.LinearProblem each iteration as it will crash.
So I would like to take the fem.petsc.LinearProblem defintion out of the loop, and modify said form “in place” this I do not know how to do (it is easy to achieve for fem.constant as in the Poisson tutorial, assigning a new value each loop with fem.constant.value, without need to recompile the problem).
In the example below I define a list of forms, and would like to iterate the linear problem solution over this list.
MWE
########## DOLFINX ##############
%matplotlib notebook
from dolfinx import *
import matplotlib.pyplot as plt
import numpy as np
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from mpi4py import MPI
from dolfinx.fem import FunctionSpace, Function
from dolfinx.fem import VectorFunctionSpace
from dolfinx.fem import petsc
from dolfinx import fem, io
from dolfinx import mesh
from dolfinx.mesh import CellType, create_rectangle
import ufl
import pyvista
pyvista.set_jupyter_backend("trame")
L = 4
domain = create_rectangle(MPI.COMM_WORLD, np.array([[0,0],[L,1]]), [30,30], cell_type=CellType.triangle)
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
# Problem parameters
thetamoy = 0.4 # target average material density
E = fem.Constant(domain, ScalarType(1))
nu = fem.Constant(domain, ScalarType(0.3))
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
f = fem.Constant(domain, ScalarType([0,-1])) # vertical downwards force
# Function space for density field
V0 = fem.FunctionSpace(domain, ("DG", 0))
# Function space for displacement
V2 = fem.VectorFunctionSpace(domain, ("CG", 2))
def clamped_boundary(x):
return np.isclose(x[0], 0)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0,0], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V2, fdim, boundary_facets), V2)
def left(x):
return np.isclose(x[0], 0)
def right(x):
return np.isclose(x[0], L)
fdim = domain.topology.dim -1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)
# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
def eps(v):
return ufl.sym(ufl.grad(v))
def sigma(v):
return (lamda*ufl.nabla_div(v)*ufl.Identity(2)+2*mu*eps(v))
# return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)
def energy_density(u, v):
return ufl.inner(sigma(u), eps(v))
# Inhomogeneous elastic variational problem
# u_ = ufl.TestFunction(V2)
# du = ufl.TrialFunction(V2)
# a = ufl.inner(sigma(u_), eps(du)) * ufl.dx
# L = ufl.dot(f, u_)*ds(2)
w = ufl.TestFunction(V2)
u = ufl.TrialFunction(V2)
a = ufl.inner(sigma(w), eps(u)) * ufl.dx
L = ufl.dot(f, w)*ds(2)
problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh= problem.solve()
form_list = [ k*ufl.inner(sigma(w), eps(u)) * ufl.dx for k in [1,2,3]]
for i in range(3):
#### modify a in place, e.g. a = form_list[i]
uh = problem.solve()
Getting into the documentation, I came up with this approach
from dolfinx.fem.forms import form as _create_form
form_list = [ k*ufl.inner(sigma(w), eps(u)) * ufl.dx for k in [1,2,3]]
for i in range(3):
#### modify a in place, e.g. a + form_list[i]
# The next line does not work
# problem._a = form_list[i]
## Tried
problem._a = _create_form(form_list[i])
uh = problem.solve()
print(uh.x.array[0])
seems create_form goes from ufl.forms to fem.forms.
Is this the way to do it? Thanks a lot