### Print or save matrices from functional or weak form

250
views
1
7 months ago by
Dear community,

I apologize for being new to fenics. I would like to print or, better, save the matrix of my system to read it in Matlab. In particular I would like to do this:

1) Starting from a functional J, obtain its Jacobian. If I write

W = FunctionSpace(mesh, feW)
w = Function(W)
J= functional
G = derivative(J,w)
bcs= boundary conditions
solve(G==0,w,bcs)

then I can solve the problem, but by doing:

Gmat = assemble(G)
print(Gmat.array() )

G seems only to be a vector (maybe of the residuals?). How can I obtain the matrix of the system/jacobian of the functional?

2) I would like to do the same with variational formulation, expressed with test functions.

3) How can I do 1) and 2) with/without strongly imposed BCS?

Thank you very much in advance. Any help would be really appreciated.

Community: FEniCS Project
Can you explicitly give your functional? Or better: Post a MWE.

If you have a bilinear Form $J\left(u,v\right)$J(u,v) then its derivative at a known point  $u_0$u0  will be a bilinear form with the direction of the derivative as its first argument:
$\mathrm{d}J=\frac{\mathrm{d}}{\mathrm{d\varepsilon}}J\left(u_0+\varepsilon\mathrm{d}u,v\right)_{\varepsilon=0}=G\left(\mathrm{d}u,v\right)$dJ=ddε J(u0+εdu,v)ε=0=G(du,v)
And after assembly it should give a matrix.
written 7 months ago by klunkean
Thank you Klunkean for your reply. In the following code, I would like to obtain the matrix related to the functional J, with and without bc, if this is possible.
from dolfin import *
from scipy.io import savemat

aL = 0.0
aR = 1.0
bB = 0.0
bT = 0.2
Na = 100
Nb = 80
C_equilibrium=1
C_constitutive=1.0
C_boundary=100.
mesh = RectangleMesh(Point(aL, bB),Point(aR, bT), Na,Nb) #,"crossed"

# some boundaries definition
def boundary_L(x):
return x[0] < aL + DOLFIN_EPS
def boundary_B(x):
return x[1] < bB + DOLFIN_EPS
def boundary_R(x):
return x[0] > aR - DOLFIN_EPS
def boundary_T(x):
return x[1] > bB - DOLFIN_EPS

# parameters
Lambda = Constant(1.0)
Mu = Constant(1.0)
beta=Constant(1.0/(2*Mu))
alpha = Constant(-beta * Lambda /(2*Lambda +2*Mu))
gamma=Constant(Lambda+2*Mu);

# boundary or volumetric forces
f = Expression(("0.0","0.0"),degree=4)
g = Expression(("0.0","0.0"),degree=4)
h = Expression((("0.0","0.0"),("0.0","0.00")),degree=2)
j = Expression((("0.0","-0.005"),("0.0","0.0")),degree=2)
kk = Expression(("0.0","-0.005"),degree=2)

# operator definition
def A(sigma):
return beta*sigma+alpha*tr(sigma)*Identity(2)
def epsilon(u):

# functional space definition
fe2DP1 = VectorElement("P",mesh.ufl_cell(),1)
fe2DRT = VectorElement("RT",mesh.ufl_cell(),1)
feW  = MixedElement([fe2DRT,fe2DP1])

W = FunctionSpace(mesh, feW)
w = Function(W)
sigma,u = split(w)

# functional definition
J =C_constitutive * inner(A(sigma)-epsilon(u),A(sigma)-epsilon(u))*dx + C_equilibrium * inner(div(sigma)+f,div(sigma)+f)*dx
G = derivative(J,w)

# bc definition
bc1 = DirichletBC(W.sub(1), g, boundary_L)
bc2 = DirichletBC(W.sub(0), h, boundary_B)
bc4 = DirichletBC(W.sub(0).sub(1), kk, boundary_T)
bc3 = DirichletBC(W.sub(0), h, boundary_R)
bcs=[bc1,bc2,bc4,bc3]

# solve the problem
solve(G==0,w,bcs)
sigma,u = w.split()

# print the matrix
Gmat = assemble(G)
print(Gmat.array() )

# save the output
ufile_pvd = File("lsfem_elasticity_u.pvd")
ufile_pvd << u
ufile_pvd = File("lsfem_elasticity_sigma.pvd")
ufile_pvd << sigma
​

written 7 months ago by Garo
Your functional J does not contain any test functions. This does not seem to me as a valid weak formulation for some problem. Since J is not a bilinear form, it won't yield a matrix after assembly and neither won't its derivative.
Maybe rethink your J. Where does this formulation come from? I can identify the balance of linear momentum and some term to weakly set the strains to some value but can't figure out how you got to this functional...
written 7 months ago by klunkean
Dear Klunkean, I consider the functional as the energy-functional. The bilinear form is its derivative. I wanted to write the functional J and then obtain the bilinear form G as its derivative to reduce the possible error in writing the bilinear form directly. So in my case I would say that G is the bilinear form. Anyhow I have also the code written in terms of bilinear form (posted below).Can you please tell me how can I obtain the matrix of my linear system, with and without bcs? thank you.
from dolfin import *
from scipy.io import savemat

aL = 0.0
aR = 1.0
bB = 0.0
bT = 0.2
Na = 80
Nb = 50
C_equilibrium=1
C_constitutive=1.0
C_boundary=100.
mesh = RectangleMesh(Point(aL, bB),Point(aR, bT), Na,Nb) #,"crossed"

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary_L(x):
return x[0] < aL + DOLFIN_EPS
def boundary_B(x):
return x[1] < bB + DOLFIN_EPS
def boundary_R(x):
return x[0] > aR - DOLFIN_EPS
def boundary_T(x):
return x[1] > bB - DOLFIN_EPS

Lambda = Constant(1.0)
Mu = Constant(1.0)
beta=Constant(1.0/(2*Mu))
alpha = Constant(-beta * Lambda /(2*Lambda +2*Mu))
gamma=Constant(Lambda+2*Mu);

fx= 0.0
fy= -0.005

f = Expression(("0.0","0.0"),degree=4)
g = Expression(("0.0","0.0"),degree=4)
h = Expression((("0.0","0.0"),("0.0","0.00")),degree=2)
bcforce = Expression(("0.0","-0.005"),degree=2)
j = Expression((("0.0","-0.005"),("0.0","0.0")),degree=2)
kk = Expression(("0.0","-0.005"),degree=2)

def A(sigma):
return beta*sigma+alpha*tr(sigma)*Identity(2)
def epsilon(u):
fe2DP1 = VectorElement("P",mesh.ufl_cell(),1)
fe2DRT = VectorElement("RT",mesh.ufl_cell(),1)
feW  = MixedElement([fe2DRT,fe2DP1])

W = FunctionSpace(mesh, feW)
w = Function(W)
sigma,u = split(w)

(sigma,u)=TrialFunctions(W)
(tau,v)=TestFunctions(W)

a =   C_equilibrium  * inner(div(sigma),div(tau))*dx
a +=  C_constitutive * inner(A(sigma),A(tau))*dx
a += -C_constitutive * inner(epsilon(u),A(tau))*dx
a += -C_constitutive * inner(A(sigma),epsilon(v))*dx
a +=  C_constitutive * inner(epsilon(u),epsilon(v))*dx

L = -C_equilibrium * inner(f,div(tau))*dx

# bcs
bc1 = DirichletBC(W.sub(1), g, boundary_L)
bc2 = DirichletBC(W.sub(0), h, boundary_B)
bc4 = DirichletBC(W.sub(0).sub(1), kk, boundary_T)
bc3 = DirichletBC(W.sub(0), h, boundary_R)
bcs=[bc1,bc2,bc4,bc3]
A, bb = assemble_system(a, L, bcs)
solver = KrylovSolver(A, "gmres")

# solve for LSFEM APPROACH
solve(a == L, w, bcs)
sigma,u = w.split()
ufile_pvd = File("assemble_lsfem_elasticity_u.pvd")
ufile_pvd << u
ufile_pvd = File("assemble_lsfem_elasticity_sigma.pvd")
ufile_pvd << sigma

​
written 7 months ago by Garo
If I run this code and just
print A.array()​
it returns a matrix. This is because the TestFunctions tau and v are arguments to your bilinear form. Maybe try to incorporate the TestFunctions in your derivative formulation.
written 7 months ago by klunkean