Add fem simulation files

This commit is contained in:
Sudarsan Balakrishnan 2026-03-12 12:27:03 -04:00
parent cbd4bf42a7
commit 9aced93bec
12 changed files with 680 additions and 0 deletions

2
anasen_fem/clean.sh Normal file
View File

@ -0,0 +1,2 @@
rm wires2d.msh
rm wires2d/*

114
anasen_fem/junk/wires.py Normal file
View File

@ -0,0 +1,114 @@
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
acolors = plt.get_cmap('tab20',24)
ccolors = plt.get_cmap('tab20',24)
k=-2*np.pi/24.
offset = 6*k + 3*k #-pi/2
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
labelsa_1 = np.array([i for i in np.arange(0,24)])
fig,ax = plt.subplots(figsize=(10,10))
ax.invert_yaxis()
ax.plot(xarra_1,yarra_1,"x",label="anode, z=-L/2")
for x,y,label in zip(xarra_1,yarra_1,labelsa_1):
ax.text(x,y,label)
kc=2*np.pi/24.
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
labelsc_1 = np.array([i for i in np.arange(0,24)])
ax.plot(xarrc_1,yarrc_1,"o",label="cathode, z=-L/2, where they are picked up")
for x,y,label in zip(xarrc_1,yarrc_1,labelsc_1):
ax.text(x,y,label)
plt.title("z=-L/2 plane, beam going into the plane along +z, (+x right, +y down)")
plt.grid()
plt.legend()
plt.savefig("plane1.png")
plt.show()
fig,ax = plt.subplots(figsize=(10,10))
ax.invert_yaxis()
offset = offset-3*k
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
labelsa_2 = np.array([i for i in np.arange(0,24)])
ax.plot(xarra_2,yarra_2,"x",label="anode, z=+L/2, where they are picked up")
for x,y,label in zip(xarra_2,yarra_2,labelsa_2):
ax.text(x,y,label)
offsetc = offsetc-3*kc
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
labelsc_2 = np.array([i for i in np.arange(0,24)])
ax.plot(xarrc_2,yarrc_2,"o",label="cathode, z=+L/2")
for x,y,label in zip(xarrc_2,yarrc_2,labelsc_2):
ax.text(x,y,label)
plt.title("z=+L/2 plane, beam going into the plane along +z, (+x right, +y down)")
plt.grid()
plt.legend()
plt.savefig("plane2.png")
plt.show()
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111,projection='3d')
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
xx_a = np.array([[x1,x2] for x1,x2 in zip(xarra_1,xarra_2)])
yy_a = np.array([[y1,y2] for y1,y2 in zip(yarra_1,yarra_2)])
zz_a = np.array([[-173.5,173.5] for x1,x2 in zip(xarra_1,xarra_2)])
for i,[xx,yy,zz] in enumerate(zip(xx_a,yy_a,zz_a)):
ax.plot(xx,yy,zz,'-',color=acolors(i/24))
for i,[x,y,label] in enumerate(zip(xarra_1,yarra_1,labelsa_1)):
ax.text(x,y,-173,"a"+str(label),color=acolors(i/24))
for i,[x,y,label] in enumerate(zip(xarra_2,yarra_2,labelsa_2)):
ax.text(x,y,+173,"a"+str(label),color=acolors(i/24))
xx_c = np.array([[x1,x2] for x1,x2 in zip(xarrc_1,xarrc_2)])
yy_c = np.array([[y1,y2] for y1,y2 in zip(yarrc_1,yarrc_2)])
zz_c = np.array([[-173.5,173.5] for x1,x2 in zip(xarrc_1,xarrc_2)])
for i,[xx,yy,zz] in enumerate(zip(xx_c,yy_c,zz_c)):
ax.plot(xx,yy,zz,'--',color=ccolors(((47-i)%24)/24))
for i,[x,y,label] in enumerate(zip(xarrc_1,yarrc_1,labelsc_1)):
ax.text(x,y,-173,"c"+str(label),color=ccolors(((25-i)%24)/24))
for i,[x,y,label] in enumerate(zip(xarrc_2,yarrc_2,labelsc_2)):
ax.text(x,y,+173,"c"+str(label),color=ccolors(((47-i)%24)/24))
ax.view_init(elev=-53, azim=-106, roll=18)
plt.tight_layout()
plt.show()
phi_qqq = np.array([[2*np.pi*(-i*16+j+0.5)/(16*4) for i in range(4)] for j in range(16)])
print(phi_qqq)
#'''
for i,[phi1,phi2,phi3,phi4] in enumerate(phi_qqq):
ax.plot([50*np.cos(phi1),100*np.cos(phi1)],[50*np.sin(phi1),100*np.sin(phi1)],[100,100],'-',color='red')
ax.text(104*np.cos(phi1),104*np.sin(phi1),100,"0_%d"%(i),color="red")
ax.plot([50*np.cos(phi2),100*np.cos(phi2)],[50*np.sin(phi2),100*np.sin(phi2)],[100,100],'-',color='green')
ax.text(104*np.cos(phi2),104*np.sin(phi2),100,"1_%d"%(i),color="green")
ax.plot([50*np.cos(phi3),100*np.cos(phi3)],[50*np.sin(phi3),100*np.sin(phi3)],[100,100],'-',color='blue')
ax.text(104*np.cos(phi3),104*np.sin(phi3),100,"2_%d"%(i),color="blue")
ax.plot([50*np.cos(phi4),100*np.cos(phi4)],[50*np.sin(phi4),100*np.sin(phi4)],[100,100],'-',color='brown')
ax.text(104*np.cos(phi4),104*np.sin(phi4),100,"3_%d"%(i),color="brown")
#'''
#coords_qqq = np.array([[50*np.cos(phi),100*np.sin(phi),100] for phi in phi_qqq]).T
plt.tight_layout()
plt.show()

View File

@ -0,0 +1,71 @@
Check Keywords Warn
Header
Mesh DB "." "wires2d"
End
Simulation
Coordinate System = Cartesian 2D
Simulation Type = Steady State
Steady State Max Iterations = 1
Output File = "elstatics.result"
Post File = "elstatics.ep"
End
Constants
Permittivity Of Vacuum = 8.8542e-12
End
Body 1
Target Bodies(1) = 13
Equation = 1
Material = 1
End
Equation 1
Active Solvers(2) = 1 2
End
Material 1
Relative Permittivity = 1
End
Solver 1
Equation = Electrostatics
Procedure = "StatElecSolve" "StatElecSolver"
Variable = Potential
Variable DOFs = 1
Calculate Electric Field = True
Calculate Electric Flux = False
Linear System Solver = Iterative
Linear System Iterative Method = CG
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-8
Linear System Preconditioning = ILU1
End
Solver 2
Equation = Electric Force
Procedure = "ElectricForce" "StatElecForce"
End
Boundary Condition 1
Target Boundaries = 10
Potential = 650
Calculate Electric Force = True
End
Boundary Condition 2
Target Boundaries = 20
Potential = 0
End

View File

@ -0,0 +1,91 @@
import numpy as np
import gmsh
gmsh.initialize()
gmsh.model.add("adaptive_mesh")
gmsh.option.setNumber('General.NumThreads', 4)
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000)
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000)
#gmsh.option.setNumber("Mesh.Adapt.MaxIter",5)
#gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3)
#gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0)
#gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
lc = 0.01
#anodes, plane 1 at -zmax/2
k=-2*np.pi/24.
offset = 6*k + 3*k #-pi/2
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
#cathodes, plane 1 at -zmax/2
kc=2*np.pi/24.
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
#anodes, plane 2 at +zmax/2
offset = offset-3*k
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
#cathodes, plane2 at +zmax/2
offsetc = offsetc-3*kc
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
pa1 = []
pa2 = []
pc1 = []
pc2 = []
wire_radius = 0.254 #mm
anode_wires = []
cathode_wires = []
aw_tags = [(3,i) for i in range(24)]
cw_tags = [(3,i+24) for i in range(24)]
for i,[xa,ya,xc,yc,xa2,ya2,xc2,yc2] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1,xarra_2,yarra_2,xarrc_2,yarrc_2)):
print(i,xa,ya,-178.3,xc,yc,-178.3,xa2,ya2,178.3,xc2,yc2,178.3)
anode_wires.append(gmsh.model.occ.addCylinder(xa,ya,-178.3,(xa2-xa),(ya2-ya),178.3*2,wire_radius,i)) #x,y,z of first face center, dx,dy,dz of the axis, then the wire radius
cathode_wires.append(gmsh.model.occ.addCylinder(xc,yc,-178.3,(xc2-xc),(yc2-yc),178.3*2,wire_radius,i+24)) #cathode tags 24-47, anode 0-23
anasen_barrel = gmsh.model.occ.addCylinder(0,0,-500,0,0,500+605,300,1234) #tag 1234
#anasen_barrel = gmsh.model.occ.addCylinder(0,0,-500,0,0,500+605,300,1234) #tag 1234
gmsh.model.occ.synchronize()
all_wires = aw_tags+cw_tags
gmsh.model.occ.fragment([(3,1234)],all_wires)
gmsh.model.occ.removeAllDuplicates()
gmsh.model.occ.synchronize()
gmsh.option.setNumber("Geometry.Tolerance", 1e-6)
gmsh.option.setNumber("Geometry.OCCFixDegenerated", 1)
gmsh.option.setNumber("Geometry.OCCFixSmallEdges", 1)
gmsh.option.setNumber("Geometry.OCCFixSmallFaces", 1)
wire_surfs = []
for w in anode_wires + cathode_wires:
wire_surfs += [s[1] for s in gmsh.model.getBoundary([(3,w)], oriented=False) if s[0] == 2]
#'''
f1 = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(f1, "FacesList", wire_surfs) # Example curves
f2 = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(f2, "InField", f1)
gmsh.model.mesh.field.setNumber(f2, "SizeMin", 0.05)
gmsh.model.mesh.field.setNumber(f2, "SizeMax", 5.)
gmsh.model.mesh.field.setNumber(f2, "DistMin", 1.)
gmsh.model.mesh.field.setNumber(f2, "DistMax", 20.)
gmsh.model.mesh.field.setAsBackgroundMesh(f2)
#'''
gmsh.option.setNumber("Mesh.Algorithm", 2)
gmsh.option.setNumber("Mesh.Algorithm3D", 1) # For 3D meshes
gmsh.model.mesh.generate(dim=3)
#gmsh.model.mesh.refine()
#gmsh.model.mesh.refine()
#gmsh.model.mesh.refine()
gmsh.fltk.run()
gmsh.finalize()

View File

@ -0,0 +1,87 @@
import numpy as np
import gmsh
gmsh.initialize()
#gmsh.model.add("adaptive_mesh")
gmsh.option.setNumber('General.NumThreads', 4)
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000)
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000)
#gmsh.option.setNumber("Mesh.Adapt.MaxIter",5)
#gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3)
#gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0)
gmsh.option.setNumber("Geometry.Tolerance", 1e-2)
#gmsh.option.setNumber("Mesh.CharacteristicLengthFromCurvature", 0)
lc = 0.04
#anodes, plane 1 at -zmax/2
k=-2*np.pi/24.
offset = 6*k + 3*k #-pi/2
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
#cathodes, plane 1 at -zmax/2
kc=2*np.pi/24.
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
#anodes, plane 2 at +zmax/2
offset = offset-3*k
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
#cathodes, plane2 at +zmax/2
offsetc = offsetc-3*kc
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
wire_radius = 0.254 #mm
anode_wires = []
cathode_wires = []
aw_tags = [(3,i) for i in range(24)]
cw_tags = [(3,i+24) for i in range(24)]
for i,[xa,ya,xc,yc,xa2,ya2,xc2,yc2] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1,xarra_2,yarra_2,xarrc_2,yarrc_2)):
print(i,xa,ya,-178.3,xc,yc,-178.3,xa2,ya2,178.3,xc2,yc2,178.3)
pa1 = gmsh.model.occ.addPoint(xa,ya,-178.3,lc)
pa2 = gmsh.model.occ.addPoint(xa2,ya2,178.3,lc)
pc1 = gmsh.model.occ.addPoint(xc,yc,-178.3,lc)
pc2 = gmsh.model.occ.addPoint(xc2,yc2,178.3,lc)
linea = gmsh.model.occ.addLine(pa1,pa2)
linec = gmsh.model.occ.addLine(pc1,pc2)
anode_wires.append(linea)
cathode_wires.append(linec)
#anasen_barrel = gmsh.model.occ.addCylinder(0,0,-500,0,0,500+605,300,1234) #tag 1234
anasen_barrel = gmsh.model.occ.addCylinder(0,0,-200, 0,0,400, 300) #tag 1234
gmsh.model.occ.synchronize()
gmsh.model.mesh.embed(1,anode_wires+cathode_wires,3,anasen_barrel)
f1 = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(f1,"CurvesList",anode_wires+cathode_wires)
f2 = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(f2,"InField",f1)
gmsh.model.mesh.field.setNumber(f2,"SizeMin",0.08)
gmsh.model.mesh.field.setNumber(f2,"SizeMax",5)
gmsh.model.mesh.field.setNumber(f2,"DistMin",1)
gmsh.model.mesh.field.setNumber(f2,"DistMax",20)
gmsh.model.mesh.field.setAsBackgroundMesh(f2)
gmsh.model.addPhysicalGroup(1, anode_wires, tag=10)
gmsh.model.setPhysicalName(1,10,"anode_wires")
gmsh.model.addPhysicalGroup(1, cathode_wires, tag=20)
gmsh.model.setPhysicalName(1,20,"cathode_wires")
gmsh.option.setNumber("Mesh.Algorithm",6)
gmsh.option.setNumber("Mesh.Algorithm3D", 10) # For 3D meshes
gmsh.model.mesh.generate(dim=3)
gmsh.model.mesh.refine()
#gmsh.model.mesh.refine()
#gmsh.model.mesh.refine()
gmsh.fltk.run()
gmsh.finalize()

BIN
anasen_fem/output.gif Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.8 MiB

62
anasen_fem/paraview_plotter.py Executable file
View File

@ -0,0 +1,62 @@
#!/home/sud/Downloads/ParaView-6.1.0-RC1-MPI-Linux-Python3.12-x86_64/bin/pvbatch
import numpy as np
import sys
from paraview.simple import *
reader = XMLUnstructuredGridReader(FileName=["wires2d/elfield_anasen_t0001.vtu"])
contour_filter = Contour(Input=reader,ContourBy = 'potential')
contour_filter.Isosurfaces = [i for i in np.arange(0,660,650/24.)]
renderView = GetActiveViewOrCreate('RenderView')
renderView.ViewSize = [800,800]
renderView.OrientationAxesVisibility = 0 # Hide axis
renderView.UseColorPaletteForBackground=0
renderView.Background = [0.1, 0.1, 0.1] # Set background to dark gray (RGB 0-1)
renderView.MultiSamples = 8 # 0 disables it, 4-8 is usually sufficient
ResetCamera()
contour_display = Show(contour_filter, renderView)
#colorbar
contour_display_potentialLUT = GetColorTransferFunction('potential', contour_display, separate=True)
contour_display_potentialLUT.ApplyPreset('Cool to Warm', True)
contour_display.SetScalarBarVisibility(renderView, True)
#axesGrid = renderView.AxesGrid
#axesGrid.Visibility = 1
#axesGrid.XTitle = "x (mm)"
#axesGrid.YTitle = "y (mm)"
# 1. Get the active view
view = GetActiveView()
# 2. Define your desired coordinate ranges (x_min, x_max, y_min, y_max, z_min, z_max)
# Example: Look at a box from -10 to 10 in all dimensions
x_min, x_max = -50.0, 50.0
y_min, y_max = -50.0, 50.0
z_min, z_max = -50.0, 50.0
# 3. Calculate Center, Position, and Parallel Scale
center = [(x_min + x_max) / 2.0, (y_min + y_max) / 2.0, (z_min + z_max) / 2.0]
# Position the camera far away along Z to look at the center
position = [center[0], center[1], z_min - 30.0]
# Parallel scale defines how much of the scene is visible.
# It is usually half the height of the viewed area.
view.CameraParallelScale = max((x_max - x_min), (y_max - y_min))/1.6
# 4. Apply settings
view.CenterOfRotation = center
view.CameraPosition = position
view.CameraFocalPoint = center
view.CameraViewUp = [0.0, 1.0, 0.0] # Y-axis is up
# 5. Enable Parallel Projection (optional, often better for exact mapping)
view.CameraParallelProjection = 1
#ResetCamera()
Render()
SaveScreenshot("contour_output.png")

15
anasen_fem/run.py Normal file
View File

@ -0,0 +1,15 @@
import os
#val=-178.3
val=17.83
count=11
while val<178.3+0.1:
print(val)
os.system("python3 wires_gmsh2d_bc.py "+str(val))
os.system("ElmerGrid 14 2 wires2d.msh")
os.system("ElmerSolver wires2d.sif")
os.system("./paraview_plotter.py")
os.system("cp contour_output.png contour_output_z_%02d_%1.4f.png"%(count,val))
val=val+17.83
count = count + 1

1
anasen_fem/scalars.dat Normal file
View File

@ -0,0 +1 @@
6.500000000000E+002

View File

@ -0,0 +1,8 @@
Metadata for SaveScalars file: ./scalars.dat
Elmer version: 26.1
Elmer compilation date: 2026-03-05
Solver input file: wires2d.sif
File started at: 2026/03/11 23:33:58
Variables in columns of matrix:
1: res: potential difference

103
anasen_fem/wires2d.sif Normal file
View File

@ -0,0 +1,103 @@
Check Keywords Warn
Header
Mesh DB "." "wires2d"
End
Simulation
Coordinate System = Cartesian 2D
Simulation Type = Steady State
Steady State Max Iterations = 1
Output File = "elstatics.result"
Post File = "elstatics.ep"
End
Constants
Permittivity Of Vacuum = 8.8542e-12
End
Body 1
Target Bodies(1) = 13
Equation = 1
Material = 1
End
Equation 1
Active Solvers(2) = 1 2
End
Material 1
Relative Permittivity = 1
End
Solver 1
Equation = Electrostatics
Procedure = "StatElecSolve" "StatElecSolver"
Variable = Potential
Variable DOFs = 1
Calculate Electric Field = True
Calculate Electric Flux = False
Linear System Solver = Iterative
Linear System Iterative Method = CG
Linear System Max Iterations = 5000
Linear System Convergence Tolerance = 1.0e-8
Linear System Preconditioning = ILU1
Calculate Vectors = Logical True
End
Solver 2
Equation = "Electric Field"
Procedure = "FluxSolver" "FluxSolver"
! Calculate from the potential
Target Variable = "Potential"
! Name of the output vector field in VTU
Flux Variable = String "Electric Field"
! Use 2D components (x, y)
Flux Coefficient = String "Permittivity"
Calculate Vectors = Logical True
End
Solver 3
Equation = Result Output
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Output File Name = elfield_anasen ! Sets prefix for output files
Output Format = Vtu
! Optional: Select specific variables to save
Scalar Field 1 = Potential
Vector Field 1 = Electric Field
End
Solver 4
Exec Solver = After All
Equation = SaveScalars
Procedure = "SaveData" "SaveScalars"
Filename = "scalars.dat"
End
Boundary Condition 1
Target Boundaries = 10
Potential = 650
Calculate Electric Force = True
End
Boundary Condition 2
Target Boundaries = 20
Potential = 0
End
!Boundary Condition 2
! Target Boundaries = 30
! Potential = 0
!End

View File

@ -0,0 +1,126 @@
import numpy as np
import gmsh,sys
gmsh.initialize()
gmsh.model.add("adaptive_mesh")
gmsh.option.setNumber('General.NumThreads', 4)
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfElements", 200000)
#gmsh.option.setNumber("Mesh.Adapt.MaxNumberOfNodes", 200000)
#gmsh.option.setNumber("Mesh.Adapt.MaxIter",5)
#gmsh.option.setNumber("Mesh.MeshSizeMin", 5e-3)
#gmsh.option.setNumber("Mesh.MeshSizeMax", 10.0)
gmsh.option.setNumber("Geometry.Tolerance", 4e-2)
#gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
lc = 0.04
#z_loc = -178.3
if len(sys.argv) < 2:
print("Usage: python3 wires_gmsh2d_bc.py <z_locus in mm>")
quit()
z_loc = float(sys.argv[1])
k=(2*np.pi/24.)
#anodes, plane 1 at -zmax/2
k=-2*np.pi/24.
offset = 6*k + 3*k #-pi/2
xarra_1 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_1 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
#cathodes, plane 1 at -zmax/2
kc=2*np.pi/24.
offsetc = -4*kc + 2*kc - np.pi/24 #-pi/4
xarrc_1 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_1 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
#anodes, plane 2 at +zmax/2
offset = offset-3*k
xarra_2 = np.array([37*np.cos(k*i+offset) for i in np.arange(0,24)])
yarra_2 = np.array([37*np.sin(k*i+offset) for i in np.arange(0,24)])
#cathodes, plane2 at +zmax/2
offsetc = offsetc-3*kc
xarrc_2 = np.array([42*np.cos(kc*i+offsetc) for i in np.arange(0,24)])
yarrc_2 = np.array([42*np.sin(kc*i+offsetc) for i in np.arange(0,24)])
direction_anodes_x = xarra_2 - xarra_1
direction_anodes_y = yarra_2 - yarra_1
direction_cathodes_x = xarrc_2 - xarrc_1
direction_cathodes_y = yarrc_2 - yarrc_1
t = (z_loc+178.3)/(2*178.3) #z=-178.3 is 0, z=+178.3 is 1
xloc_a = xarra_1 + t*direction_anodes_x
yloc_a = yarra_1 + t*direction_anodes_y
xloc_c = xarrc_1 + t*direction_cathodes_x
yloc_c = yarrc_1 + t*direction_cathodes_y
wire_radius = 0.254 #mm
anode_wires = []
cathode_wires = []
aw_tags = [(3,i) for i in range(24)]
cw_tags = [(3,i+24) for i in range(24)]
#for i,[xa,ya,xc,yc] in enumerate(zip(xarra_1,yarra_1,xarrc_1,yarrc_1)):
for i,[xa,ya,xc,yc] in enumerate(zip(xloc_a,yloc_a,xloc_c,yloc_c)):
print(i,xa,ya,-178.3,xc,yc,-178.3)
adisk = gmsh.model.occ.addDisk(xa,ya,0,wire_radius,wire_radius)
cdisk = gmsh.model.occ.addDisk(xc,yc,0,wire_radius,wire_radius)
anode_wires.append(adisk)
cathode_wires.append(cdisk)
anasen_barrel = gmsh.model.occ.addDisk(0,0,0,500,500)
#gmsh.model.occ.synchronize()
#gmsh.model.mesh.embed(1,anode_wires+cathode_wires,2,anasen_barrel)
gmsh.option.setNumber("Geometry.Tolerance", 1e-6)
gmsh.option.setNumber("Geometry.OCCFixDegenerated", 1)
gmsh.model.occ.synchronize()
awire_surfs = []
for w in anode_wires:
awire_surfs += [s[1] for s in gmsh.model.getBoundary([(2,w)], oriented=False) if s[0] == 1]
cwire_surfs = []
for w in cathode_wires:
cwire_surfs += [s[1] for s in gmsh.model.getBoundary([(2,w)], oriented=False) if s[0] == 1]
gmsh.model.mesh.embed(1,cwire_surfs+awire_surfs,2,anasen_barrel)
for s in gmsh.model.getBoundary([(2,w)],oriented=False):
if s[0] == 1:
anasen_bdry=s[1]
f1 = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(f1,"CurvesList",cwire_surfs+awire_surfs)
f2 = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(f2,"InField",f1)
gmsh.model.mesh.field.setNumber(f2,"SizeMin",0.1)
gmsh.model.mesh.field.setNumber(f2,"SizeMax",5)
gmsh.model.mesh.field.setNumber(f2,"DistMin",1)
gmsh.model.mesh.field.setNumber(f2,"DistMax",20)
gmsh.model.mesh.field.setAsBackgroundMesh(f2)
gmsh.model.addPhysicalGroup(1, awire_surfs, tag=10)
gmsh.model.setPhysicalName(1,10,"anode_wires")
gmsh.model.addPhysicalGroup(1, cwire_surfs, tag=20)
gmsh.model.setPhysicalName(1,20,"cathode_wires")
#gmsh.model.addPhysicalGroup(1, [anasen_bdry], tag=30)
#gmsh.model.setPhysicalName(1,30,"barrel_boundary")
gmsh.model.addPhysicalGroup(2,[anasen_barrel],tag=13)
gmsh.model.setPhysicalName(1,13,"gas")
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(dim=2)
gmsh.model.mesh.refine()
gmsh.write("wires2d.msh")
#gmsh.fltk.run()
gmsh.finalize()