# Assembling global LHS from smaller matrices

I want to assemble a global LHS and RHS from multiple equations, but it gives error.
The code is attached:

``````a_01 = -inner(p,div(v))*dx
a_10 = -inner(div(u),q)*dx
L_0 = inner(f,v)*dx
L_1 = inner(Constant(1.0), q)*dx

def jump(phi, n):
return outer(phi("+"), n("+")) + outer(phi("-"), n("-"))

lmbda = conditional(gt(dot(u_n, n), 0), 1, 0)
u_uw = lmbda("+") * u("+") + lmbda("-") * u("-")
a_00 =(inner(u / k, v) * dx
- inner(u, div(outer(v, u_n))) * dx
+ inner((dot(u_n, n))("+") * u_uw, v("+")) * dS
+ inner((dot(u_n, n))("-") * u_uw, v("-")) * dS
+ inner(dot(u_n, n) * lmbda * u, v) * ds)
a_00 = assemble(a_00)
a_01 = assemble(a_01)
a_10 = assemble(a_10)
none_matrix = np.zeros_like(a_00) * np.nan
a = [[a_00,a_01], [a_10,none_matrix]]
a = assemble(a)
L_0 += inner(u_n/k , v) * dx - inner(dot(u_n,n)*(1-lmbda)*u_0,v) * ds
L_0 = assemble(L_0)
L_1 = assemble(L_1)
L = ([L_0,L_1])
L = assemble(L)
``````

How can I resolve this problem?

This is not a reproducible example. It is also unclear if you are using legacy dolfin or dolfinx.

Sure. Here is the complete code. By global LHS I mean the a matrix which is constructed from a_00, a_01 and a_10 in the code. the none is supposed to be zero matrix.
and by global RHS I mean the L vector.

``````from __future__ import print_function
from fenics import *
import dolfin
import numpy as np
import matplotlib
from matplotlib import pyplot
import matplotlib.pyplot as plt

T = 0.1           # final time
num_steps = 50    # number of time steps
dt = T / num_steps # time step size
mu = 1.0             # kinematic viscosity
rho = 1            # density

# Create mesh and define function spaces
P = 100
mesh = UnitSquareMesh(P, P)
u_0 = Expression(('-cos(2*pi*x[0])*sin(2*pi*x[1])', 'cos(2*pi*x[1])*sin(2*pi*x[0])'), degree = 1)
#u_0 = Expression(('sin(2*pi*x[1])', 'cos(2*pi*x[0])'), degree = 1)
class PeriodicBoundary(SubDomain):

# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return (near(x[0], 0) or near(x[1],0)) and on_boundary

# Map right boundary (H) to left boundary (G)
def map(self, x, y):
if near(x[0], 1):
y[0] = x[0] - 1
y[1] = x[1]
elif near(x[1], 1):
y[0] = x[0]
y[1] = x[1] - 1
else:
y[0] = x[0]
y[1] = x[1]
# Define boundary condition
pbc = PeriodicBoundary()
V = FunctionSpace(mesh, 'RT', 2, constrained_domain=pbc)
Q = FunctionSpace(mesh, 'DG', 1,constrained_domain=pbc)

# # Define boundaries
# inflow  = 'near(x[0], 0)'
# outflow = 'near(x[0], 1)'
# walls   = 'near(x[1], 0) || near(x[1], 1)'
pressure_PC = '(near(x[0],0) and near(x[1],0)) || (near(x[0],0) and near(x[1],1)) || (near(x[0],1) and near(x[1],0)) || (near(x[0],1) and near(x[1],1))'
# # Define boundary conditions
# bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
# bcp_inflow  = DirichletBC(Q, Constant(8), inflow)
# bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcp = DirichletBC(Q, Constant(0), pressure_PC)
# bcu = [bcu_noslip]
bcp = [bcp]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_nn = Function(V)
u_nn = project(u_0, V)
u_n = project(u_0, V)
u_  = Function(V)
p_n = Function(Q)
p_  = Function(Q)

# Define expressions used in variational forms
n   = FacetNormal(mesh)
f   = Constant((0, 0))
k   = Constant(dt)
mu  = Constant(mu)
rho = Constant(rho)

# Define variational problem for step 2
####################
a_01 = -inner(p,div(v))*dx
a_10 = -inner(div(u),q)*dx
L_0 = inner(f,v)*dx
L_1 = inner(Constant(1.0), q)*dx

def jump(phi, n):
return outer(phi("+"), n("+")) + outer(phi("-"), n("-"))

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
t += dt
lmbda = conditional(gt(dot(u_n, n), 0), 1, 0)
u_uw = lmbda("+") * u("+") + lmbda("-") * u("-")
a_00 =(inner(u / k, v) * dx
- inner(u, div(outer(v, u_n))) * dx
+ inner((dot(u_n, n))("+") * u_uw, v("+")) * dS
+ inner((dot(u_n, n))("-") * u_uw, v("-")) * dS
+ inner(dot(u_n, n) * lmbda * u, v) * ds)
a_00 = assemble(a_00)
a_01 = assemble(a_01)
a_10 = assemble(a_10)
none_matrix = np.zeros_like(a_00) * np.nan
a = [[a_00,a_01], [a_10,none_matrix]]
a = assemble(a)
L_0 += inner(u_n/k , v) * dx - inner(dot(u_n,n)*(1-lmbda)*u_0,v) * ds
L_0 = assemble(L_0)
L_1 = assemble(L_1)
L = ([L_0,L_1])
L = assemble(L)
solve(a, u_.vector(), L)

# Plot solution
plot(u_)

print('max u:', u_.vector().max())
``````

I think you would have to use cbc block if you are using legacy fenics:

or maybe you could use blockvector and blockmatrix:
https://bitbucket.org/fenics-project/dolfin/src/1c52e837eb54cc34627f90bde254be4aa8a2ae17/python/demo/undocumented/block-matrix/demo_block-matrix.py?at=master#demo_block-matrix.py-40
see: Bitbucket

I added these lines(shown by arrow) in the beginning of the code:

``````from __future__ import print_function
from fenics import *
import dolfin
------- >import ufl_legacy as ufl
------- >from block import *
import numpy as np
import matplotlib
from matplotlib import pyplot
import matplotlib.pyplot as plt
``````

but it gives error:

``````Traceback (most recent call last):
File "/content/NSperiodic-KimMoin-modified.py", line 14, in <module>
from block import *
File "/usr/local/lib/python3.10/dist-packages/block/__init__.py", line 19, in <module>
from .block_bc import block_bc
File "/usr/local/lib/python3.10/dist-packages/block/block_bc.py", line 6, in <module>
from .splitting import split_bcs
File "/usr/local/lib/python3.10/dist-packages/block/splitting.py", line 3, in <module>
from ufl.corealg.map_dag import MultiFunction
File "/usr/local/lib/python3.10/dist-packages/ufl/__init__.py", line 1, in <module>
[20]``````

Seems like cbc block has not been updated since the change from ufl to ufl_legacy.

It would mean that you either have to roll back your dolfin installation, or fix the import statement in cbc block, as shown in

I don’t know how to fix the import statement. it says import ufl_legacy as ufl. I did so but it still gives error.
I also checked this demo of block matrices:

https://bitbucket.org/fenics-project/dolfin/src/946dbd3e268dc20c64778eb5b734941ca5c343e5/python/demo/undocumented/block-matrix/demo_block-matrix.py?at=master#demo_block-matrix.py-1,39:40,46:47,51,54

but the problem is, The final matrix is not an acceptable type by PETSc solvers.
what should I do?

The easiest would be to move to dolfinx, as it has support for block and nest assembly with petsc.
See
https://docs.fenicsproject.org/dolfinx/v0.7.3/python/demos/demo_stokes.html?highlight=nest#nested-matrix-solver

I’m trying to solve a pde on a bi-periodic domain and since It is not possible to use mpc along with DG elements since dS measure is not implemented in mpc yet, The only way I could solve this problem is by fenics.
So converting a blockmatrix type to another type which is recognized by PETSc solver would be much easier.
I would appreciate it if you could somehow help me with it

You need to change the cbc.block source code for it to work, or downgrade your version of dolfin to something prior to the ufl legacy release.