from dolfin import *
from ufl import indices
import numpy as np
from dolfin import PETScKrylovSolver
nx = 10
mesh = UnitCubeMesh(nx, nx, nx)
(ε11, ε22, ε33, ε12, ε13, ε23)
V = VectorFunctionSpace(mesh, “CG”, 1, dim=6)
def as_sym_tensor(v):
return as_tensor([[v[0], v[3], v[4]],
[v[3], v[1], v[5]],
[v[4], v[5], v[2]]])
u = TrialFunction(V)
v = TestFunction(V)
u_tensor = as_sym_tensor(u)
v_tensor = as_sym_tensor(v)
def curl2(A):
i, j, k, l = indices(4)
return as_tensor(Levi[j, k, l] * A[i, l].dx(k), (i, j))
def curl1(A):
i, j, k, l = indices(4)
return as_tensor(Levi[i, k, l] * A[l, j].dx(k), (i, j))
a = inner(curl2(u_tensor), curl1(v_tensor)) * dx
zero_vec = Constant((0.0, 0.0, 0.0, 0.0, 0.0, 0.0))
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, zero_vec, boundary)
So far, is my definition correct? And which solver should I use next to solve it?