### Print stiffness matrix from bilinear form

155
views
-1
3 months ago by
Hi All,

I am Trying to find stiffness matrix of a bilinear form as follows;

a((Z,Q) , (yup, R)) = <<Q, grad(yup)>> + <<z, yup>> + <<grad(Z) , grad(yup)>> + <<Q , R>> + <<div(Q) , div(R)>>

$\left(Z,yup\right)\in H^1$(Z,yup)H1
$\left(Q,R\right)\in H^d$(Q,R)Hd

Which << , >> refers to Both the L2-inner product of vector fields and the L2-inner product of tensor fields.
Each of Q , Z , yup and R has three components for example: div(Q)
\begin{equation*}
\mathbf{div}\,\boldsymbol{Q}=\mathbf{div}\left[ \begin{array}{c} \boldsymbol{Q}_{1} \\ \boldsymbol{Q}_{2} \\ \boldsymbol{Q}_{3} \end{array}\right] = \left[ \begin{array}{c} \mathrm{div}\,\boldsymbol{Q}_{1} \\ \mathrm{div}\,\boldsymbol{Q}_{2} \\ \mathrm{div}\,\boldsymbol{Q}_{3} \end{array}\right],
\end{equation*}

here you can find script:

from dolfin import *
import time

start_time = time.time()
# Create mesh
mesh = UnitSquareMesh(1, 1)
# Define function spaces
RT = FiniteElement("RT", mesh.ufl_cell(), 1)
CG = FiniteElement("CG", mesh.ufl_cell(), 1)
W = CG * RT
V = FunctionSpace(mesh,W)
#Define DOFs
Gdof = V.dim()
#Ldof = V.dofmap().local_dimension("unowned") # For both owned and unowned dofs
# Define trial and test functions
(Q1 ,Z1) = TrialFunctions(V)
(Q2 ,Z2) = TrialFunctions(V)
(Q3 ,Z3) = TrialFunctions(V)
(yup1 ,R1) = TestFunctions(V)
(yup2 ,R2) = TestFunctions(V)
(yup3 ,R3) = TestFunctions(V)

# Define bilinear form and right hand side
+ Q1*R1 + Q2*R2 + Q3*R3)*dx

#####################################################
matA = assemble(a)
stringA = str(matA.array())
outf = open("stiffmat.txt","w")
outf.write(array2string(matA.array(),max_line_width=sys.maxint))
outf.close()

exit()

print('Global_DOF',Gdof)
print('---%s seconds ---' %(time.time() - start_time))
#Hold plot
#interactive()
​

The error is:
matA = assemble(a)
^
SyntaxError: invalid syntax
Community: FEniCS Project
2
I have resolved the problem:

Here you can find the corrected script.

from dolfin import *
from petsc4py import PETSc
from numpy import array2string
import sys
import time

start_time = time.time()
# Create mesh
mesh = UnitSquareMesh(1, 1)
# Define function spaces
RT = FiniteElement("RT", mesh.ufl_cell(), 1)
CG = FiniteElement("CG", mesh.ufl_cell(), 1)
W = CG * RT
V = FunctionSpace(mesh,W)
w = Function(V)
#Define DOFs
Gdof = V.dim()

# Define trial and test functions
(Z1 ,Q1) = TrialFunctions(V)
(Z2 ,Q2) = TrialFunctions(V)
(Z3 ,Q3) = TrialFunctions(V)
(yup1 ,R1) = TestFunctions(V)
(yup2 ,R2) = TestFunctions(V)
(yup3 ,R3) = TestFunctions(V)

# Define bilinear form and right hand side
+ dot(Q1, R1) + dot(Q2, R2) + dot(Q3, R3) \
+ ((div(Q1) * div(R1))) + ((div(Q2) * div(R2))) + ((div(Q3) * div(R3))))*dx
#####################################################

print('Global_DOF',Gdof)

matA = assemble(a)
stringA = str(matA.array())
outf = open("stiffmat.txt","w")
outf.write(array2string(matA.array(),max_line_width=sys.maxint))
outf.close()
print('---%s seconds ---' %(time.time() - start_time))
exit()​

written 3 months ago by Ali
(Q1 ,Z1) = TrialFunctions(V)
(Q2 ,Z2) = TrialFunctions(V)
(Q3 ,Z3) = TrialFunctions(V)
(yup1 ,R1) = TestFunctions(V)
(yup2 ,R2) = TestFunctions(V)
(yup3 ,R3) = TestFunctions(V)

is probably not doing what you want. Q1 and Q2 are the same. Those are the same trial functions on V.

written 3 months ago by Jan Blechta
Then, how can I resolve that?
written 3 months ago by Ali
I guess you want a mixed space. Refer to class MixedElement.
written 3 months ago by Jan Blechta

0
3 months ago by
Hello Ali Gerami,

Don't you have to integrate the expression over the domain? Like:

matA = assemble(a * dx) ???

I think there is a similar question here:
1
Sorry, Now I saw that you are integrating it.

I did a quick test, and it looks like that you have rank dimension problem:

ufl.log.UFLException: Dot product requires non-scalar arguments, got arguments with ranks 0 and 1.

Taking a small piece of the code, this operation has rank dimension error

a = dot(Q1, grad(yup1)) * dx
written 3 months ago by Ricardo D. Lahuerta
-1
3 months ago by
Based on your error, I create a simple example to test it using linear elasticity, and it works, I think you are not defining your problem properly...

Here is the code sample:
from __future__ import print_function
from fenics import *

import numpy as np
import sys

L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'CG', 2)

# Define boundary condition
tol = 1E-14
d = 3  # space dimension

class ClampBound(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and on_boundary

class RightBound(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], L)

clamped_bc = ClampBound()
right_boundary = RightBound()

boundary_part = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_part.set_all(0)

clamped_bc.mark(boundary_part, 2)
right_boundary.mark(boundary_part, 3)

ds = Measure("ds")[boundary_part]

# plot bcs
ds_mark = File("ds_mark.pvd")
ds_mark << boundary_part

# boundary
bc = DirichletBC(V, Constant((0, 0, 0)), clamped_bc)

def epsilon(u):

# solving the elasticity problem
u = Function(V)
v = TestFunction(V)
t = TrialFunction(V)

T = Constant((0.05, 0., 0.))

Psi = (lambda_/2)*(tr(epsilon(u)))**2+mu*tr(epsilon(u)*epsilon(u))

Pi = Psi*dx - inner(T, u)*ds(3)

dPi = derivative(Pi, u, v)
d2Pi = derivative(dPi, u, t)

matA = assemble(d2Pi)
stringA = str(matA.array())

outf = open("stiffmat.txt", "w")
outf.write(np.array2string(matA.array(), max_line_width=sys.maxint))
outf.close()​