Calculating flux from simulation results  Improving Interpolation

Hi everyone,
I am trying to calculate the magnetic flux through a small, circular aperture using python scripts. My current process is to extract the values for the Bfield from the simulation results, then add all the Bfields in a radius and multiply by the aperture area. The code for this is below.import numpy import s4l_v1.analysis as analysis import s4l_v1.document as document import s4l_v1.model as model import s4l_v1.units as units from s4l_v1 import ReleaseVersion from s4l_v1 import Unit # Creating the analysis pipeline # Adding a new SimulationExtractor simulation = document.AllSimulations["Self Inductance TMS Device"] simulation_extractor = simulation.Results() # Adding a new EmSensorExtractor em_sensor_extractor = simulation_extractor["Overall Field"] em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All" document.AllAlgorithms.Add(em_sensor_extractor) B_field = em_sensor_extractor.Outputs["B(x,y,z,f0)"] B_field.Update() out = B_field.Data.Field(0) # Extract and correct dimensionality dims = np.array(B_field.Data.Grid.Dimensions) if np.prod(dims) != len(out): dims = dims  1 # get x component of Bfield. Reshape to correspond to volume # Creates an array of [x,y,z] values out_Bx = out[:,0].reshape(dims, order = 'F') # Get cooordinates for cell nodes (not cell centres) xaxis_nodes = field.Data.Grid.XAxis yaxis_nodes = field.Data.Grid.YAxis zaxis_nodes = field.Data.Grid.ZAxis # Extrapolate coordinate points xaxis = (xaxis_nodes[:1] + xaxis_nodes[1:])/2 yaxis = (yaxis_nodes[:1] + yaxis_nodes[1:])/2 zaxis = (zaxis_nodes[:1] + zaxis_nodes[1:])/2 # Define circular apperture centre = np.array([0, 0, 15])*1e3 radius = 1*1e3 # Find index closest to xcoordinate of centre centre_idx_x = (np.abs(xaxis  centre[0])).argmin() # Calculate flux through apperture flux = 0.0 for y_idx, ycoord in enumerate(yaxis): for z_idx, zcoord in enumerate(zaxis): if ((ycoord  centre[1])**2 + (zcoord  centre[2])**2) <= radius**2: # print("y and z coords are {} and {}".format(ycoord, zcoord)) # print("Value being added to flux is %f" % str(out_Bx[centre_idx_x, ycoord, zcoord])) flux = flux + out_Bx[centre_idx_x, y_idx, z_idx] flux = flux * np.pi * radius**2 print("Final flux is {}".format(flux))
In addition, I also attempted using the flux evaluator tool. However, when I create circle entities, I am unable to extract any surfaces. Instead, I had to use an infinitely thin cylinder and apply suitable transforms to position it appropriately.
import numpy as np import s4l_v1.analysis as analysis import s4l_v1.document as document import s4l_v1.model as model import s4l_v1.units as units from s4l_v1 import ReleaseVersion from s4l_v1 import Unit from s4l_v1 import Vec3 # Define the version to use for default values ReleaseVersion.set_active(ReleaseVersion.version5_2) cylinder2 = model.CreateSolidCylinder(Vec3(0,0,0),Vec3(0,0,0),0.5) rotation = model.Rotation(Vec3(0,1,0), np.pi/2) translation = model.Translation(Vec3(0,0,15)) cylinder2.ApplyTransform(rotation) cylinder2.ApplyTransform(translation) cylinder2.Name = "Cylinder 2" # Creating the analysis pipeline # Adding a new SimulationExtractor simulation = document.AllSimulations["Self Inductance TMS Device"] simulation_extractor = simulation.Results() # Adding a new ModelToGridFilter inputs = [] model_to_grid_filter = analysis.core.ModelToGridFilter(inputs=inputs) model_to_grid_filter.Name = "Cylinder 2" model_to_grid_filter.Entity = model.AllEntities()["Cylinder 2"] model_to_grid_filter.MaximumEdgeLength = 7.07106769085e05, units.Meters model_to_grid_filter.UpdateAttributes() document.AllAlgorithms.Add(model_to_grid_filter) # Adding a new EmSensorExtractor em_sensor_extractor = simulation_extractor["Overall Field"] em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All" document.AllAlgorithms.Add(em_sensor_extractor) # Adding a new FieldInterpolationFilter inputs = [em_sensor_extractor.Outputs["B(x,y,z,f0)"], model_to_grid_filter.Outputs["Surface"]] field_interpolation_filter = analysis.core.FieldInterpolationFilter(inputs=inputs) field_interpolation_filter.UpdateAttributes() document.AllAlgorithms.Add(field_interpolation_filter) # Adding a new SurfaceFieldFluxEvaluator inputs = [field_interpolation_filter.Outputs["B(x,y,z,f0)"]] surface_field_flux_evaluator = analysis.core.SurfaceFieldFluxEvaluator(inputs=inputs) surface_field_flux_evaluator.UpdateAttributes() surface_field_flux_evaluator.Update() document.AllAlgorithms.Add(surface_field_flux_evaluator) flux_data = surface_field_flux_evaluator.Outputs["Total Flux(B(x,y,z,f0))"].Data mag_flux = flux_data.GetComponent(0)[0] print("magnetic flux is {}".format(mag_flux))
I prefer the second method, since it interpolates values, providing a more accurate measure of the flux. I wanted to post this code to provide some more references which new users of sim4life can draw from, and also to get feedback on my coding techniques.
Cheers,
Jexyll 
This looks quite good, thank you for providing your scripts.
Note that, for the second method, you could create your cylinder directly at the desired position:cylinder2 = model.CreateSolidCylinder(Vec3(0,1,15),Vec3(0,1,15),0.5)
There is also the
XCoredModeling.CoverWires()
function that allows you to make a surface (i.e. face) out of a circle entity.Last, but not least, note that the line
model_to_grid_filter.MaximumEdgeLength
is ultimately what determines the tradeoff between accuracy and computational cost of the interpolation (since it defines the resolution of the triangulated mesh on which the interpolation is done).