### Print stiffness matrix from bilinear form

155

views

-1

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$(

$\left(Q,R\right)\in H^d$(

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:

The error is:

matA = assemble(a)

^

SyntaxError: invalid syntax

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`,`y``u``p`)∈`H`^{1}$\left(Q,R\right)\in H^d$(

`Q`,`R`)∈`H`^{d}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
a = (dot(Q1, grad(yup1)) + dot(Q2, grad(yup2)) + dot(Q3, grad(yup3)) \
+ Z1*yup1 + Z2*yup2 + Z3*yup3 + dot(grad(Z1), grad(yup1)) + dot(grad(Z2), grad(yup2)) + dot(grad(Z3), grad(yup3)) \
+ 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 Answers

0

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:

https://www.allanswered.com/post/erexl/#kgnbz

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

matA = assemble(a * dx) ???

I think there is a similar question here:

https://www.allanswered.com/post/erexl/#kgnbz

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:

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

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:

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):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
# 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()
```

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

Here you can find the corrected script.

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

`MixedElement`

.