How to apply particular bc on imported mesh?


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)
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)

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]


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”

Thank you in advance for your time

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.

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

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