GMSH: Defining BC for physical volume and applying load on edge of cantilever beam mesh: Linear Elasticity

Question

Hi,

I have been though the following link and even participated on the discussion in there, but I am confused as to how to enforce boundary condition on physical volume.

http://fenicsproject.org/qa/2986/how-to-define-boundary-condition-for-mesh-generated-by-gmsh

So, I have a unit cube in gmsh with the following physical surfaces and volumes:

Physical Surface(0) = {6};     /top face
Physical Surface(1) = {23};   /bottom face
Physical Surface(2) = {28};   /left (fixed)
Physical Surface(3) = {15};   /right..so on and so forth
Physical Surface(4) = {27}; 
Physical Surface(5) = {19};
Physical Volume(0) = {1};    /Volume of cube
Physical Line(6) = {34};      /Upper edge of right face of cantilever beam 

I believe that to impose Drichlet BC on a physical surface and volume, I can do the following:

mesh = Mesh("beam.xml")
boundaries = MeshFunction("size_t", mesh, "beam_facet_region.xml")
//Should I use Facetfunction in the above statement?

subdomains = MeshFunction("size_t", mesh, "beam_physical_region.xml")
//Should I use CellFunction in the above statement?

bc1 = DirichletBC(V, u1, boundaries, 2)   #boundary condition for boundary (face of cube)
bc2 = DirichletBC(V, u2, subdomains, 0) #boundary condition for volume 

1. Is this the correct way?

2. Please shed some light on inbuilt function DomainBoundary() in Fenics.
   I saw this being implemented in linear elasticity (Fig 26.2) of Fenics Book.

   bc2 = DirichletBC(V, u2, DomainBoundary())

Does this return the entire boundary of the domain, i.e, all edges in case of 2D and the entire volume in case of 3D)?

1. If I want to impose a load on the top EDGE of the beam (which I have defined as a physical edge above); can I code it like this?

   boundaries = MeshFunction("size_t", mesh, "beam_facet_region.xml")
       ds = Measure("ds")[boundaries]
       L = - inner(g_T, v)*ds(6)    <---- Pls. Note (6) here
       a = inner(grad(v),sigma(U))*dx
   

Comments

Sorry to revive an old thread. When importing a mesh from GMSH using dolfin-convert, it looks like the requirement of manually renumbering physical surfaces and volumes is not necessary anymore, as it was implied in the link above. Is this correct?

Answer 1

Hello,

Here is how I do this with gmsh. Maybe it helps you. It is a two-dimensional code, however.

To apply neumann conditions (i.e. forces), I do

 boundaries = MeshFunction("uint", mesh, mesh_input_file_base          +"_facet_region.xml     ")                                                                                               
force = Constant((1.0, 1.0))                 
# '4' is the number from the gmsh .geo file                    
L =  inner(force, w)*ds(4, domain=mesh,  subdomain_data=boundaries)

For Dirichlet conditions, your approach seems ok.

Is there a reason, why you want to apply boundary conditions to a volume. This is ok from a mathematical point of view, but I'm not sure why you want to do it. Also, in your example .geo file, you only define a single volume, i.e. the whole cube. So you would be applying specific displacements to the whole cube, which makes using a finite element code kindof pointless, because there is nothing left to solve.

In my opinion, it would make sense to first figure out if you actually need volume boundary conditions. I sometimes do, but never from gmsh input. If you decide that you realy need them (from gmsh input), then you should try to define something other than the whole cube as your volume. You should be able to simply open your .geo file in gmsh, add a few faces in one of the corners and then use those to define a volume.

If you need more help, I will try. But firstI suggest that you split your problem up into several minimum examples, each containing only one of the problems that you are working on, i.e. one example for neumann conditions, one example for volume conditions, etc. If you then post those files here, it will be much easier for me to help.

regards

Comments

Thanks a lot multigrid202 for taking the time to help me out!

No, I am not applying any BC on Volume. As you said, It wouldn't make any sense in case of just 1 volume. I just wanted to ensure what I thought was the right way to code it.

It would be awesome if you could guide on the second question as well and on the following too:

boundaries = MeshFunction("size_t", mesh, "beam_facet_region.xml")
//Should I use Facetfunction in the above statement?

subdomains = MeshFunction("size_t", mesh, "beam_physical_region.xml")
//Should I use CellFunction in the above statement?

---

Hi there, no problem
I think I answered question (2) in one of you other posts. The short version: I'm not 100% sure, but it seems that DomainBoundary() returns the boundary (an object of type SubDomain), i.e. all the edges in case of 2D and all the faces in case of 3D.

For the boundaries, MeshFunction() works fine for me.

For the volumes, unfortunately I don't know the answer, but I also think that it is unlikely that you ever need it. If you need it, then I need a minimum code example to further help you. That means I would need your .geo file and the fenics code.

regards

---

Hello

for some reason the formatting of the webpage cuts off your code. Can you please post it again here?

Also, what exactly are you trying to achieve? If I understand you correctly, then you want to find out if you can apply boundary conditions to a volume, is that correct?

---

Hi multigrid202,

The following are my python code and .geo file for cantilever beam.

https://github.com/ChaitanyaGoyal/3D_Linear_Elasticity_Cantilever_Beam

I came across THIS example of solving Poisson eqn. with subdomains:

From what I learnt from the variational formulation shown above, for the elasticity problem, I would need to define boundaries as 'ds' on the RHS, and tag them to apply force or displacement selectively on a boundary. So we would represent this in variational form as:

a = dot(grad(v), sigma(U))*dx
L = dot(v, f)*ds(*gmsh tag*)

I understand that DOLFIN predefines the “measures” dx, ds and dS representing integration over cells, exterior facets (that is, facets on the boundary) and interior facets, respectively.

I used this approach successfully for a beam mesh defined in Fenics itself (code posted in repository). However, things are not clear when importing the mesh from gmsh.

2) To apply a point load on free end -- I defined one of the corner vertices on the free end of beam as a 'physical point' in gmsh.

I also defined the edge on free end as a 'physical line' to try loading an edge. Now, according to what I understand, this should have worked, but it doesn't.

ds = Measure("ds")[boundaries]
a = inner(grad(v),sigma)dx
L= dot(v,f)ds(point tag) <-- i.e. the number for point tag or edge tag

To sum up the story, I want to understand how to successfully imply load on a point, or an edge of the cantilever beam when the code is imported from gmsh. All codes including geo file are posted. Thanks a lot!!!!

---

I know my last comment was very confusing and lengthy. Apologize fr that. I have edited it to express the problem precisely now. Thanks!!

---

hello, sorry I'm taking some time but need to finish some project work at uni atm. I will get back to you when that is done. Sorry again for the delay.