Using Gmsh for defining a mesh.

5 months ago by
I am using Gmsh to define a cylindrical solid with a mesh on it. I used the following .geo file:
// Gmsh project created on Fri Feb 23 19:03:54 2018
Cylinder(1) = {0, 0, 0, 0, 0, 1, 1, 2*Pi};
Physical Surface(1) = {1, 2, 3};
Physical Volume(2) = {1};​

then I used the command "dolfin-convert cylinder.msh cylinder.xml". After that I used the following code:

mesh = Mesh('cylinder.xml')
V = VectorFunctionSpace(mesh, 'P', 1)

around =  CompiledSubDomain('near(x[0]*x[0]+x[1]*x[1] ,1) && on_boundary')  
up = CompiledSubDomain('near(x[2], side) && on_boundary', side = 1.0)
down = CompiledSubDomain('near(x[2], side) && on_boundary', side = 0.0)

in_mesh = MeshFunction('size_t' , mesh , 'cylinder_physical_region.xml')  
bnd_mesh = MeshFunction('size_t', mesh , 'cylinder_facet_region.xml') 

around.mark(bnd_mesh, 1)
down.mark(bnd_mesh, 3)  
up.mark(bnd_mesh, 2)

ds = Measure('ds', subdomain_data=bnd_mesh)
dx = Measure('dx', subdomain_data=in_mesh)

v  = TestFunction(V)             
u  = Function(V) 

T  = Expression(('-1*x[0]/sqrt(x[0]*x[0]+x[1]*x[1])' , 
                '-1*x[1]/sqrt(x[0]*x[0]+x[1]*x[1])' , 
                 '0.0'), degree=1)

Pi = dot(T,u)*ds(1)

Suppose I just want to integrate dot(T,u) on the surface of the cylinder excluding top and bottom. Are my subdomain definitions correct? Is Pi a surface integral that is applied only on the surface without top and bottom?

Community: FEniCS Project
Just as an info: you can plot your meshfunctions or write them to files for debugging purposes. Then you can check if your surfaces have the right marker.
written 5 months ago by klunkean  
Thank you for your reply.
Sorry for asking trivial questions. But I am new to this.
Could you please tell me how I can plot bnd_mesh for example?
written 5 months ago by Nima  
Apparently you can't plot a MeshFunction directly using the built-in plot functionality. But consider this example:
from dolfin import *

mesh = UnitCubeMesh(4,4,4)

class LeftHalf(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]<=.5
class RightHalf(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]>=.5

left = LeftHalf()
right = RightHalf()

domains = MeshFunction("size_t", mesh, 3)
left.mark(domains, 1)
right.mark(domains, 2)

File('meshfun.pvd') << domains​

The file meshfun.pvd can be viewed in ParaView for example. The relevant line for you is

File('meshfun.pvd') << domains​
written 5 months ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »