Including a line in a compound surface?

ArjunNarayanan commented 5 years ago

Hi All,

I am trying to specify a curve separating two regions of a compound surface as being part of the surface. A challenge that I am facing is that I cannot obtain a valid Surface ID for the compound surface so that I can use something like:

geom.add_raw_code("Line{l} In Surface{s};")

As an example, consider that I am trying to define a disk inside an annular ring, and I wish to specify the curve separating the disk and the ring as a physical line. So you can imagine the domain as an inner "core" surrounded by an outer "shell" and I wish to specify the curve separating these two regions as an "interface":

import pygmsh, meshio

geom = pygmsh.opencascade.Geometry()
d1 = geom.add_disk([0,0,0], 1, char_length = 0.1)
core = geom.add_disk([0,0,0], 0.5, char_length = 0.1)
shell = geom.boolean_difference([d1], [core], delete_other = False)

domain = geom.boolean_union([core, shell], delete_other = False)

interface_line_loop = geom.add_line_loop([core])
s = 'Line{%s} In Surface{%s};' % (interface_line_loop.id, domain.id)
geom.add_raw_code(s)

geom.add_physical_surface([shell], label = "shell")
geom.add_physical_surface([core], label = "core")
geom.add_physical_line([interface_line_loop], label = "interface")

points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom, geo_filename = "Disk.geo")
meshio.write_points_cells("disk.vtu", points, cells)

Executing this via Python gives me a gmsh syntax error. By inspecting the associated .geo file I find that the error is because domain.id of the compound surface domain is, in fact, something like bo1[] which I guess is not a valid surface id in gmsh parlance.

Perhaps my whole approach is a little too convoluted? Any tips or suggestions would be welcome!

nschloe commented 5 years ago

Hm, weird. This

// This code was created by pygmsh v4.3.6.
SetFactory("OpenCASCADE");
s0 = news;
Disk(s0) = {0, 0, 0, 1};
pts_s0[] = PointsOf{Surface{s0};};
Characteristic Length{pts_s0[]} = 0.1;
s1 = news;
Disk(s1) = {0, 0, 0, 0.5};
pts_s1[] = PointsOf{Surface{s1};};
Characteristic Length{pts_s1[]} = 0.1;
bo1[] = BooleanDifference{ Surface{s0}; Delete; } { Surface{s1}; };
bo2[] = BooleanUnion{ Surface{s1}; Delete; } { Surface{bo1[]}; };
ll0 = newll;
Line Loop(ll0) = {s1};
// Line{ll0} In Surface{bo2[]};
Physical Surface("shell") = {bo1[]};
Physical Surface("core") = {s1};
Physical Line("interface") = {ll0};

works fine in Gmsh 3.0.6. What version are you using?

ArjunNarayanan commented 5 years ago

Ah, I am on 4.0.2

nschloe commented 5 years ago

Yeah alright, looks like an incompatibility there then. Will be a while before I have 4.* running.

ArjunNarayanan commented 5 years ago

Alright! Thanks for your input.

Since It looks like I can't get gmsh v3 from the official website, do you see any workaround that might work for me? For example, can I specify a specific surface ID to be used by pygmsh?

Thanks for your time.

ArjunNarayanan commented 5 years ago

So, I managed to locate the gmsh binaries and switched to version 3.0.6

I wound up getting the same error

Error : 'Disk.geo', line 15 : syntax error (])

Then, I noticed that the offending line is commented out in your example above which is why it worked for you.

nschloe commented 5 years ago

Yes, I'm not familiar with this syntax at all. What works in gmsh again?

nschloe / pygmsh

Including a line in a compound surface? #232