How to apply particular bc on imported mesh?

1

HI,

Please i would like to know how to apply pointwise bc and a gaussian load over a subdomain on a generated mesh with gmsh.

Here is a code which reproduce my error:


from dolfin import *
from mshr import *
from math import *
import numpy as np
import sympy
import ufl

# mechanical parameters
E, nu = Constant(210e9), Constant(0.3)
mu, lmbda = E/2/(1+nu), E*nu/(1+nu)/(1-2*nu)
def sigma(v):
	return lmbda*tr(eps(v))*Identity(2)+2.0*mu*eps(v)
def eps(v):
	return sym(grad(v))

# mesh import and principal parameters
mesh = Mesh("gmsh_cube.xml")
mm = 10^(-3)
L = 20*mm
H = 20*mm
mean_mu = 0.5
standard_deviation = 1/np.sqrt(2*pi)
rho_g = 4400*9.81

# subdomains
class upper_set(SubDomain):
	def inside(self, x, on_boundary):
		return (x[0]>L/2-5*standard_deviation) and (x[0]<L/2+5*standard_deviation)
class lower_set_left(SubDomain):
	def inside(self, x, on_boundary):
		return near(x[0],0.,DOLFIN_EPS) and near(x[1],0.,DOLFIN_EPS)
class lower_set_right(SubDomain):
        def inside(self, x, on_boundary):
                return near(x[0],L,DOLFIN_EPS) and near(x[1],0.,DOLFIN_EPS)
boundary_top, boundary_left, boundary_right = upper_set(), lower_set_left, lower_set_right
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
subdomains.set_all(0)
boundary_top.mark(subdomains, 1)
ds=Measure('ds', domain=mesh, subdomain_data=subdomains)

# variational problem and bc
f = Constant((0,-rho_g))
g = Expression(("0","(1/(sqrt(2*pi*standard_deviation)))*exp(- pow((x[0]-mean_mu),2)/(2*pow(standard_deviation,2)))"), degree=0, mean_mu=mean_mu, standard_deviation=standard_deviation)

V=VectorFunctionSpace(mesh,'CG',degree=2)
u = Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du),eps(u_))*dx
l = inner(f,u_)*dx+inner(g,u_)*ds(1)

bc_left = DirichletBC(V, Constant((0.,0.)), boundary_left, method="pointwise")
bc_right = DirichletBC(V.sub(1), Constant(0.), boundary_right, method="pointwise")
bc_u = [bc_left, bc_right, boundary_top]

solve(a==l,u,bc_u)

and here is the gmsh code :

//cube 20x20

mm = 10^(-3);
mesh_size = 2*mm;

L = 20*mm;
H = 20*mm;

Point(1) = {0, 0, 0, mesh_size};
Point(2) = {L, 0, 0, mesh_size};
Point(3) = {L, H, 0, mesh_size};
Point(4) = {0, H, 0, mesh_size};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

Line Loop(5) = {1, 2, 3, 4};

Plane Surface(6) = {5};

Physical Surface(7) = {6};

Here is how i proceed from the gmsh file : 1. name the gmsh code as “gmsh_cube.geo” , 2. use the commands “gmsh -2 -format msh22 gmsh_cube.geo” and “dolfin-convert gmsh_cube.msh gmsh_cube.xml” to mesh then to convert it on xml 3. run the fenics code “python code.py”

Thank you in advance for your time

3

Your first issue is the code quoted above, as it has to be changed to:

boundary_top, boundary_left, boundary_right = upper_set(), lower_set_left(), lower_set_right()

The second issue is that you are trying to send boundary_top into the solver as a DirichletBC, when it is a SubDomain.

4

Hi, thank you very much for your answer. The first issue is a typo when i tried to make a minimal code to represent my problem, i forgot to write the “()”.

Please could you tell me 1. how to correctly apply my particular “g” load over my subdomain then ? it will be correct if i just write “bc_u = [bc_left, bc_right]” and 2. are “bc_left” and “bc_right” correctly defined without taking into account the “()” that i forgot ?

Thank you for you time

5

As the top boundary condition is imposed through a Neumann condition,

you only need

bc_u = [bc_left, bc_right]

with the corrections (adding () as described above).
Please note that you could use GMSH to mark boundaries, as described in Transitioning from mesh.xml to mesh.xdmf, from dolfin-convert to meshio - #181 by rkailash

1 Like