Hi,
I’m trying to run a simple 2D problem on fenics based on the 3D example here. Currently, my code works because I’ve set all the z-components for the boundary conditions and traction values to 0. However, I am not sure if this is the best way to go about this problem. Currently, if I change the Dirichlet condition to a 2D vector I get the following:
*** Error: Unable to create Dirichlet boundary condition.
*** Reason: Illegal value dimension (2), expecting (3).
*** Where: This error was encountered inside DirichletBC.cpp.
A similar problem arises if I set the traction force as a 2D vector:
ufl.log.UFLException: Dimension mismatch in dot product.
MWE:
#!/usr/bin/env python3
from fenics import *
from ufl import nabla_div
from ufl import nabla_grad
import numpy as np
import os
import meshio
mesh_file = './mesh' #Mesh folder name
model_name = 'Half_cy' #Mesh model to load
msh = meshio.read(mesh_file+'/'+model_name+'.msh')
#read mesh
meshio.write(mesh_file+'/'+model_name+'.xdmf', meshio.Mesh(points=msh.points, cells={"triangle":msh.cells["triangle"]}))
#mesh function (mf) to pick the boundaries (Pick up line)
meshio.write(mesh_file+'/mf'+model_name+'.xdmf', meshio.Mesh(points=msh.points, cells={"line":msh.cells["line"]}, \
cell_data={"line":{"name_to_read":msh.cell_data["line"]["gmsh:physical"]}}))
#cell function (cf) to pick out the different regions (Pick up surface)
meshio.write(mesh_file+'/cf'+model_name+'.xdmf', meshio.Mesh(points=msh.points, cells={"triangle":msh.cells["triangle"]}, \
cell_data={"triangle":{"name_to_read":msh.cell_data["triangle"]["gmsh:physical"]}}))
mesh = Mesh()
with XDMFFile(mesh_file+'/'+model_name+'.xdmf') as infile:
infile.read(mesh) #Feeds data from xdmf to mesh variable
mvc = MeshValueCollection("size_t", mesh, 1) #Container for mesh values
with XDMFFile(mesh_file+'/mf'+model_name+'.xdmf') as infile:
infile.read(mvc, "name_to_read") #Feeds boundary data to mvc
mvc2 = MeshValueCollection("size_t", mesh, 1) #Container for mesh values
with XDMFFile(mesh_file+'/cf'+model_name+'.xdmf') as infile:
infile.read(mvc2, "name_to_read") #Feeds region data into mvc2
boundaries = MeshFunction("size_t", mesh, mvc) #Contains information regarding the boundaries
regions = MeshFunction("size_t", mesh, mvc2) #Constains information regarding the regions
V = VectorFunctionSpace(mesh, 'P', 1) #Vector function space for the equations
L = 1
W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma
#Based on gmsh Line and sufurce marker
Fixed_boundary_marker = 1
Traction_boundary_marker = 2
Region_marker = 3
boundary_conditions = {
1 : {"Dirichlet" : Constant((0, 0, 0))} #Should it be 3D?
}
bcs = []
for i in boundary_conditions:
if "Dirichlet" in boundary_conditions[i]:
#Here V.sub(1) refers to pressure while 0 refers to velocity
bc = DirichletBC(V, boundary_conditions[i]["Dirichlet"], boundaries, i)
bcs.append(bc)
#Strain
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
#Stress
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Constant((0, rho*g, 0.0))
T = Constant((0, 0, 0.0)) # Should it be 3D???
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds
# Compute solution
u = Function(V)
solve(a == L, u, bcs)
# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
print('min u: '+str(u_magnitude.vector()[:].min()))
print('max u: '+str(u_magnitude.vector()[:].max()))
The gmsh geo file:
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0.5, 0, 0, 1.0};
//+
Point(3) = {1.5, 0, 0, 1.0};
//+
Point(4) = {2.0, 0, 0, 1.0};
//+
Point(5) = {1.0, 0, 0, 1.0};
//+
Circle(1) = {3, 5, 2};
//+
Circle(2) = {4, 5, 1};
//+
Line(3) = {1, 2};
//+
Line(4) = {3, 4};
//+
Curve Loop(1) = {1, -3, -2, -4};
//+
Plane Surface(1) = {1};
//+
Physical Curve("1") = {3};
//+
Physical Curve("2") = {4};
//+
Physical Surface("3") = {1};
Worth noting that the gmsh version is: 4.4.1
Thanks
