### How can I calculate SVD in fenics

58

views

0

Hi All,

I am Trying to figure out how can I calculate the smallest singular value in fenics.

I have mentioned an example in below that I was trying to print out the smallest singular value of stiffness matrix A.

But It is not working.

Thanks,

~Ali

I am Trying to figure out how can I calculate the smallest singular value in fenics.

I have mentioned an example in below that I was trying to print out the smallest singular value of stiffness matrix A.

```
""" Solve the Inf-Sup in fenics.
TODO: consider to convert that to 3D form.
"""
import dolfin
from dolfin import *
import scipy
from scipy.linalg import sqrtm
from petsc4py import PETSc
from numpy import array2string
import numpy as np
import sys
import mshr
#Create mesh
mesh = UnitSquareMesh(12, 12)
#Define function spaces
BDM = VectorElement("BDM", mesh.ufl_cell(), 2, dim=2)
CG = VectorElement("CG", mesh.ufl_cell(), 1, dim=2)
element = MixedElement([CG, BDM])
V = FunctionSpace(mesh, element)
w = Function(V)
#Define DOFs
Gdof = V.dim()
# Define trial and test functions
Z ,Q = TrialFunctions(V)
yup ,R = TestFunctions(V)
# Define bilinear form and right hand side
a = (dot(Q[0], grad(yup[0])) + dot(Q[1], grad(yup[1])) \
+ dot(Z[0], yup[0]) + dot(Z[1], yup[1]) + dot(grad(Z[0]), grad(yup[0])) + dot(grad(Z[1]), grad(yup[1])) \
+ dot(Q[0], R[0]) + dot(Q[1], R[1]) + ((div(Q[0]) * div(R[0]))) + ((div(Q[1]) * div(R[1]))))*dx
#####################################################
print('Global_DOF',Gdof)
matA = assemble(a)
#print out stiffness matrix
stringA = str(matA.array())
outf = open("stiffmat.txt","w")
outf.write(array2string(matA.array(),max_line_width=sys.maxint))
outf.close()
exit()
#calculate SVD of stiffness matrix
S = SLEPc.SVD()
S.create()
S.setOperator(matA)
S.setType(S.Type.LANCZOS)
S.solve()
nconv = S.getConverged()
if nconv > 0:
v, u = matA.getVecs()
for i in range(nconv):
sigma = S.getSingularTriplet(i, u, v)
print sigma
```

But It is not working.

Thanks,

~Ali

Community: FEniCS Project

Please login to add an answer/comment or follow this question.