diff --git a/anasen_fem/clean.sh b/anasen_fem/clean.sh
new file mode 100644
index 0000000..f0c8051
--- /dev/null
+++ b/anasen_fem/clean.sh
@@ -0,0 +1,2 @@
+rm wires2d.msh
+rm wires2d/*
diff --git a/anasen_fem/junk/wires.py b/anasen_fem/junk/wires.py
new file mode 100644
index 0000000..4ba5f5e
--- /dev/null
+++ b/anasen_fem/junk/wires.py
@@ -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()
+
diff --git a/anasen_fem/junk/wires2d_test.sif b/anasen_fem/junk/wires2d_test.sif
new file mode 100644
index 0000000..5e2788f
--- /dev/null
+++ b/anasen_fem/junk/wires2d_test.sif
@@ -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
diff --git a/anasen_fem/junk/wires_gmsh.py b/anasen_fem/junk/wires_gmsh.py
new file mode 100644
index 0000000..10317c4
--- /dev/null
+++ b/anasen_fem/junk/wires_gmsh.py
@@ -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()
diff --git a/anasen_fem/junk/wires_gmsh_bc.py b/anasen_fem/junk/wires_gmsh_bc.py
new file mode 100644
index 0000000..a631070
--- /dev/null
+++ b/anasen_fem/junk/wires_gmsh_bc.py
@@ -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()
diff --git a/anasen_fem/output.gif b/anasen_fem/output.gif
new file mode 100644
index 0000000..abb4b2e
Binary files /dev/null and b/anasen_fem/output.gif differ
diff --git a/anasen_fem/paraview_plotter.py b/anasen_fem/paraview_plotter.py
new file mode 100755
index 0000000..1959251
--- /dev/null
+++ b/anasen_fem/paraview_plotter.py
@@ -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")
diff --git a/anasen_fem/run.py b/anasen_fem/run.py
new file mode 100644
index 0000000..2b4309c
--- /dev/null
+++ b/anasen_fem/run.py
@@ -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
+
diff --git a/anasen_fem/scalars.dat b/anasen_fem/scalars.dat
new file mode 100644
index 0000000..14ef6f3
--- /dev/null
+++ b/anasen_fem/scalars.dat
@@ -0,0 +1 @@
+   6.500000000000E+002
diff --git a/anasen_fem/scalars.dat.names b/anasen_fem/scalars.dat.names
new file mode 100644
index 0000000..b9970f2
--- /dev/null
+++ b/anasen_fem/scalars.dat.names
@@ -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
diff --git a/anasen_fem/wires2d.sif b/anasen_fem/wires2d.sif
new file mode 100644
index 0000000..6bb38f6
--- /dev/null
+++ b/anasen_fem/wires2d.sif
@@ -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
diff --git a/anasen_fem/wires_gmsh2d_bc.py b/anasen_fem/wires_gmsh2d_bc.py
new file mode 100644
index 0000000..0023838
--- /dev/null
+++ b/anasen_fem/wires_gmsh2d_bc.py
@@ -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()