Solution of Mixed formulation

1

Hi,
I’m trying to know the value of the solution sigma. I would like to view its values. Can someone help me please.

from dolfin import *
import matplotlib.pyplot as plt

n = 4
mesh = UnitSquareMesh(n, n)
BDM = FiniteElement(“RT”, mesh.ufl_cell(), 1)
DG = FiniteElement(“DG”, mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, BDM * DG)
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
f = Expression(“10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)”, degree=2)
a = (dot(sigma, tau) + div(tau)*u + div(sigma)v)dx
L = - fvdx

class BoundarySource(UserExpression):
def init(self, mesh, **kwargs):
self.mesh = mesh
super().init(**kwargs)
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = sin(5x[0])
values[0] = gn[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)

G = BoundarySource(mesh, degree=2)

def boundary(x):
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
bc = DirichletBC(W.sub(0), G, boundary)

w = Function(W)
solve(a == L, w, bc)
(sigma, u) = w.split()

print(sigma)
print(u)

2

Please encapsulate your code with ``` to format it properly.
When you say view the solution, what do you imply?

  1. look at the numpy array of dof values
sigma, u= w.split(deepcopy=True)
print(sigma.vector().get_local()
  1. visualize the values on the mesh:
plot(sigma)
plt.show()

You can also use vedo for this.
3. Save it to file and visualize it in Paraview:

File("sigma.pvd") << sigma
1 Like
3

Thank you this was very helpful.
One question, now if my sigma is a matrix and i’m interested to know the trace of the solution how can I do this??

4

Project the trace to a suitable function space and visualize it as above.

5

I just need to check that the values of the trace of the matrix( ie my solution) that it sums to zero or not.

What I have is
if I
print(sigma)
print(s.vector().get_local())

I get the output:

[
[f_33[0], f_33[1]],
[f_33[2], f_33[3]]
]

AttributeError: ‘ListTensor’ object has no attribute ‘vector’

How can I overcome this please to get the trace of above matrix (my solution)??

Thanks

6

Please provide a minimal code example that is properly formatted and reproduces your error. As I stated above, encapsulate your code with ``` to get proper formatting.

7

Can you give me an example on how to project the solution, which is a tensor, ??
I’ll try to encapsulate the code,

8

You still havent encapsulated the code.

An example of this behavior is:

from dolfin import *

mesh = UnitSquareMesh(10,10)
degree = 2
V = VectorFunctionSpace(mesh, "CG", degree)
x = SpatialCoordinate(mesh)
u = project(as_vector((x[0], cos(x[1]))), V)
Q = TensorFunctionSpace(mesh, "DG", degree-1)
q = project(grad(u), Q)

Please familarize yourself with the demos.

9

I don’t know how to encapsulate the code as I’m not the developer.

10

I already tried to explain it to you above. You have to encapsulate the code with ```.
To be even more explicit, see: https://github.com/adam-p/markdown-here/wiki/Markdown-Cheatsheet#code-and-syntax-highlighting

11

I thank you for your prompt response.
Please can you help me in the last part… viewing the sigma at the nodes to check the trace…
I hope I did it right the way you told me, but i’m not sure

from future import print_function
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
import petsc4py as PETSc

n = 2

mesh = UnitSquareMesh(n,n,“right”)

RTel = VectorElement(“RT”, mesh.ufl_cell(), 1)
CGel = VectorElement(“CG”, mesh.ufl_cell(), 1)

element = MixedElement([RTel, CGel])
W = FunctionSpace(mesh,element)

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

bcs = DirichletBC(W.sub(1),Constant((0.0, 0.0)) , DomainBoundary())

SS = pow(2,-1)*(sigma-pow(2,-1)*tr(sigma)Identity(2))
TT = pow(2,-1)(tau-pow(2,-1)*tr(tau)*Identity(2))
a = (inner(SS,TT)+dot(div(sigma),div(tau))-inner(SS,sym(grad(u))))*dx
b = -dot(u,div(tau))*dx
A = PETScMatrix()
B = PETScMatrix()

assemble_system(a, bcs, A_tensor=A)
assemble_system(b, bcs, A_tensor=B)
bcs.zero(B)

solver.solve()

num_of_converged_eig_val = solver.get_number_converged()

computed_eigenvalues = []
for i in range(num_of_converged_eig_val):
r = solver.get_eigenvalue(i)
computed_eigenvalues.append(complex(r[0],r[1]))

computed_eigenvalues_list = np.array(computed_eigenvalues)

_, _, x_real, _=solver.get_eigenpair(0)

v = Function(W,x_real)

sigma, u = split(v)
‘’‘I would like to view the values of the sigma on each node to check its trace
print(sigma)
‘’’
plot(u)
plt.show()

12

You still do not encapsulate your code with ```.
You should not use VectorElement on an RT space, as the dofs are integral moments in normal direction over each facet. I would suggest to FiniteElement for RT. I am not sure what you would like to do.

13

It didn’t work when I switched it.

I’ll try to read what you sent me to encapsulate the code.

Thanks a lot… appreciate your help.

14

Your code should be encapsulated between triple backticks (the key with just on the left of 1 (not the keypad but near the Q ) on the QWERTY keyboard).

You can also see this for markdown usage. Preformatting helps others copy-paste your code and run it to help you debug…

To get started, copy-paste this string as-is,

```
import fenics
import numpy as np
```
15

Thank you, Now I get it,