[Gmsh] Bug? Points in Volume only meshed for first volume

[Gregor Mitscha-Baude]

Dear all,

I was happy to see that the new gmsh version (2.14.1) includes the 
long-missing "Point In Volume" feature.
However, if multiple volumes are created in the .geo file, with points 
added to each volume by "Point In Volume", then only the points lying in 
the *first* volume are meshed in the output .msh file.

I attached a simple example file input.geo with two cubes stacked on top 
of each other, divided by the z=0 plane. vol1 is the lower cube (z<0) 
and vol2 the upper one (z>0). A series of points of the form (0, 0, z) 
with z negative and positive is added with Point in Volume. The output 
file was obtained by running

 >> gmsh -3 input.geo -o out.msh

and is also attached. As can be seen, all the embedded points with 
positive z values are missing from the .msh file. Hence, only the points 
embedded in vol1 were meshed. I observed this also on more complicated 
examples with many different volumes: only the points embedded in volume 
1 are meshed.
Is this a bug? Is there any workaround to mesh embedded points in 
multiple volumes?

Cheers, and thanks a lot,
Gregor Mitscha-Baude
-------------- next part --------------
A non-text attachment was scrubbed...
Name: out.msh
Type: model/mesh
Size: 9323 bytes
Desc: not available
URL: <http://onelab.info/pipermail/gmsh/attachments/20161114/9b134cff/attachment.msh>
-------------- next part --------------
// This file was created by python4gmsh v0.1.
// Copyright (c) 2013, Nico Schlömer <nico.schloemer at gmail.com>
// All rights reserved.
//
// The latest updates can be retrieved from
//     https://github.com/nschloe/python4gmsh
// where you can also make suggestions 
// and help improve python4gmsh.
//
p1 = newp;
Point(p1) = {-1, -1, -1, 1};
p2 = newp;
Point(p2) = {-1, 1, -1, 1};
l1 = newl;
Line(l1) = {-p1, p2};
p3 = newp;
Point(p3) = {1, -1, -1, 1};
p4 = newp;
Point(p4) = {1, 1, -1, 1};
l2 = newl;
Line(l2) = {-p3, p4};
l3 = newl;
Line(l3) = {-p1, p3};
l4 = newl;
Line(l4) = {-p2, p4};
ll1 = newll;
Line Loop(ll1) = {-l1,l2,l3,-l4};
surf1 = news;
Plane Surface(surf1) = {ll1};
p5 = newp;
Point(p5) = {-1, -1, 0, 1};
p6 = newp;
Point(p6) = {-1, 1, 0, 1};
l5 = newl;
Line(l5) = {-p5, p6};
p7 = newp;
Point(p7) = {1, -1, 0, 1};
p8 = newp;
Point(p8) = {1, 1, 0, 1};
l6 = newl;
Line(l6) = {-p7, p8};
l7 = newl;
Line(l7) = {-p5, p7};
l8 = newl;
Line(l8) = {-p6, p8};
ll2 = newll;
Line Loop(ll2) = {-l5,l6,l7,-l8};
surf2 = news;
Plane Surface(surf2) = {ll2};
l9 = newl;
Line(l9) = {-p1, p5};
l10 = newl;
Line(l10) = {-p3, p7};
ll3 = newll;
Line Loop(ll3) = {-l9,l10,l3,-l7};
surf3 = news;
Plane Surface(surf3) = {ll3};
l11 = newl;
Line(l11) = {-p2, p6};
l12 = newl;
Line(l12) = {-p4, p8};
ll4 = newll;
Line Loop(ll4) = {-l11,l12,l4,-l8};
surf4 = news;
Plane Surface(surf4) = {ll4};
ll5 = newll;
Line Loop(ll5) = {-l9,l11,l1,-l5};
surf5 = news;
Plane Surface(surf5) = {ll5};
ll6 = newll;
Line Loop(ll6) = {-l10,l12,l2,-l6};
surf6 = news;
Plane Surface(surf6) = {ll6};
surfloop1 = newsl;
Surface Loop(surfloop1) = {-surf1,+surf2,+surf3,-surf4,-surf5,+surf6};
vol1 = newv;
Volume(vol1) = surfloop1;
p9 = newp;
Point(p9) = {-1, -1, 1, 1};
p10 = newp;
Point(p10) = {-1, 1, 1, 1};
l13 = newl;
Line(l13) = {-p9, p10};
p11 = newp;
Point(p11) = {1, -1, 1, 1};
p12 = newp;
Point(p12) = {1, 1, 1, 1};
l14 = newl;
Line(l14) = {-p11, p12};
l15 = newl;
Line(l15) = {-p9, p11};
l16 = newl;
Line(l16) = {-p10, p12};
ll7 = newll;
Line Loop(ll7) = {-l13,l14,l15,-l16};
surf7 = news;
Plane Surface(surf7) = {ll7};
l17 = newl;
Line(l17) = {-p5, p9};
l18 = newl;
Line(l18) = {-p7, p11};
ll8 = newll;
Line Loop(ll8) = {-l17,l18,l7,-l15};
surf8 = news;
Plane Surface(surf8) = {ll8};
l19 = newl;
Line(l19) = {-p6, p10};
l20 = newl;
Line(l20) = {-p8, p12};
ll9 = newll;
Line Loop(ll9) = {-l19,l20,l8,-l16};
surf9 = news;
Plane Surface(surf9) = {ll9};
ll10 = newll;
Line Loop(ll10) = {-l17,l19,l5,-l13};
surf10 = news;
Plane Surface(surf10) = {ll10};
ll11 = newll;
Line Loop(ll11) = {-l18,l20,l6,-l14};
surf11 = news;
Plane Surface(surf11) = {ll11};
surfloop2 = newsl;
Surface Loop(surfloop2) = {-surf2,+surf7,+surf8,-surf9,-surf10,+surf11};
vol2 = newv;
Volume(vol2) = surfloop2;
Physical Volume(1) = vol1;
Physical Volume(2) = vol2;
Physical Surface(1) = {surf1,surf7,surf3,surf8,surf4,surf9,surf5,surf10,surf6,surf11};
p13 = newp;
Point(p13) = {0, 0, -0.8, 1};
Point{p13} In Volume{vol1};
p14 = newp;
Point(p14) = {0, 0, -0.6, 1};
Point{p14} In Volume{vol1};
p15 = newp;
Point(p15) = {0, 0, -0.4, 1};
Point{p15} In Volume{vol1};
p16 = newp;
Point(p16) = {0, 0, -0.2, 1};
Point{p16} In Volume{vol1};
p17 = newp;
Point(p17) = {0, 0, 0.2, 1};
Point{p17} In Volume{vol2};
p18 = newp;
Point(p18) = {0, 0, 0.4, 1};
Point{p18} In Volume{vol2};
p19 = newp;
Point(p19) = {0, 0, 0.6, 1};
Point{p19} In Volume{vol2};
p20 = newp;
Point(p20) = {0, 0, 0.8, 1};
Point{p20} In Volume{vol2};
General.ExpertMode = 1;