Skip to content
  • Search
Skins
  • Light
  • Cerulean
  • Cosmo
  • Flatly
  • Journal
  • Litera
  • Lumen
  • Lux
  • Materia
  • Minty
  • Morph
  • Pulse
  • Sandstone
  • Simplex
  • Sketchy
  • Spacelab
  • United
  • Yeti
  • Zephyr
  • Dark
  • Cyborg
  • Darkly
  • Quartz
  • Slate
  • Solar
  • Superhero
  • Vapor

  • Default (No Skin)
  • No Skin
Collapse

ZMT zurich med tech

  1. Home
  2. Sim4Life
  3. Python API
  4. Calculating flux from simulation results - Improving Interpolation

Calculating flux from simulation results - Improving Interpolation

Scheduled Pinned Locked Moved Python API
interpolationpythoncodeexamplesanalysis
2 Posts 2 Posters 890 Views
  • Oldest to Newest
  • Newest to Oldest
  • Most Votes
Reply
  • Reply as topic
Log in to reply
This topic has been deleted. Only users with topic management privileges can see it.
  • J Offline
    J Offline
    Jexyll_S
    wrote on last edited by
    #1

    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

    1 Reply Last reply
    0
    • SylvainS Offline
      SylvainS Offline
      Sylvain
      ZMT
      wrote on last edited by
      #2

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

      1 Reply Last reply
      0
      Reply
      • Reply as topic
      Log in to reply
      • Oldest to Newest
      • Newest to Oldest
      • Most Votes


      • Login

      • Don't have an account? Register

      • Login or register to search.
      • First post
        Last post
      0
      • Search