Transferring my scripts from dolfin to dolfin-x: ufl.abs() is not accepted by ufl

1

Dear community,
While transferring my scripts from dolfin to dolfin-x, I have get some errors. Here is a (more or less) minimal example of my code.

import dolfinx
import ufl
import numpy as np
from mpi4py import MPI
#create mesh
mesh = dolfinx.RectangleMesh(MPI.COMM_WORLD, np.array([[0,0,0],[20, 200,0]]), [30,30], cell_type=dolfinx.cpp.mesh.CellType.triangle)
# create facet-tags
boundaries = [(1, lambda x: np.isclose(x[0], 0.0)), (2, lambda x: np.isclose(x[1], 200.0)), (3, lambda x: np.isclose(x[1], 0.0))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
   facets = dolfinx.mesh.locate_entities(mesh, fdim, locator)
   facet_indices.append(facets)
   facet_markers.append(np.full(len(facets), marker))
facet_tags = dolfinx.MeshTags(mesh, fdim, np.array(np.hstack(facet_indices), dtype=np.int32),
                            np.array(np.hstack(facet_markers), dtype=np.int32))

# define element and FunctionSpace
e0 = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
V = dolfinx.FunctionSpace(mesh, e0)                 # V = dolfinx.VectorFunctionSpace(mesh, ("CG", 1))
e1 = ufl.FiniteElement("CG", mesh.ufl_cell(), 1)
W = dolfinx.FunctionSpace(mesh,e1)                  #  W = dolfinx.FunctionSpace(mesh,("CG",1))
WW = dolfinx.FunctionSpace(mesh, ('DG', 0))

## Boundary conditions
bc_u =[]
bc_p =[]

fdim = mesh.topology.dim - 1

V0 = V.sub(0).collapse()
u_left = dolfinx.Function(V0)
with u_left.vector.localForm() as u_left_loc:
   u_left_loc.set(0)
dofs_left = dolfinx.fem.locate_dofs_topological((V.sub(0), V0), fdim, facet_tags.indices[facet_tags.values == 1]) # x=0 left_tag=1
bc_left = dolfinx.DirichletBC(u_left, dofs_left, V.sub(0))
bc_u.append(bc_left)  # add bondary condition bc_left

V1 = V.sub(1).collapse()
u_bottom = dolfinx.Function(V1)
with u_bottom.vector.localForm() as u_bottom_loc:
   u_bottom_loc.set(0)

dofs_bottom = dolfinx.fem.locate_dofs_topological((V.sub(1), V1), fdim, facet_tags.indices[facet_tags.values ==3]) # y=0 bottom_tag=3
bc_bottom = dolfinx.DirichletBC(u_bottom, dofs_bottom, V.sub(1))
bc_u.append(bc_bottom)  # add bondary condition bc_left

class MyExpression:
   def __init__(self):
       self.t = 0.0
   def eval(self, x):
       return np.full(x.shape[1], self.t)    # Added some spatial variation here. Expression is sin(t)*x

u_ex = MyExpression()    
u_top = dolfinx.Function(V1)
u_top.interpolate(u_ex.eval) # interpolate expression into function 
u_top.x.scatter_forward()    # sending (scattering) the ghost values in the underlying data structure of u_top   
dofs_top = dolfinx.fem.locate_dofs_topological((V.sub(1), V1), fdim, facet_tags.indices[facet_tags.values == 2]) # y=200 top_tag=2
bc_top = dolfinx.DirichletBC(u_top, dofs_top, V.sub(1))
bc_u.append(bc_top)  # add bondary condition bc_left  # bc_u = [bc_left,bc_bottom,bc_top]

# define measure for surface integration
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags)
n = ufl.FacetNormal(mesh)  # normal

# Define material properties
young = 70.0e3    # 210.0e3 N/mm/mm
nu = 0.22        # 0.3
lambda_ = young/((1.0+nu)*(1.0-2.0*nu)) #  N/mm/mm
mu = young/(2.0*(1.0+nu))              #   N/mm/mm
l, Gc =0.01, 0.007                              #   N/mm
#Variational formulation
def epsilon(u):
   return 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
   return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)
# Define the traction 𝑇 , and body force 
T = dolfinx.Constant(mesh, (0, 0))       # load force  N/mm/mm
f = dolfinx.Constant(mesh, (0, 0))       # body force  N/mm/mm/mm
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)   # the trial and test function for
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
#Solve the linear variational problem
problem = dolfinx.fem.LinearProblem(a, L, bcs=bc_u, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})

uold, unew = dolfinx.Function(V)
pnew, pold, Hold = dolfinx.Function(W), dolfinx.Function(W), dolfinx.Function(W)

def psi(u):
   return 0.5*(lambda_+mu)*(0.5*(ufl.tr(epsilon(u))+abs(ufl.tr(epsilon(u)))))**2+\
          mu*ufl.inner(ufl.dev(epsilon(u)),ufl.dev(epsilon(u)))
# def psi(u):
#     return 0.5*(lambda_+mu)*(0.5*(ufl.tr(epsilon(u))+ufl.abs(ufl.tr(epsilon(u)))))**2+\
#            mu*ufl.inner(ufl.dev(epsilon(u)),ufl.dev(epsilon(u)))
def H(uold,unew):
   return ufl.conditional(ufl.lt(psi(uold),psi(unew)),psi(unew),psi(uold))
p, q = ufl.TrialFunction(V), ufl.TestFunction(V)

phi_a = (Gc*l*ufl.inner(ufl.grad(p),ufl.grad(q)) + ((Gc/l)+2.0*H(uold,unew))*ufl.inner(p,q)-2.0*H(uold,unew)*q)*ufl.dx

t = 0 # Start time
T = 0.01 # End time
num_steps = 5 # Number of time steps
dt = (T-t)/num_steps # Time step size
for i in range(num_steps):
   # Update Diriclet boundary condition 
   u_ex.t+=dt
   u_top.interpolate(u_ex.eval)
   u_top.x.scatter_forward()
   uh = problem.solve()
   force = ufl.inner(sigma(uh)* n, n)*ds(2)
   energy = psi(uh)*ufl.dx
   R = dolfinx.fem.assemble_scalar(force)
   all_e = dolfinx.fem.assemble_scalar(energy)
   print(R,all_e)

the main issue is: When setting

phi_a = (Gc*l*ufl.inner(ufl.grad(p),ufl.grad(q)) + ((Gc/l)+2.0*H(uold,unew))*ufl.inner(p,q)-2.0*H(uold,unew)*q)*ufl.dx

I got a error: ERROR:UFL:Transposed is only defined for rank 2 tensors. I can not fix it. Can someone give me an idea how to fix it ?

A second issue is: ufl.abs() is not accepted by ufl, I get: AttributeError: module ‘ufl’ has no attribute ‘abs’. I have no idea why this is happening.

Thanks!

Best wishes,
Zhihai Wang

2

You are missing the definititon of unew here, i.e.
uold, unew = dolfinx.Function(V), dolfinx.Function(V)

Secondly, as the error stated, there is no function called ufl.abs.
There is the function ufl.algebra.Abs (ufl package — Unified Form Language (UFL) 2021.1.0 documentation), but is there any reason for not using the built-in `abs function?

1 Like