Hi @dokken ,
Thanks for responding to my question, and sorry for tagging you in a general question.
GMSH script:
tmpDir=sys.argv[1]
vals=list(map(eval,sys.argv[2:]))
rmax,S,na,nref=vals[:4]
v=array(vals[4:])
v.shape=na,3
xi=v[:,0]
yi=v[:,1]
zi=v[:,2]
#####
geoFile="%s/tmp.geo" % (tmpDir,)
gf=open(geoFile,"w")
msString="""Function Mysphere;
p1 = newp; Point(p1) = {x, y, z} ;
p2 = newp; Point(p2) = {x+r,y, z} ;
p3 = newp; Point(p3) = {x, y+r,z} ;
p4 = newp; Point(p4) = {x, y, z+r} ;
p5 = newp; Point(p5) = {x-r,y, z } ;
p6 = newp; Point(p6) = {x, y-r,z } ;
p7 = newp; Point(p7) = {x, y, z-r} ;
c1 = newreg; Circle(c1) = {p2,p1,p7};
c2 = newreg; Circle(c2) = {p7,p1,p5};
c3 = newreg; Circle(c3) = {p5,p1,p4};
c4 = newreg; Circle(c4) = {p4,p1,p2};
c5 = newreg; Circle(c5) = {p2,p1,p3};
c6 = newreg; Circle(c6) = {p3,p1,p5};
c7 = newreg; Circle(c7) = {p5,p1,p6};
c8 = newreg; Circle(c8) = {p6,p1,p2};
c9 = newreg; Circle(c9) = {p7,p1,p3};
c10 = newreg; Circle(c10) = {p3,p1,p4};
c11 = newreg; Circle(c11) = {p4,p1,p6};
c12 = newreg; Circle(c12) = {p6,p1,p7};
l1 = newreg; Line Loop(l1) = {c5,c10,c4}; Ruled Surface(newreg) = {l1};
l2 = newreg; Line Loop(l2) = {c9,-c5,c1}; Ruled Surface(newreg) = {l2};
l3 = newreg; Line Loop(l3) = {c12,-c8,-c1}; Ruled Surface(newreg) = {l3};
l4 = newreg; Line Loop(l4) = {c8,-c4,c11}; Ruled Surface(newreg) = {l4};
l5 = newreg; Line Loop(l5) = {-c10,c6,c3}; Ruled Surface(newreg) = {l5};
l6 = newreg; Line Loop(l6) = {-c11,-c3,c7}; Ruled Surface(newreg) = {l6};
l7 = newreg; Line Loop(l7) = {-c2,-c7,-c12};Ruled Surface(newreg) = {l7};
l8 = newreg; Line Loop(l8) = {-c6,-c9,c2}; Ruled Surface(newreg) = {l8};
theloops[t] = newreg;
Surface Loop(theloops[t]) = {l8+1,l5+1,l1+1,l2+1,l3+1,l7+1,l6+1,l4+1};
Return
Field[1]=ExternalProcess;
Field[1].CommandLine = "../../size.py";
Background Field = 1;"""
print ( msString,file=gf)
print ( "x=%f ; y=%f ; z=%f; r=%f ; t= %d ;" % (0,0,0,rmax,0),file=gf)
print ( "Call Mysphere;",file=gf)
print ( """
Volume (186)={theloops[]};
Physical Volume(999)= 186;""",file=gf)
print ("Mesh.MshFileVersion = 2;",file=gf)
for i,ri in enumerate(zip(xi,yi,zi)):
x,y,z=ri
print ("p11=newp ;Point(p11)={%f,%f,%f};" % (x,y,z),file=gf)
print ("Point{p11} In Volume{186};",file=gf)
gf.close()
# replace invalid double minus "--" by +
gf=open(geoFile,"r")
tmp=gf.read()
newfile=tmp.replace("--","+")
gf.close()
gf=open(geoFile,"w")
gf.write(newfile)
gf.close()
os.system("gmsh -3 %s/tmp.geo" % (tmpDir,))
for i in range(nref):
os.system("gmsh -format msh2 -refine %s/tmp.msh" % (tmpDir,))
size.py
while(True):
xyz = struct.unpack("ddd", os.read(0,24))
x,y,z=xyz
if math.isnan(x):
break
f = sizefield(x,y,z)
os.write(1,struct.pack("d",f))
sizefield.py
hm=0.1
S=load("S.npy")[0]
ni=loadtxt("nucCoords.dat")
def sizefield(x,y,z):
p=array([x,y,z])
di=[linalg.norm(p-v) for v in ni]
dm=min(di)
size=S*sqrt(dm**2+hm**2)
return size
Please know that I have tried it using both FEnics and FEnicsx, and I am getting the same results.
Input in the main code:
def convert_msh_to_xdmf(msh_file, xdmf_file):
mesh = meshio.read(msh_file)
print(f"Read mesh with {len(mesh.points)} vertices")
print(f"Found {sum(len(cells.data) for cells in mesh.cells)} cells")
meshio.write(xdmf_file, mesh)
print(f"Converted {msh_file} to {xdmf_file}")
convert_msh_to_xdmf("./tmp.msh", "mesh.xdmf")
with io.XDMFFile(comm, "mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid")
Please assist!