Initialize in subdomain from gmsh

1

I solve diffusion equation in my program:
\frac{\rho f}{dt} \, + \, \nabla \cdot (\rho f v) \, = \, \nabla \cdot (\rho D \nabla f)
I create my geometry in gmsh:

Y-tube
Y-tube1443×366 13.9 KB

I must definition subdomain and initialize F=0 and F=1. I imported my geometry, use this code:

import meshio
msh = meshio.read("tube.msh")
def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    if prune_z:
        out_mesh.prune_z_0()
    return out_mesh

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)

line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("facet_mesh.xdmf", line_mesh)

from dolfin import * 
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)

And imported facet domains in this code:

dim = mesh.topology().dim()
mvc = MeshValueCollection("size_t", mesh, dim-1) 
with XDMFFile("facet_mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
facet_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)

My .geo file:

Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {0, 0.000150, 0, 1.0};
//+
Point(3) = {0, -0.000150, 0, 1.0};
//+
Point(4) = {-0.000750, -0.000150, 0, 1.0};
//+
Point(5) = {-0.00075, -0.000300, 0, 1.0};
//+
Point(6) = {-0.00075, 0.000150, 0, 1.0};
//+
Point(7) = {-0.00075, 0.000300, 0, 1.0};
//+
Point(8) = {0.002, 0, 0, 1.0};
//+
Point(9) = {0.002, 0.000150, 0, 1.0};
//+
Point(10) = {0.002, -0.000150, 0, 1.0};//+
Line(1) = {5, 3};
//+
Line(2) = {4, 1};
//+
Line(3) = {6, 1};
//+
Line(4) = {7, 2};
//+
Line(5) = {7, 6};
//+
Line(6) = {4, 5};
//+
Line(7) = {3, 10};
//+
Line(8) = {1, 8};
//+
Line(9) = {2, 9};
//+
Line(10) = {9, 8};
//+
Line(11) = {8, 10};
//+
Curve Loop(1) = {5, 3, -2, 6, 1, 7, -11, -10, -9, -4};
//+
Plane Surface(1) = {1};
//+
Physical Surface("sur", 12) = {1};
//+
Physical Curve("A", 13) = {5};
//+
Physical Curve("B", 14) = {6};
//+
Physical Curve("outlet", 15) = {10, 11};
//+
Physical Curve("wall", 16) = {4, 9, 7, 1, 2, 3};
//+
Curve Loop(2) = {5, 3, 8, -10, -9, -4};
//+
Plane Surface(2) = {2};
//+
Curve Loop(3) = {6, 1, 7, -11, -8, -2};
//+
Plane Surface(3) = {3};
//+
Physical Surface("sur_a", 17) = {2};
//+
Physical Surface("sur_b", 18) = {3};

I create those two surface in gmsh (sur_a, sur_b). How i can import subdomain?
My subdomain must interact each other.
I see other topic for my theme, but didn’t found answer.

2

SImply use:

dim = mesh.topology().dim()
mvc = MeshValueCollection("size_t", mesh, dim)
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mvc, "name_to_read")
surface_domains = cpp.mesh.MeshFunctionSizet(mesh, mvc)
File("surface.pvd") << surface_domains
3

Thanks, and how can initialize f field thereafter?

4

See for instance: How to define different materials / importet 3D geometry (gmsh) - #2 by dokken

5

Thanks, In your instance, you used:

V = dolfin.FunctionSpace(mesh, "DG", 0)

But in my problem f:

F = FunctionSpace(mesh, 'P', 1)

I try use “DG” elements, but my subdomain did’t interact each other
And I still don’t understand how I can extract from surface_domains the subdomain

6
  1. what do you imply when you say interact?
  2. the function you have sketched is 1 in one domain, 0 in the other. It is discontinuous across cell (i.e one data value per cell). For such data DG 0 is suitable.