Hello,
I’m having some problems working with Fenics solver and I don’t know ho to solve it (I think that is more due to syntaxe). I have the following code : (minimum to reproduce the error)
from __future__ import print_function
from dolfin import *
from mshr import *
from dolfin import File
import matplotlib.pyplot as plt
import numpy as np
from dolfin import plot
# Initialize circular subdomain
small_square = Rectangle(Point(-0.1, -0.1), Point(0.1, 0.1))
domain1 = Rectangle(Point(-0.5, -0.5), Point(0.5, 0.5)) - Rectangle(Point(-0.5, -0.0001), Point(0, 0.0001))
domain = domain1 - small_square
mesh=generate_mesh(domain,25)
mesh = refine(mesh)
class SubDomainBoundary(SubDomain):
def inside(self,x,on_boundary):
return (near(x[0], -0.1) or near(x[0], 0.1) or near(x[1], -0.1) or near(x[1], 0.1)) and on_boundary
# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
# Introduce manually the material parameters
Gc = 2.7
l = 0.015
lmbda = 121.1538e3
mu = 80.7692e3
# Constituive functions
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))
def psi(u):
return 0.5*(lmbda+mu)*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2+mu*inner(dev(epsilon(u)),dev(epsilon(u)))
def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)
# Boundary conditions and crack propagation
sub_domain_boundary = SubDomainBoundary()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
sub_domain_boundary.mark(boundaries, 1)
load = Expression("t", t = 0.0, degree=1)
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(1), load, top)
bc_u = [bcbot, bctop]
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
neumann_bc = Expression(("0.0", "0.0"), degree=1)
def Crack(x):
return abs(x[1]) < 1e-03 and x[0] <= 0.0
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
# Variational form
unew, uold = Function(W), Function(W)
pnew, pold, Hold = Function(V), Function(V), Function(V)
E_du = ((1.0-pold)**2)*inner(grad(v),sigma(u))*dx
E_phi = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
*inner(p,q)-2.0*H(uold,unew,Hold)*q)*dx + inner(neumann_bc, q) * ds(1)
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew, bc_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)
And I got the following error:
Shapes do not match: <Coefficient id=127173665399296> and <Argument id=127173241947808>.
Traceback (most recent call last):
File "/home/jandui/Documents/MATH code/First FEM approx/brittle_csquare.py", line 79, in <module>
*inner(p,q)-2.0*H(uold,unew,Hold)*q)*dx + inner(neumann_bc, q) * ds(1)
File "/usr/lib/python3/dist-packages/ufl_legacy/operators.py", line 158, in inner
return Inner(a, b)
File "/usr/lib/python3/dist-packages/ufl_legacy/tensoralgebra.py", line 147, in __new__
error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))
File "/usr/lib/python3/dist-packages/ufl_legacy/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl_legacy.log.UFLException: Shapes do not match: <Coefficient id=127173665399296> and <Argument id=127173241947808>.
Why is not correct do the following:
+ inner(neumann_bc, q) * ds(1)
I saw an unsolved problem in the forum and I’m trying to addapt. Is it correct the way I define the subomain and neumann zero boundary conditions in the subdomain?