FEniCS
diff --git a/‎python/demo/demo_axis.py‎
Lines changed: 1 addition & 0 deletions b/‎python/demo/demo_axis.py‎
Lines changed: 1 addition & 0 deletions
diff --git a/‎python/demo/demo_biharmonic.py‎
Lines changed: 1 addition & 0 deletions b/‎python/demo/demo_biharmonic.py‎
Lines changed: 1 addition & 0 deletions
diff --git a/‎python/demo/demo_navier-stokes.py‎
Lines changed: 34 additions & 59 deletions b/‎python/demo/demo_navier-stokes.py‎
Lines changed: 34 additions & 59 deletions
diff --git a/‎python/demo/demo_pml.py‎
Lines changed: 1 addition & 0 deletions b/‎python/demo/demo_pml.py‎
Lines changed: 1 addition & 0 deletions
diff --git a/‎python/demo/demo_poisson.py‎
Lines changed: 1 addition & 0 deletions b/‎python/demo/demo_poisson.py‎
Lines changed: 1 addition & 0 deletions
diff --git a/‎python/demo/demo_scattering_boundary_conditions.py‎
Lines changed: 1 addition & 0 deletions b/‎python/demo/demo_scattering_boundary_conditions.py‎
Lines changed: 1 addition & 0 deletions
@@ -676,6 +676,7 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):
         },
     )
     Esh_m = problem.solve()
+    assert isinstance(Esh_m, fem.Function)
     assert problem.solver.getConvergedReason() > 0, "Solver did not converge!"
 
     # Scattered magnetic field
 
@@ -228,6 +228,7 @@
 
 problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
 uh = problem.solve()
+assert isinstance(uh, fem.Function)
 
 # The solution can be written to a  {py:class}`XDMFFile
 # <dolfinx.io.XDMFFile>` file visualization with ParaView or VisIt
 
@@ -181,7 +181,6 @@
 
 import ufl
 from dolfinx import default_real_type, fem, io, mesh
-from dolfinx.fem.petsc import apply_lifting, assemble_matrix, assemble_vector, set_bc
 
 try:
     from petsc4py import PETSc
@@ -195,6 +194,7 @@
     print("This demo requires petsc4py.")
     exit(0)
 
+from dolfinx.fem.petsc import LinearProblem
 
 if np.issubdtype(PETSc.ScalarType, np.complexfloating):
     print("Demo should only be executed with DOLFINx real mode")
@@ -307,7 +307,6 @@ def jump(phi, n):
 a -= ufl.inner(p, ufl.div(v)) * ufl.dx
 a -= ufl.inner(ufl.div(u), q) * ufl.dx
 
-a_blocked = fem.form(ufl.extract_blocks(a))
 
 f = fem.Function(W)
 u_D = fem.Function(V)
@@ -317,43 +316,39 @@ def jump(phi, n):
     + (alpha / h) * ufl.inner(ufl.outer(u_D, n), ufl.outer(v, n)) * ufl.ds
 )
 L += ufl.inner(fem.Constant(msh, default_real_type(0.0)), q) * ufl.dx
-L_blocked = fem.form(ufl.extract_blocks(L))
 
 # Boundary conditions
+msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
 boundary_facets = mesh.exterior_facet_indices(msh.topology)
 boundary_vel_dofs = fem.locate_dofs_topological(V, msh.topology.dim - 1, boundary_facets)
 bc_u = fem.dirichletbc(u_D, boundary_vel_dofs)
 bcs = [bc_u]
 
 
 # Assemble Stokes problem
-A = assemble_matrix(a_blocked, bcs=bcs)
-A.assemble()
-
-b = assemble_vector(L_blocked, kind=PETSc.Vec.Type.MPI)
-bcs1 = fem.bcs_by_block(fem.extract_function_spaces(a_blocked, 1), bcs)
-apply_lifting(b, a_blocked, bcs=bcs1)
-b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
-bcs0 = fem.bcs_by_block(fem.extract_function_spaces(L_blocked), bcs)
-set_bc(b, bcs0)
-
-# Create and configure solver
-ksp = PETSc.KSP().create(msh.comm)
-ksp.setOperators(A)
-ksp.setType("preonly")
-ksp.getPC().setType("lu")
-ksp.getPC().setFactorSolverType("mumps")
-opts = PETSc.Options()  # type: ignore
-opts["mat_mumps_icntl_14"] = 80  # Increase MUMPS working memory
-opts["mat_mumps_icntl_24"] = 1  # Option to support solving a singular matrix (pressure nullspace)
-opts["mat_mumps_icntl_25"] = 0  # Option to support solving a singular matrix (pressure nullspace)
-opts["ksp_error_if_not_converged"] = 1
-ksp.setFromOptions()
-
-# Solve Stokes for initial condition
-x = A.createVecRight()
+solver_options = {
+    "ksp_type": "preonly",
+    "pc_type": "lu",
+    "pc_factor_mat_solver_type": "mumps",
+    "mat_mumps_icntl_14": 80,  # Increase MUMPS working memory
+    "mat_mumps_icntl_24": 1,  # Option to support solving a singular matrix (pressure nullspace)
+    "mat_mumps_icntl_25": 0,  # Option to support solving a singular matrix (pressure nullspace)
+    "ksp_error_if_not_converged": 1,
+}
+u_h = fem.Function(V)
+p_h = fem.Function(Q)
+p_h.name = "p"
+stokes_problem = LinearProblem(
+    ufl.extract_blocks(a),
+    ufl.extract_blocks(L),
+    u=[u_h, p_h],
+    bcs=bcs,
+    kind="mpi",
+    petsc_options=solver_options,
+)
+
 try:
-    ksp.solve(b, x)
+    stokes_problem.solve()
 except PETSc.Error as e:  # type: ignore
     if e.ierr == 92:
         print("The required PETSc solver/preconditioner is not available. Exiting.")
@@ -362,15 +357,6 @@ def jump(phi, n):
     else:
         raise e
 
-# Split the solution
-u_h = fem.Function(V)
-p_h = fem.Function(Q)
-p_h.name = "p"
-offset = V.dofmap.index_map.size_local * V.dofmap.index_map_bs
-u_h.x.array[:offset] = x.array_r[:offset]
-u_h.x.scatter_forward()
-p_h.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
-p_h.x.scatter_forward()
 # Subtract the average of the pressure since it is only determined up to
 # a constant
 p_h.x.array[:] -= domain_average(msh, p_h)
@@ -406,37 +392,26 @@ def jump(phi, n):
     + ufl.inner((ufl.dot(u_n, n))("-") * u_uw, v("-")) * ufl.dS
     + ufl.inner(ufl.dot(u_n, n) * lmbda * u, v) * ufl.ds
 )
-a_blocked = fem.form(ufl.extract_blocks(a))
 
 L += (
     ufl.inner(u_n / delta_t, v) * ufl.dx
     - ufl.inner(ufl.dot(u_n, n) * (1 - lmbda) * u_D, v) * ufl.ds
 )
-L_blocked = fem.form(ufl.extract_blocks(L))
+
+navier_stokes_problem = LinearProblem(
+    ufl.extract_blocks(a),
+    ufl.extract_blocks(L),
+    u=[u_h, p_h],
+    bcs=bcs,
+    kind="mpi",
+    petsc_options=solver_options,
+)
 
 # Time stepping loop
-bcs1 = fem.bcs_by_block(fem.extract_function_spaces(a_blocked, 1), bcs)
 for n in range(num_time_steps):
     t += delta_t.value
 
-    A.zeroEntries()
-    fem.petsc.assemble_matrix(A, a_blocked, bcs=bcs)  # type: ignore
-    A.assemble()
-
-    with b.localForm() as b_loc:
-        b_loc.set(0)
-    assemble_vector(b, L_blocked)
-    apply_lifting(b, a_blocked, bcs=bcs1)
-    b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
-    set_bc(b, bcs0)
-
-    # Compute solution
-    ksp.solve(b, x)
-
-    u_h.x.array[:offset] = x.array_r[:offset]
-    u_h.x.scatter_forward()
-    p_h.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
-    p_h.x.scatter_forward()
+    navier_stokes_problem.solve()
     p_h.x.array[:] -= domain_average(msh, p_h)
 
     u_vis.interpolate(u_h)
 
@@ -721,6 +721,7 @@ def create_eps_mu(
     },
 )
 Esh = problem.solve()
+assert isinstance(Esh, fem.Function)
 assert problem.solver.getConvergedReason() > 0, "Solver did not converge!"
 # -
 
 
@@ -158,6 +158,7 @@
 # +
 problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
 uh = problem.solve()
+assert isinstance(uh, fem.Function)
 # -
 
 # The solution can be written to a {py:class}`XDMFFile
 
@@ -640,6 +640,7 @@ def curl_2d(f: fem.Function):
 gdim = mesh_data.mesh.geometry.dim
 V_dg = fem.functionspace(mesh_data.mesh, ("Discontinuous Lagrange", degree, (gdim,)))
 Esh_dg = fem.Function(V_dg)
+assert isinstance(Esh, fem.Function)
 Esh_dg.interpolate(Esh)
 
 with io.VTXWriter(mesh_data.mesh.comm, "Esh.bp", Esh_dg) as vtx:
Original file line number	Diff line number	Diff line change
`@@ -676,6 +676,7 @@ def create_eps_mu(pml, rho, eps_bkg, mu_bkg):`
`676`	`676`	`},`
`677`	`677`	`)`
`678`	`678`	`Esh_m = problem.solve()`
	`679`	`+ assert isinstance(Esh_m, fem.Function)`
`679`	`680`	`assert problem.solver.getConvergedReason() > 0, "Solver did not converge!"`
`680`	`681`
`681`	`682`	`# Scattered magnetic field`
Original file line number	Diff line number	Diff line change
`@@ -721,6 +721,7 @@ def create_eps_mu(`
`721`	`721`	`},`
`722`	`722`	`)`
`723`	`723`	`Esh = problem.solve()`
	`724`	`+assert isinstance(Esh, fem.Function)`
`724`	`725`	`assert problem.solver.getConvergedReason() > 0, "Solver did not converge!"`
`725`	`726`	`# -`
`726`	`727`