### How can I calculate SVD in fenics

58
views
0
8 weeks ago by
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.
""" 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

+ 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