Problem with surfboolean in python

[Scott Pesenti]

Hi doctor Fang,

I am having some trubles with the function surfboolean trying to unite two separate meshes. 

The code correctly generates the two independent meshes, but when it unites them something goes astray.

Here is the Python code with the outputted mesh:

import iso2mesh as i2m #tool per lavorare con le mesh
import pmmc #ha bisogno del pacchetto sparse-numba
import numpy as np #array
import matplotlib.pyplot as plt #per plottare
import pmcx

#Geometrical paramters (NB all quantities are mm based)
t1 = 100 #height
c0 = np.array([0, 0, 0]) #center of the face below
c1 = np.array([0, 0, t1]) #center of the face above
r1 = 10 #radius
p0 = np.array([-0.5, -r1, 0]) #edge of the box (min)
p1 = np.array([-t1, r1, t1]) #edge of the box (max)

#parameters medium 1
att = 2 #thickness of the layer
mua1 = 0
mus_red1 = 1 #NB it's mus'
g1 = 0
mus1 = mus_red1/(1 - g1) #convertion of mus' into mus for the calculations
n1 = 1.40

#parameters medium 2
mua2 = 0.05
mus_red2 = 1
g2 = 0
mus2 = mus_red2/(1 - g2)
n2 = 1.40

#mesh parameters
r2 = r1 - att
tsize = 5 #maximum area of the triangles
maxvol = 5 #maximum volume of the tethraedra
ndiv = 40 #number of edges of the two faces
keepratio = 1 #how many of the surface triangles we are keeping between the operation steps
ext1 = np.array([-t1, -r1, 0]) #edge of the volume (min)
ext2 = np.array([r1, r1, t1]) #edge of the volume (max)

#creation of the mesh
node1, face1 , elem1 = i2m.meshacylinder(c0, c1, r1, tsize, maxvol, ndiv) #cylinder
#node1, face1 = i2m.removeisolatednode(node1, face1)
node2, face2, elem2 = i2m.meshabox(p0, p1, tsize, maxvol) #parallelepipe
#node2, face2 = i2m.removeisolatednode(node2, face2)

newnode, newface = i2m.surfboolean(node1, face1, 'union', node2, face2) #merge the surfaces
newnode, newface = i2m.meshcheckrepair(newnode, newface, 'dup') #correction of eventual errors
newnode, newface = i2m.remeshsurf(newnode, newface, 1)

newnode, newelem, newface = i2m.surf2mesh(newnode, newface, ext1, ext2, keepratio, maxvol) #creation of the volume mesh from the surface one
i2m.plotmesh(newnode, newface)
plt.show()
#distinction of the two layers
centers = np.mean(newnode[newelem[:, :4]-1], axis=1)
r = np.linalg.norm(centers[:, :2], axis=1) #average radial position of each tethraedra
newelem[:, 4] = 1
newelem[r < r2, 4] = 2 #each internal point is labeled as medium 2