# Linear Equation from FFC generated from variational formulation **URL:** \ **Category:** variational formulation\ **Created:** [March 14, 2023, 12:31am UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675 "2023-03-14T00:31:41Z")\ **Posts on this page:** 14\ **Page:** 1
**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\ **Post date:** [March 14, 2023, 12:31am UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/1 "2023-03-14T00:31:41Z")
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](https://comet-fenics.readthedocs.io/en/latest/demo/periodic_homog_elas/periodic_homog_elas.html) ```auto 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. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 14, 2023, 5:12am UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/2 "2023-03-14T05:12:37Z")
Try ```python A = assemble(a) ``` ---
**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\ **Post date:** [March 14, 2023, 3:08pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/3 "2023-03-14T15:08:59Z")
Thanks @dokken for your response. I am getting as object type A, and I am not able to use or extract any information from A. ```auto A=assemble(a) print(A) ``` Also, for general view, when we got this hexadecimal location. How can we get the object data/ value? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 14, 2023, 3:28pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/4 "2023-03-14T15:28:59Z")
Following [Transform dolfinx assemble matrix to numpy array - #11 by Samuel\_Groth](https://fenicsproject.discourse.group/t/transform-dolfinx-assemble-matrix-to-numpy-array/2589/11), but adapted to legacy DOLFIN, I would do something along the lines of ```python from scipy.sparse import csr_matrix ai, aj, av = as_backend_type(A).mat().getValuesCSR() Asp = csr_matrix((av, aj, ai)) ``` to get a SCIPY CSR matrix. ---
**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\ **Post date:** [March 14, 2023, 4:58pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/5 "2023-03-14T16:58:19Z")
Thanks @dokken for your response. ```auto A=assemble(a) from scipy.sparse import csr_matrix ai, aj, av = as_backend_type(A).mat().getValuesCSR() Asp = csr_matrix((av, aj, ai)) Asp <19851x19851 sparse matrix of type ' with 535905 stored elements in Compressed Sparse Row format> ``` Can you tell me little more about ai,aj,av. I read about `csr_matrix(data,(row,column)) `and I was trying to relate with ai, aj and av based on their size. ```auto len(ai) 19852 len(aj) 535905 len(av) 535905 ``` ```auto Asp = csr_matrix((av, aj, ai)).toarray() # giving Asp as matrix of size 19851 X 19851 ``` But, I am getting difficulty of what they try to interpret. I tried with mesh data. ```auto mesh.num_vertices() 6771 mesh.num_cells() 13232 ``` Can you tell me how Asp is interpreted with respect to form a (bilinear variation form used in Assemble(a))? Also, in as\_backend\_type, Does Asp gives coefficent matrix (K) from Kx=f (say in general with x as fluctuation function to be evaluated) ? and can we evaluate the vector f, since, as\_backend\_type also return a vector. Thanks for your help. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 14, 2023, 5:43pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/6 "2023-03-14T17:43:46Z")
> [@bagla0](#): > > Can you tell me little more about ai,aj,av. I read about `csr_matrix(data,(row,column)) ` Please read the scipy docs: [scipy.sparse.csr\_matrix — SciPy v1.11.4 Manual](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html) > [@bagla0](#): > > ```auto > mesh.num_vertices() > 6771 > mesh.num_cells() > 13232 > > ``` You are using a second order Lagrange space, with the addition of a real space dof. You cannot relate those to mesh entities. > [@bagla0](#): > > Can you tell me how Asp is interpreted with respect to form a (bilinear variation form used in Assemble(a))? > Also, in as\_backend\_type, Does Asp gives coefficent matrix (K) from Kx=f (say in general with x as fluctuation function to be evaluated) ? > and can we evaluate the vector f, since, as\_backend\_type also return a vector. > > Thanks for your help. ```auto b =assemble(L).get_local() x = w.vector().get_local() ``` aligns as Ax=b. Note that his would only work in serial, as numpy/scipy based matrices do not have parallel support (as in being distributed). ---
**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\ **Post date:** [March 15, 2023, 12:05am UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/7 "2023-03-15T00:05:53Z")
Thanks @dokken for your response. It helped me to identify the linear system. The Asp (coefficient matrix) formed is of very large size even with linear elements. ```auto A=assemble(a) from scipy.sparse import csr_matrix ai, aj, av = as_backend_type(A).mat().getValuesCSR() Asp = csr_matrix((av, aj, ai)).toarray() # Asp: 19851x19851 sparse matrix of type ' #with 535905 stored elements in Compressed Sparse Row format> ``` Is there any suitable method/ function in FEniCS to invert Asp (coefficient) matrix efficiently and accurately ? (Mainly the numerical method) I tried the following: ```auto from scipy.sparse import csc_matrix from scipy.sparse.linalg import inv A = csc_matrix(Asp) inv(A) ``` But, it crashes the kernel ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 15, 2023, 5:48am UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/8 "2023-03-15T05:48:03Z")
The whole point of using dolfin is that it invert this matrix (using petsc), and gives you a solution to the problem. You will never want to store the actual inverse, as this is a dense system if you use a direct solver (LUSolver), one can store the LU factorization, which makes the solve step quick. See for instance: [https://fenicsproject.org/qa/6045/repeated-use-of-matrix-with-lu-decomposition/](https://fenicsproject.org/qa/6045/repeated-use-of-matrix-with-lu-decomposition/) [https://fenicsproject.org/pub/tutorial/html/.\_ftut1017.html](https://fenicsproject.org/pub/tutorial/html/._ftut1017.html) ---
**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\ **Post date:** [March 15, 2023, 4:51pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/9 "2023-03-15T16:51:47Z")
Thanks @dokken for your response. I have few doubts: 1. The linear system which FFC gives is of form Aw=b, where FFC performs some calculation to get A^{-1} to solve for w=A^{-1}b. Is there any way to work on FFC code / PETSc solver to get that A^{-1} ? 2. Can I know the intermediate linear form which comes just after minimization of variational statement (a), which later FFC convert to Aw=b form. For example : \int \:\sigma(w)\:\ grad(v)\:dvol \: where w : trial function (fluctuation function), v: test function, dvol: volume gives the minimization output of U\:w\:=\:P\:\epsilon ....(1) analytically, where epsilon is given macro-strain [6X1] vector given as Eps. [in periodic homogenization code (FEniCs example)]. However, the output given by FFC through sparse csr is A w =b; instead of initial U\:w\:=\:P\:\epsilon , where U= A and P \epsilon=b How can I get the intermediate term P ? Can we know the form generated by FFC eqn …(1) instead of direct A w =b ? Any help is appreciated. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 15, 2023, 5:23pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/10 "2023-03-15T17:23:50Z")
> [@bagla0](#): > > The linear system which FFC gives is of form Aw=bAw=bAw=b, where FFC performs some calculation to get A^{-1}A−1A^{-1} to solve for > w=A^{-1}bw=A−1bw=A^{-1}b. > Is there any way to work on FFC code / PETSc solver to get that A^{-1}A−1A^{-1} ? Lets make some things clear here: 1. FFC does not create the in inverse. FFC takes a ufl form, and generates code for integrating over a cell/facet to return the local tensor (scalar, vector, matrix) for a single cell/facet. 2. A= assemble(a) generates the sparse matrix by calling the kernels from step 1 for all cells/facets specified in the form. 3. this matrix is sent to a linear algebra backend, such as PETSc, that uses a direct or iterative solver to get the solution to the problem Ax=b. At no point is the inverse matrix actually formed. See: [KSP: Linear System Solvers — PETSc 3.20.1 documentation](https://petsc.org/release/docs/manual/ksp/) For instance it states > To solve a linear system with a direct solver (currently supported by PETSc for sequential matrices, and by several external solvers through PETSc interfaces, see [Using External Linear Solvers](https://petsc.org/release/docs/manual/ksp/#sec-externalsol)) one may use the options `-ksp_type` `preonly` (or the equivalent `-ksp_type` `none`) `-pc_type` `lu` (see below). For instance you can use the LU factorizer UMFPACK ([MATSOLVERUMFPACK — PETSc 3.20.1 documentation](https://petsc.org/release/docs/manualpages/Mat/MATSOLVERUMFPACK/)) or mumps ([MATSOLVERMUMPS — PETSc 3.20.1 documentation](https://petsc.org/release/docs/manualpages/Mat/MATSOLVERMUMPS/)) An LU factorization never forms the full inverse, see: [LU decomposition - Wikipedia](https://en.m.wikipedia.org/wiki/LU_decomposition) ---
**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\ **Post date:** [March 15, 2023, 9:20pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/11 "2023-03-15T21:20:02Z")
Thanks @dokken for your response. I am facing problem with integrating float type tensor over all the cells. ```auto Em = 50e3 num = 0.2 Er = 210e3 nur = 0.3 material_parameters = [(Em, num), (Er, nur)] nphases = len(material_parameters) ``` ```auto def Com(i): E, nu = material_parameters[i] lmbda = E*nu/(1+nu)/(1-2*nu) mu = E/2/(1+nu) C1=lmbda+2*mu return Constant([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)]) ``` ```auto D= np.zeros((6,6)) for j in range(6): for k in range(6): D[j,k]= assemble(sum([(Com(i)[j,k])*dx(i) for i in range(nphases)]))/vol This integral is missing an integration domain. --------------------------------------------------------------------------- UFLException Traceback (most recent call last) /tmp/ipykernel_455/2734975225.py in 2 for j in range(6): 3 for k in range(6): ----> 4 D[j,k]= assemble(sum([(Com(i)[j,k])*dx(i) for i in range(nphases)]))/vol /tmp/ipykernel_455/2734975225.py in (.0) 2 for j in range(6): 3 for k in range(6): ----> 4 D[j,k]= assemble(sum([(Com(i)[j,k])*dx(i) for i in range(nphases)]))/vol /usr/lib/python3/dist-packages/ufl/measure.py in __rmul__ (self, integrand) 430 domain, = domains 431 elif len(domains) == 0: --> 432 error("This integral is missing an integration domain.") 433 else: 434 error("Multiple domains found, making the choice of integration domain ambiguous.") /usr/lib/python3/dist-packages/ufl/log.py in error(self, *message) 156 "Write error message and raise an exception." 157 self._log.error(*message) --> 158 raise self._exception_type(self._format_raw(*message)) 159 160 def begin(self, *message): UFLException: This integral is missing an integration domain. ``` Can you help me to compute the assemble command ? I am getting error with multiple attempts like changing Com function `return as_tensor([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)]) ` or `return ([(C1,lmbda,lmbda,0,0,0),(lmbda,C1,lmbda,0,0,0),(lmbda,lmbda,C1,0,0,0),(0,0,0,mu,0,0),(0,0,0,0,mu,0),(0,0,0,0,0,mu)]) ` Any help is greatly appreciated. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 15, 2023, 9:22pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/12 "2023-03-15T21:22:33Z")
The code you have posted is not complete, as you have not defined `dx` and it is unclear what `dx(i)` is supposed to represent. Please see:[Read before posting: How do I get my question answered? - #4 by nate](https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21/4) for examples following the forum guidelines ---
**Author:** ![bagla0](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/bagla0/32/4677_2.png) [@bagla0](https://fenicsproject.discourse.group/u/bagla0)\ **Post date:** [March 15, 2023, 9:59pm UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/13 "2023-03-15T21:59:46Z")
> [@bagla0](#): > > `dx = Measure('dx')(subdomain_data=subdomains)` dx is taken from the above step and subdomain is similar to the periodic homogenization example. ```auto from __future__ import print_function from dolfin import * import numpy as np import matplotlib.pyplot as plt %matplotlib notebook a = 1. # unit cell width b = sqrt(3.)/2. # unit cell height c = 0.5 # horizontal offset of top boundary R = 0.2 # inclusion radius vol = a*b # unit cell volume # we define the unit cell vertices coordinates for later use vertices = np.array([[0, 0.], [a, 0.], [a+c, b], [c, b]]) fname = "hexag_incl" mesh = Mesh(fname + ".xml") subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml") facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml") plt.figure() plot(subdomains) plt.show() ``` I apologize for the incomplete step. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [March 16, 2023, 5:47am UTC](https://fenicsproject.discourse.group/t/linear-equation-from-ffc-generated-from-variational-formulation/10675/14 "2023-03-16T05:47:40Z")
Then change > [@bagla0](#): > > `dx = Measure('dx')(subdomain_data=subdomains)` To `dx = Measure('dx')(domain=mesh, subdomain_data=subdomains)`