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
	# 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 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
	radius = 1*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
	# 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.07106769085e-05, 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).