# 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 B-field from the simulation results, then add all the B-fields 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
simulation = document.AllSimulations["Self Inductance TMS Device"]
simulation_extractor = simulation.Results()

em_sensor_extractor = simulation_extractor["Overall Field"]
em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All"

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 B-field. 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])*1e-3

# Find index closest to x-coordinate 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
simulation = document.AllSimulations["Self Inductance TMS Device"]
simulation_extractor = simulation.Results()

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.07106769085e-05, units.Meters
model_to_grid_filter.UpdateAttributes()

em_sensor_extractor = simulation_extractor["Overall Field"]
em_sensor_extractor.FrequencySettings.ExtractedFrequency = u"All"

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()

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()

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).