from contextlib import ExitStack
import numpy as np
#from dolfinx import *
from dolfinx import la
from dolfinx.cpp.function import near
from dolfinx.fem import (Expression, Function, FunctionSpace,
VectorFunctionSpace, dirichletbc, form,
locate_dofs_topological, Constant)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import (CellType, GhostMode, create_box, create_rectangle,
locate_entities_boundary)
import ufl
from ufl import dx, grad, inner, sym, Identity, TrialFunction, TestFunction
from mpi4py import MPI
from petsc4py import PETSc
comm = MPI.COMM_SELF
L = 25.
H = 1.
Nx = 250
Ny = 10
mesh = create_rectangle(comm=comm, points=((0.0, 0.0), (L, H)), n=(Nx, Ny),
cell_type=CellType.triangle)
def eps(v):
return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)
#def sigma(v):
# return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def sigma(v):
return lmbda*eps(v)*Identity(2) + 2.0*mu*eps(v)
rho_g = 1e-3
f = Constant((0, -rho_g))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx
def left(x, on_boundary):
return near(x[0], 0.)
bc = dirichletbc(V, Constant((0.,0.)), left)
u = Function(V, name="Displacement")
#solve(a == l, u, bc)
problem = PETSc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
#plot(1e3*u, mode="displacement")
print('done...')
I moved this code into dolfinx-0.6.0 from:
https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/2D_elasticity.py.html
So far there seems to be a difficulty with import dolfinx.cpp.funtion.near. I see there is a way to implement a linear solver from fenicsx tutorials if that is the same type of solver as was meant by the author of 2D_elasticity.html tutorial. So far the function tr() is undefined. Could it be stated what that tr() was originally meant to compute it?