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

250

views

1

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.

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

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):
return sym(grad(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...

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):
return sym(grad(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

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

If you have a bilinear Form $J\left(u,v\right)$

J(u,v) then its derivative at a known point $u_0$u_{0}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(u_{0}+εdu,v)_{ε=0}=G(du,v)And after assembly it should give a matrix.