Hi FEniCS family,

I am working on on the following variational statement taken from

https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html

```
Ve = VectorElement("CG", mesh.ufl_cell(), 2)
Re = VectorElement("R", mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([Ve, Re]), constrained_domain=PeriodicBoundary(vertices, tolerance=1e-10))
V = FunctionSpace(mesh, Ve)
v_,lamb_ = TestFunctions(W)
dv, dlamb = TrialFunctions(W)
w = Function(W)
dx = Measure('dx')(subdomain_data=subdomains)
Eps = Constant(((0, 0), (0, 0)))
F = sum([inner(sigma(dv, i, Eps), eps(v_))*dx(i) for i in range(nphases)])
a, L = lhs(F), rhs(F)
a += dot(lamb_,dv)*dx + dot(dlamb,v_)*dx
```

For the variational statement `a, L = lhs(F), rhs(F)`

, the FEniCS find a linear equation (using FFC ) and after solve `(a==L,w,[])`

give output of fluctuation function, w for each node.

I wanted to find the linear system of equation generated by FEniCS (like Ax=b, with x=w being the fluctuation function). I want to use coefficient matrix A for my post processing.

I know that, FFC forms a C++ code from symbolic variational statement, and after solution gives back, function w in python format. So, there must be some linear equation which is generated to solve the variational form (a==L). Kindly tell me, how I can get access to the linear system of equation (coeff matrix, A being assembled over all elements of mesh, through guassian quadrature) and interpret it for usage.

Any help is greatly appreciated.