Basic queries from the fenics book, chapter 1, Set #2

Uncategorized
1

1.) How do I save the degrees of unknowns in a file or display in a command prompt in order to have a look?
The following does not work, I receive runtime error

1a.)
u_nodal_values = u.vector()
u_array = u_nodal_values.get_local()
print(u_array)

1b.)

u_nodal_values = u.vector()
u_array = u_nodal_values.get_local()
File('poisson/u_array') << u_array

2.) Briefly, how to relate the u_array to the global mesh coordinates number? Basically I would like to learn the one-to-one correspondence of global index i for an array u.vector().get_local() with the coordiantes returned by the command mesh. coordiantes()

u_nodal_values = u.vector()
u_array = u_nodal_values.get_local()

coor = mesh.coordinates()

print(coor)
print(u_array)

3.) This is supposed to diplay a

str(mesh)

a short “pretty print” description of the mesh. But str(mesh) does not work. What is the updated function name? The following does not work.
print(mesh)
print(str(mesh))

p. 13 of the fenics book.

4.) How lok look for documentation in a terminal? The following does not work.

pydoc ddolfin.Mesh

5.) How the degree can be 1 in the Expression object with exponential function?

page.32, fenics tutorial

p = Expression(’4exp(-pow(beta, 2)(pow(x[0], 2) + pow(x[1] - R0, 2)))’,
degree=1, beta=beta, R0=R0)

6.) Why there is an error in print statement, the syntax seems to be correct.

print(‘sigma = %8g: max deviation = %8g’, %(sigma, dev))

7.) How to plot multiple figures like in matlab, figure(1),…figure(n). When I run in a linux terminal, I have to close one plot to see the next one

plot(f1)
plt.show()

plot(f2)
plt.show()

plot(f3)
plt.show()

  1. Where to look for the elevate attribute? I get error like this.
    AttributeError: ‘TriContourSet’ object has no attribute ‘elevate’

import time
viz_w = plot(w,
wireframe = False,
title = “Scaled membrane deflection”,
rescale = False,
axes = True,
)

viz_w.elevate(-65) # tilt camera -65 degrees (latitude dir)
viz_w.set_min_max(0, 0.5*max_w)
viz_w.update(w) # bring settings above into action
viz_w.write_png(‘membrane_deflection.png’)
viz_w.write_ps(‘membrane_deflection’, format = ‘eps’)

9.) Where to look for the functionality of u.rename?

u.rename(‘u’, ‘solution field’)

This appears in program d5_p2D.py
line 36

It does not substitue the title value when I plot it

plot(u, title = u.name())

10.) I suppose the following is intended to show the summary of variables in a terminal, but it does not show

Quick summary print of key variables

print mesh
print u

2

Can you supply a minimal example ? This should work

import numpy as np
from dolfin import *
V = FunctionSpace(UnitSquareMesh(10, 10), "CG", 1) 
u = Function(V)
u.vector()[:] = np.arange(V.dim())
print(u.vector().get_local())

It’s a numpy array. You can simply use numpy.savetxt() or for larger arrays numpy.save() or pickle or even h5py.

See this

from dolfin import *
help(Mesh)

I may have to check, but you can try and increase the degree to see what happens :slight_smile:

print(r"sigma = %8g: max deviation = %8g"%(sigma, dev)) #this is a python question

You can look at examples in matplotlib. Basically create subplots, and then traverse the axes list. See this for a minimal example

Again, I would suggest create a matplotlib.pyplot figure/axis and then play around with the view_angles

See this. Although it’s old, but still should give you what you want.

I know google results for functions will more often than not lead you to older releases of dolfin but searching the documentation and this forum should lead you to answers with most of the basic questions.

3

Try:

from dolfin import *

mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh, 'CG', 1)

u = Function(V)
print(u.vector()[:])
with XDMFFile('u.xdmf') as file:
    file.write(u)

Try:

for v in vertices(mesh):
    print([v.x(0), v.x(1)], V.dofmap().entity_dofs(mesh, 0, [v.index()]))

This is definetly suboptimal, maybe look into dof_to_vertex_map(). Anyway, both only work for FunctionSpace with all dofs on the vertices of the mesh, e.g. CG1.

Maybe try plotting the mesh. try vtkplotter:

from dolfin import *
from vtkplotter.dolfin import plot

mesh = UnitSquareMesh(5,5)
plot(mesh)

Expressions are interpolated onto finite element function spaces of the specified degree (CG) or a specified element type. If you need higher accuracy for quadrature try defining a UFL-expression like this:

x = SpatialCoordinate(mesh)
p = 4 * exp(-beta**2 * x[0]**2 + (x[1] - R0)**2)

About all plotting issues I recommend looking into vtkplotter / vedo. The original plot commands from earlier dolfin versions are quite outdated and mostly do not work anymore.