Record the use process of fenics

Occasionally I need to use fenics To calculate , But all kinds of unfamiliar , Every time I look online , There is a lot of information , and fenics The updated version , Many functions report errors , So here is a record of your own use process .

At the same time, it is necessary to refer to fenics Introduction to basic types in : Docs » User manual » Form language).

1. After defining the grid and space , Get the number of degrees of freedom of the space 、 The value of the function at the degree of freedom 、

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(4, 4)

#  Coordinate values at mesh nodes 
mesh_coor = mesh.coordinates()

#  Building scalar and vector spaces 
space_S1 = FunctionSpace(mesh, 'CG', 1)  # 1  Express  Lagrange 1  Dimension 
space_S3 = VectorFunctionSpace(mesh, 'CG', 1, dim=3)

#  Get the number of degrees of freedom 
dof_f1 = np.array(df.dof_to_vertex_map(space_S1), dtype=int)
dof_f3 = np.array(df.dof_to_vertex_map(space_S3), dtype=int)

#  establish  space_S1  Interpolation on  f1,  And get the value of the function at the degree of freedom 
f1 = interpolate(df.Expression("1.0+0*x[1]", degree=1), space_S1)
val_dof_f1 = f1.vector().get_local()  

#  Calculate the function value at the mesh node 
v_vertex = f1.compute_vertex_values(mesh)

#  Take the component of the gradient  
Df0 = grad(f1)[0]
Df1 = grad(f1)[1]
# #  stay  Docs » User manual » Form language  Mid search  grad  You can see :
# # In UFL, the following pairs of declarations are equivalent:
# Dfi = grad(f)[i]
# Dfi = f.dx(i)
# Dvi = grad(v)[i, j]
# Dvi = v[i].dx(j)
# DAi = grad(A)[..., i]
# DAi = A.dx(i)
# # for a scalar expression f, a vector expression v, and a tensor expression A of arbitrary rank.

above And get the value of the function at the degree of freedom It should be noted that fenics It started with f1.vector().array(), But this latest version is no longer available ( See Replace vector().array() to vector().get_local()).

2. take Two functions as components form a vector function

from dolfin import *

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

u = Function(V)
u.interpolate(Expression("1+0*x[0]", degree=1))
v = Function(V)
v.interpolate(Expression("2+0*x[0]", degree=1))
W = VectorFunctionSpace(mesh, 'CG', 1, dim=2)
w = Function(W)
ua = u.vector().get_local()
va = v.vector().get_local()
wa = w.vector().get_local()
wa[::2] = ua
wa[1::2] = va
w.vector().set_local(wa)

3. The resulting stiffness matrix is transformed into numpy.array() Format

a = inner(grad(u), grad(v)) * dx
aa = assemble(a)
aM = aa.array()

L = inner(f, q) * dx
ll = assemble(L)
l1[:]  #  This lists all the values of the right-hand end item 

obtain Constant Value

aa = Constant(1e-3)
print(aa) # Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 49)
aa_ = float(aa)
print(aa_)

Remove from the solver “Solving linear variational problem” Tips

add to set_log_active(False) You can turn off all prompts ( Include ‘Error’, ‘Warning’ …)

Get the value of the interpolation function at a certain point

For example, it defines fenics Space Vh, And given True solution Phi = Expression("x[0] + x[1]", degree=3) ( namely Phi = x + y), then Phi Interpolate to Vh In the space : Phi_h = interpolate(Phi, Vh), We want to be spot P=(0.3, 1.2) The value of the function at

Phi_h = interpolate(Phi, Vh)
P = (0.3, 1.2)
val = Phi_h(P)