Overview

• Introduction to the API

  This document provides a comprehensive guide to using our Python API for Intact.Simulation. The API allows users to define geometries, material properties, components, assemblies, and boundary conditions. Users can specify units, create simulation scenarios, execute the simulation, and query the results. The simulations can be run with different solver types and the results can be written to a file for visualization. The document also provides Python code examples and snippets for each step of the simulation process.

• Key features and capabilities

  • Physics Support

    • Stress

    • Modal

    • Thermal

    • Linear Buckling

  • Geometry Support

    • Brep: Surface mesh

    • VDB volume

  • Level set topology optimization

Getting Started

• A simple example of running a structural simulation. Any geometry files used in these examples are available in the examples/ folder.

  • Import Intact Module

  from PyIntact import *
  

  • Create Geometry

  beam_geometry = MeshModel("beam.stl")
  beam_geometry.instance_id = "simple_cantilever"
  beam_geometry.refine(0.02) # refine stl for smoother visualization
  

  • Create Material (Properties are in MKS unit system)

  material = IsotropicMaterialDescriptor() # steel
  material.density = 7800.0 # kg/m^3
  material.poisson_ratio = 0.3
  material.youngs_modulus = 2.1e11 # Pa
  

  • Create Restraints and Loads

  # Create a fixed boundary condition at one end of the beam
  fixed_boundary = FixedBoundaryDescriptor()
  fixed_boundary.boundary = MeshModel("restraint.stl")
  
  # Create a pressure load at the free end of the beam
  pressure_load = PressureForceDescriptor()
  pressure_load.units = UnitSystem.MeterKilogramSecond
  pressure_load.magnitude = 1000  # Applying pressure in Pascals
  pressure_load.boundary = MeshModel("load.stl")
  

  • Setup simulation and solve

  # Setup the simulation scenario descriptor
  scenario = LinearElasticScenarioDescriptor()
  scenario.materials = {"Steel": material}
  scenario.metadata.resolution = 10000
  scenario.metadata.units = UnitSystem.MeterKilogramSecond
  scenario.boundary_conditions = [fixed_boundary, pressure_load]
  
  # Associate the material with the model
  component = MaterialDomain(beam_geometry, "Steel", scenario)
  assembly = [component]
  
  # Initialize and run the simulator
  simulator = StressSimulator(assembly, scenario)
  simulator.solve()
  

  • Query result and create output file to visualize

  # Query results for stress distribution
  results = QueryResult(assembly)
  stress_query = FieldQuery(Field.VonMisesStress)
  simulator.sample(stress_query, results)
  
  # Write results to a file (unit system is saved as a metadata in the vtu file)
  results.writeVTK("cantilever_beam_results.vtu", UnitSystem.MeterKilogramSecond)
  

How To

Creating Simulation Input

Units

The default Unit System is MKS or MeterKilogramSecond. We allow the following customizations:

• Specify custom unit for geometry through the Scenario unit (default MKS)

• Specify custom unit for loads (default MKS)

• Specify custom unit for materials (default MKS)

Unit System Keywords	Mass	Force	Stress	Length
MeterKilogramSecond	kg	N	Pa	m
CentimeterGramSecond	g	dyne	dyne/cm²	cm
MillimeterMegagramSecond	Mg	N	MPa	mm
FootPoundSecond	slug	lbf	lbf/ft²	ft
InchPoundSecond	lbf s²/in	lbf	psi	in

Geometry

Geometry is specified through a Model object

• MeshModel from surface mesh can be defined in two ways.

  • filename to define mesh from a file (STL, PLY)

    • instance_id a unique name

    mesh_geometry = MeshModel(filename="beam.stl")
    mesh_geometry.instance_id = "beam"
    

  • Constructing mesh face-by-face

    # define an empty mesh
    mesh_geometry = MeshModel()
    mesh_geometry.instance_id = "fixed_boundary"
    
    # add vertices to the mesh and get the vertex id as output
    v0 = mesh_geometry.addVertex([0.0, 0.0, 0.0])
    v1 = mesh_geometry.addVertex([1.0, 0.0, 0.0])
    v2 = mesh_geometry.addVertex([0.0, 1.0, 0.0])
    
    # add faces defined by the vertex id
    mesh_geometry.addFacet(v0, v1, v2)
    

  • Mesh refinement (usually used for finer interpolation of results during visualization). Uses MeshModel.refine with refinement_level as an input which is the largest allowed triangle edge as a ratio of the bounding box diagonal.

    # mesh geometry created above is further refined
    mesh_geometry.refine(refinement_level=0.02)
    

• VDB Model

  • Using a .vdb file

    vdb_geometry= VDBModel("box.vdb")
    

  • In memory using a VDB object

    import pyopenvdb as vdb
    
    #create a vdb double grid object
    cube = vdb.DoubleGrid()
    cube.fill(min=(100, 100, 100), max=(199, 199, 199), value=1.0)
    
    #create an Intact Model from VDB object
    vdb_geometry = PyIntact.VDBModel(cube)
    
    #create a tesselation of the vdb to get a MeshModel and its list of vertices for sampling
    mesh = vdb_geometry.tesselate() # mesh is a MeshModel()
    num_vertices = mesh.vertexCount()
    print(num_vertices)
    

Material Properties

• Structural

  • IsotropicMaterialDescriptor to create an Isotropic material. (Unit is MeterKilogramSecond)

    • density is the material density

    • youngs_modulus is the elastic modulus

    • poisson_ratio is the Poisson ratio

    steel_structural = IsotropicMaterialDescriptor()
    
    # Set the density to 7845 kg/m^3, modulus to 200 GPa and poison ratio to 0.29
    steel_structural.density = 7845
    steel_structural.youngs_modulus = 200e9
    steel_structural.poisson_ratio = 0.29
    

  • OrthotropicMaterialDescriptor to create an Orthotropic material. Orthotropic materials are often used to represent composites and wood, where the material properties along one axis are significantly different from the properties along the other axes. (Unit is MeterKilogramSecond)

  • density is the material density

  • Ex is the elastic modulus in the material’s x-direction

  • Ey is the elastic modulus in the material’s y-direction

  • Ez is the elastic modulus in the material’s z-direction

  • Gxy is the shear modulus in the xy plane

  • Gxz is the shear modulus in the xz plane

  • Gyz is the shear modulus in the yz plane

  • vxy is the poisson ratio in the xy plane

  • vxz is the poisson ratio in the xx plane

  • vyz is the poisson ratio in the yz plane

  • transform is the optional material transform to arbitrarily rotate the axes along which the material properties are defined. Specified as a list of 9 floating point numbers corresponding to a 3x3 matrix. Note the first 3 floating point numbers corresponds to the first row.

  # Create an orthotropic material (Red Pine)
  material = OrthotropicMaterialDescriptor()
  material.density = 460.0             # kg/m³
  material.Ex = 11.2e9                 # Pa
  material.Ey = 492.8e6                # Pa
  material.Ez = 985.6e6                # Pa
  material.Gxy = 907.2e6               # Pa
  material.Gxz = 1.08e9                # Pa
  material.Gyz = 123.2e6               # Pa
  material.vxy = 0.315
  material.vxz = 0.347
  material.vyz = 0.308
  material.transform = [1.0, 0.0, 0.0,
                      0.0, 1.0, 0.0,
                      0.0, 0.0, 1.0] # optional in this case
  

• Thermal

  • ThermalMaterialDescriptor to create a thermal material. (Default unit is MeterKilogramSecond)

    • density is the material density

    • conductivity is the conductivity of the material

    • specific_heat is the specific heat of the material

    thermal_material = ThermalMaterialDescriptor()
    
    # Set the density, conductivity, and specific heat of the material
    thermal_material.density = 2700.0 # kg/m^3
    thermal_material.conductivity = 170  # W/(m·K)
    thermal_material.specific_heat = 500  # J/(kg·K)
    

Component and Assembly

• A Component is created as a MaterialDomain and takes three inputs

  • A Model object (MeshModel or VDBModel specifically)

  • A material name

  • A ScenarioDescriptor object

  bar1_structural = MaterialDomain(mesh_geometry, "Steel", structural_scenario)
  bar2_structural = MaterialDomain(vdb_geometry, "Steel", structural_scenario)
  

• Assembly is a list of components that are bonded together

  structural_assembly = [bar1_structural, bar2_structural]
  

Boundary Conditions

• Restraints

  • Fixed Boundary A “Fixed Boundary” restraint fixes the selected geometry in all directions.

    The fixed boundary only has one input:

    • the boundary surface to be fixed

    fixed = FixedVectorDescriptor()
    
    # Set the geometry "restraint.ply" to be fixed
    fixed.boundary = MeshModel("restraint.ply")
    

  • Fixed Vector A “Fixed Vector” restraint allows for each direction to be optionally set to a specified displacement value, 0 being fixed and ‘None’ being un-restrained. Note that a structural problem must have all three directions restrained somewhere to be valid. Also, note that this boundary condition doesn’t take custom units and it’s always in the units of the scenario metadata.

    A fixed vector requires four inputs:

    • the boundary surface to be fixed

    • the x_value to set the displacement for the x-axis direction (optional)

    • the y_value to set the displacement for the y-axis direction (optional)

    • the z_value to set the displacement for the z-axis direction (optional)

    fixed_vector = FixedVectorDescriptor()
    
    # Set the geometry "restraint.stl" to be restrained
    fixed_vector.boundary = MeshModel("restraint.stl")
    
    # Set the displacements in each direction
    fixed_vector.x_value = 0.2  # X-direction displacement of 0.2 m
    fixed_vector.y_value = None # un-restrained in the Y-direction at this surface
    fixed_vector.z_value = 0.0  # fixed Z-direction
    

  • Sliding Restraint A “Sliding Restraint” allows for motion tangential to a specified surface, but fixes motion normal to the specified surface.

    A sliding restraint only has one input:

    • the boundary surface for the sliding restraint condition

    sliding_boundary = SlidingBoundaryDescriptor()
    
    # Set the geometry "restraint.stl" for the sliding restraint condition
    sliding_boundary.boundary = MeshModel("restraint.stl")
    

• Structural Loads (TractionDescriptor)

  • Vector Force

    “Vector Force” load is a surface load applied to a face in a specified direction. An example of this load is pressing on the top of a book to push it across a table.

    A vector load requires four inputs:

    • the boundary surfaces where the load is applied

    • the direction vector of the force

    • the UnitSystem that applies to the magnitude

    • the magnitude of the force.

    load = VectorForceDescriptor()
    
    # Set the geometry "load.stl" the vector load is applied to
    load.boundary = MeshModel("load.stl")
    
    # Set the vector load direction to be in the -Z direction with a magnitude of 100 lbf
    load.direction = [0, 0, -1]
    load.units = UnitSystem.InchPoundSecond
    load.magnitude = 100 # lbf
    

  • Torque

    “Torque” load is a surface load that applies a twisting force around an axis. The direction of the torque is determined using the right-hand rule: using your right hand, point your thumb in the direction of the axis. A positive torque value applies a torque acting in the direction the fingers of your right hand would wrap around the axis. The torque load is applied among the load faces with a distribution that varies linearly from zero at the axis.

    A Torque load requires five inputs:

    • the boundary surfaces where the load is applied

    • the axis of rotation about which the torque acts

    • origin is the starting point of the axis of rotation

    • the UnitSystem that applies to the magnitude

    • magnitude of the torque.

    torque_load = TorqueForceDescriptor()
    
    # Set the geometry "load.stl" the torque load is applied to
    torque_load.boundary = MeshModel("load.stl")
    
    # Set the axis of the torque axis
    torque_load.origin = (10, 1, 1)
    torque_load.axis = [1, 0, 0] # Set the torque to be about the +X (right-hand rule)
    torque_load.units = UnitSystem.MeterKilogramSecond
    torque_load.magnitude = 10 # Set the torque magnitude to 10 N*m
    

  • Pressure

    A “Pressure” load is a surface load specified in terms of force per unit area. Positive pressures ‘push’ into the surface, and negative pressures ‘pull’.

    A Pressure Load requires three inputs:

    • the boundary surfaces where the load is applied

    • the UnitSystem that applies to the magnitude

    • the magnitude of the pressure.

    pressure_load = PressureForceDescriptor()
    
    # Set the geometry "load.stl" the pressure load is applied to
    pressure_load.boundary = MeshModel("load.stl")
    
    # Set the pressure magnitude to 10 Pa
    pressure_load.units = UnitSystem.MeterKilogramSecond
    pressure_load.magnitude = 10
    

  • Bearing Force

    A “Bearing Force” is a surface load applied to a (typically) cylindrical face to approximate the effects of a shaft pressing against the side of a hole. The applied force gets converted to a varying pressure distribution on the portion of the face experiencing compressive pressure. The pressure distribution is computed automatically to achieve the specified overall bearing force.

    A Bearing Force requires four inputs:

    • the boundary surfaces where the load is applied

    • the direction vector of the bearing force

    • the UnitSystem that applies to the magnitude

    • the magnitude of the force

    bearing_load = BearingForceDescriptor()
    
    # Set the geometry "load.stl" the bearing load is applied to
    bearing_load.boundary = MeshModel("load.stl")
    
    # Set the loading direction to be in the -Z
    bearing_load.direction = [0, 0, -1]
    
    # Set the magnitude of the load to 100 N
    bearing_load.units = UnitSystem.MeterKilogramSecond
    bearing_load.magnitude = 100
    

    • Hydrostatic Load

    A “Hydrostatic” load is a spatially varying pressure on the surfaces submerged in a fluid due to the weight of that fluid. The pressure at any point depends on the height and density of the fluid, increasing from zero at the fluid surface to a maximum at the deepest point.

    A Hydrostatic load requires four inputs:

    • the boundary surfaces where the load is applied

    • the height of the fluid surface

    • the density of the fluid

    • the UnitSystem that applies to the fluid height and density

    # Create a hydrostatic load on a face of the box
    hydrostatic_load = HydrostaticForceDescriptor()
    
    # Set the geometry "load.stl" the hydrostatic load is applied to
    hydrostatic_load.boundary = MeshModel("load.stl")
    
    # Set the fluid height to 78.0 cm with a fluid density of 1.00 g/cm^3
    hydrostatic_load.units = UnitSystem.CentimeterGramSecond
    hydrostatic_load.height = 78
    hydrostatic_load.density = 1.0
    

• Thermal Loads

  • Fixed Boundary (Fixed Temperature)

    A “Fixed Boundary” load fixes the selected geometry to a specified temperature when the value is specified and non-zero.

    The fixed boundary has two inputs:

    • the boundary surface to be fixed

    • the value of the temperature to fix at the boundary surface

    fixed = FixedBoundaryDescriptor()
    
    # Set the geometry "fixed_temp.ply" to set a fixed temperaure of 320 K
    fixed_boundary.boundary = MeshModel("fixed_temp.ply")
    fixed_boundary.value = 320  # Kelvin
    

  • Convection

    A “Convection” load specifies the transfer of heat from a surrounding medium.

    A thermal convection boundary condition requires three inputs:

    • the boundary surface(s) where the convection is applied.

    • the heat transfer coefficient

    • the environment_temperature of the surrounding medium

    • the units (default MKS)

    # Create an instance of ConvectionDescriptor
    convection = ConvectionDescriptor()
    convection.units = UnitSystem.MeterKilogramSecond
    
    # Set the geometry "face.ply" the convection is applied to
    convection.boundary = MeshModel("face.ply")
    
    # Set the heat transfer coefficient to 25 W/m^2K and environment temperature
    convection.coefficient = 25  # W/m^2K
    convection.environment_temperature = 300  # Kelvin
    

  • Surface Flux

    Surface “Thermal or Heat Flux” specifies the heat flow per unit of surface area.

    A surface thermal flux requires two inputs:

    • the boundary surface(s) where the flux is applied.

    • the magnitude of the heat flux

    • the units (default MKS)

    # Create an instance of ConstantFluxDescriptor
    flux = ConstantFluxDescriptor()
    flux.units = UnitSystem.MeterKilogramSecond
    
    # Set the geometry "face.ply" the constant flux is applied to
    flux.boundary = MeshModel("face.ply")
    
    # Set the flux magnitude to 500 W/m^2
    flux.magnitude = 500  # W/m^2
    

• Body Loads/Internal Conditions

  Add body load to the scenario as shown below. (see Scenario Setup section for more details)

  scenario.internal_conditions = [rotational_load, gravity_load]
  

  • Body Load (Linear Acceleration, Gravity)

    Body loads comprise forces that are distributed over a solid volume. They are specified by entering the components of a linear acceleration vector field in which the body is immersed. The material in the body will tend to be pulled in the direction of the acceleration vector.

    The inputs to the body load are:

    • the direction vector of the acceleration field

    • the UnitSystem that applies to the magnitude

    • the magnitude of acceleration

    # Example for a "gravity load"
    body_load = BodyLoadDescriptor()
    
    # Set the direction vector to be downward (-Z)
    body_load.direction = [0, 0, -1]
    
    # Set the magnitude to 9.80655 m/s^2
    body_load.units = UnitSystem.MeterKilogramSecond
    body_load.magnitude = 9.80665
    

  • Rotational Load

    Rotational body loads simulate the effect of a body rotating around an axis. Two contributions are considered in a rotational body load: angular velocity and angular acceleration. The angular velocity term simulates the centrifugal effects that tend to throw a body’s material away from the axis of rotation. The angular acceleration term simulates the effect of a rotational acceleration field around the axis of rotation. A positive angular acceleration tends to drag the body’s material in the positive rotational direction according to the right-hand rule.

    A rotational body load has 4 inputs:

    • origin point for the axis of rotation

    • vector defining the axis of rotation

    • angular_velocity (in radians/sec)

    • angular_acceleration (in radians/sec²)

    rotational_load = RotationalLoadDescriptor()
    
    # Set the origin at the coordinate system origin
    rotational_load.origin = (0, 0, 0)
    
    # Set the axis of rotation about the y-axis
    rotational_load.axis = [0, 1, 0]
    
    # Set angular velocity to 10 rad/s and angular acceleration to 0.5 rad/s^2
    rotational_load.angular_velocity = 10
    rotational_load.angular_acceleration = 0.5
    

• Thermal Loads

  • Constant Heat

    A “Constant Heat” or body heat flux load applies uniform heat generation over a specified volume.

    Constant heat flux has 2 inputs

    • the instance_id of the components which are producing heat flux

    • the magnitude of the body heat flux

    • the units (default MKS)

    constant_heat = ConstantHeatDescriptor()
    constant_heat.units = UnitSystem.MeterKilogramSecond
    
    # Set the heat generation to -200,000 W/m^3 for a beam component
    constant_heat.instance_id = "beam"
    constant_heat.magnitude = -200000.0  # W/m^3
    

Simulation Scenario Setup and Solution

• Scenario Setup using ScenarioDescriptor

  • Linear Elasticity scenario is created using LinearElasticScenarioDescriptor. It takes the following inputs:

    • boundary_conditions, a list of boundary conditions

    • internal_conditions, a list of internal conditions (body loads)

  • Scenario Metadata, metadata, a set of solver parameters to control the accuracy and speed of the simulation:

    • resolution is the target number of finite elements. An iterative process determines a cell_size that achieves approximately the specified number of elements. You can directly specify the cell_size instead of resolution. Note that decreasing cell_size can quickly result in large numbers of finite elements and long solve times. resolution is recommended for most cases.

    • units sets the unit of the scenario (geometry and results will be in this unit) for example, a geometry in ‘MKS’ that is 6 m long would be 6 mm long when set to ‘MMS’.

    • basis_order is the type of finite element used. The default is 1, which would be linear elements. 2 is for quadratic elements.

    • solver_type is the solver type used in the simulation with the following options:

      • MKL_PardisoLDLT (default) is the direct solver and is typically faster for lower resolution (< 200K cells on 32 GB memory), but gets slower and consumes more memory at higher resolutions

      • AMGCL_amg_rigid_body is the iterative solver which is typically faster at high resolution.

    scenario = LinearElasticScenarioDescriptor()
    scenario.materials = {"Aluminum 6061-T6": material}
    scenario.boundary_conditions = [fixed_boundary, vector_load, torque_load]
    scenario.internal_conditions = [body_load]
    
    scenario.metadata.resolution = 10000
    scenario.metadata.units = UnitSystem.MeterKilogramSecond
    
    # Optional settings
    scenario.metadata.basis_order = 2 # 2 = quadratic elements
    scenario.metadata.solver_override = Solver.AMGCL_amg_rigid_body # iterative solver
    

  • Similarly, Modal scenario is created using ModalScenarioDescriptor which takes boundary_conditions, a list of boundary conditions, as input.

    • The metadata is the same as for the linear elastic scenario except there is necessary input modal_metadata.desired_eigenvalues that takes in an integer for the number of eigenvalues.

    # Setup the modal simulation scenario
    modal_scenario = ModalScenarioDescriptor()
    modal_scenario.materials = {"Aluminum 6061-T6": material}
    modal_scenario.boundary_conditions = [fixed_boundary]
    
    modal_scenario.metadata.resolution = 10000
    modal_scenario.modal_metadata.desired_eigenvalues = 10
    
    # Optional settings
    modal_scenario.metadata.basis_order = 2 # 2 = quadratic elements
    modal_scenario.metadata.solver_override = Solver.AMGCL_amg_rigid_body  #iterative solver
    

  • Similarly, Linear Buckling scenario is created using LinearBucklingScenarioDescriptor which takes boundary_conditions, a list of boundary conditions, as input. A linear buckling scenario first performs a linear elasticity simulation, and then solves a generalized eigenvalue problem to determine the critical load factors.

    • The metadata is the same as for the linear elastic scenario except there is necessary input buckling_metadata.desired_eigenvalues that takes in an integer for the number of eigenvalues.

    # Setup the linear buckling simulation scenario
    buckling_scenario = LinearBucklingScenarioDescriptor()
    buckling_scenario.materials = {"Aluminum 6061-T6": material}
    buckling_scenario.boundary_conditions = [fixed_boundary]
    
    buckling_scenario.metadata.resolution = 10000
    # Generally, the first critical load factor from a buckling analysis is most useful,
    # a small number of additional buckling modes may be useful to examine.
    buckling_scenario.buckling_metadata.desired_eigenvalues = 3
    
    # Optional settings
    buckling_scenario.metadata.basis_order = 2 # 2 = quadratic elements
    # Use iterative solver for the linear elasticity simulation
    buckling_scenario.metadata.solver_override = Solver.AMGCL_amg_rigid_body
    

  • Similarly, a Thermal scenario is created using StaticThermalScenarioDescriptor which takes a list of boundary_conditions and a list of internal_conditions as input. Note that the material must be a thermal material for this scenario type.

    • The metadata is the same as for the linear elastic scenario with additional input thermal_metadata.environment_temperature which specifies the temperature of the environment.

    # Setup the thermal simulation scenario
    thermal_scenario = StaticThermalScenarioDescriptor()
    thermal_scenario.materials = {"Aluminum 6061-T6": material}
    thermal_scenario.boundary_conditions = [fixed_temp, convection]
    thermal_scenario.internal_conditions = [constant_heat]
    thermal_scenario.thermal_metadata.environment_temperature = 0.0
    
    thermal_scenario.metadata.resolution = 10000
    
    # Optional settings
    thermal_scenario.metadata.basis_order = 2 # 2 = quadratic elements
    thermal_scenario.metadata.solver_override = Solver.AMGCL_amg # iterative solver
    

• Finalize Simulation Setup and Execute

  • Stress Simulation is created using StressSimulator. It takes the following inputs:

    • An assembly or the list of MaterialDomains

    • LinearElasticScenarioDescriptor that describes the simulation scenario

    # Initialize and run the linear elastic scenario
    simulator = StressSimulator(assembly, scenario)
    simulator.solve()
    

  • Modal Simulation is created using ModalSimulator. It takes the following inputs:

    • An assembly or the list of MaterialDomains

    • ModalScenarioDescriptor that describes the simulation scenario

    # Initialize and run the modal scenario
    modal_simulator = ModalSimulator(assembly, modal_scenario)
    modal_simulator.solve()
    

  • Linear Buckling Simulation is created using LinearBucklingSimulator. It takes the following inputs:

    • An assembly or the list of MaterialDomains

    • LinearBucklingScenarioDescriptor that describes the simulation scenario

    # Initialize and run the linear buckling scenario
    buckling_simulator = LinearBucklingSimulator(assembly, buckling_scenario)
    buckling_simulator.solve()
    

  • Static Thermal Simulation is created using StaticThermalSimulator. It takes the following inputs:

    • An assembly or the list of MaterialDomains

    • StaticThermalScenarioDescriptor that describes the simulation scenario

    # Initialize and run the static thermal scenario
    thermal_simulator = StaticThermalSimulator(assembly, thermal_scenario)
    thermal_simulator.solve()
    

Query and Result Output

• Create a QueryResult object that specifies the domain to be sampled. It takes a single component or an assembly as input

  assembly = [component_1, component_2]
  
  results = QueryResult(component_1) # sample only the defined component
  results = QueryResult(assembly) # sample on the full assembly
  

• Specify a query class & type that defines what quantity you are querying. The query classes and types available are as follows :

  • Global Query is for querying quantities defined for the entire domain or structure. The available GlobalQueryType``s are ``Compliance for linear elasticity simulations, Frequency for modal simulations, CriticalLoadFactor for linear buckling simulations, and ThermalCompliance for static thermal simulations.

  # GlobalQuery example
  # Create a GlobalQuery to query the first mode
  results = QueryResult(assembly)
  frequency_query = GlobalQuery(GlobalQueryType.Frequency, DiscreteIndex(0))
  
  # Sample the simulation for the specified query
  freq_1 = simulator.sample(mode1_query, results)
  print(freq_1.get(0,0)) # Print out the first mode in Hz
  
  # Example use to get the first 10 modes
  # Note, the simulator would need this:
  # modal_scenario.modal_metadata.desired_eigenvalues = 10
  
  for i in range(10):
      # Create the query for the current mode
      query = GlobalQuery(GlobalQueryType.Frequency, DiscreteIndex(i))
  
      # Sample the query using the simulator
      r = simulator.sample(query, results)
  
      # Append the value obtained from result.get(0, 0) to a list
      values.append(r.get(0, 0))
      print(values) # list of first 10 modes in Hz
  

  • Field Query is for querying field quantities defined at any point in the domain. This depends on the physics type.

    • For Stress Simulation, the following FieldType inputs are available

      • Displacement tuple of dimension 3 for each point

        [x-displacement, y-displacement, z-displacement]
        

      • Strain and Stress are each a tuple of dimension 6 for each point

        [stress_xx, stress_yy, stress_zz, stress_yz, stress_xz, stress_xy]
        [strain_xx, strain_yy, strain_zz, strain_yz, strain_xz, strain_xy]
        

      • TopologicalSensitivity, StrainEnergyDensity, and VonMisesStress are each a tuple of dimension 1 for each point

    • For Modal Simulation and Linear Buckling Simulation, the only FieldType available is Displacement, which corresponds to the mode shape information and requires an index for the specified mode.

      # Create a field query for the displacement/mode shape of mode 1
      mode_num = 0 # first mode
      modal_field_query = FieldQuery(f=Field.Displacement, scheme=DiscreteIndex(mode_num))
      

    • For Thermal Simulation, the available FieldType are * Temperature which is a tuple of dimension 1

      • HeatFlux which is a tuple of dimension 3

    • Field Query Interface usage

      # FieldQuery example
      # Create a FieldQuery to query displacement
      results = QueryResult(assembly)
      displacement_query = FieldQuery(f=Field.Displacement)
      
      # Create a FieldQuery with optional input to query the y-component
      y_displacement_query = FieldQuery(f=Field.Displacement, component=1)
      
      # FieldQuery with optional input to query the norm of displacement
      displacement_magnitude_query = FieldQuery(f=Field.Displacement, norm=True)
      
      # Sample the simulation with the given query to create a VectorArray
      # * more info on VectorArrays are provided in the subsequent section *
      displacement_VectorArray1 = simulator.sample(displacement_query, results)
      displacement_VectorArray2 = simulator.sample(y_displacement_query, results)
      displacement_VectorArray3 = simulator.sample(displacement_magnitude_query, results)
      
      displacement1 = displacement_VectorArray1.get(i, j); # jth component of the displacement (0 = X, 1 = Y, or 2 = Z) at the ith sampling point
      displacement2 = displacement_VectorArray2.get(1, 0); # y-displacement value at the second sampling point
      displacement3 = displacement_VectorArray3.get(1, 0); # displacement magnitude at the second sampling point
      

Sample and Export Results

• Sampling results in-memory

  • Methods to sample results of the simulation and store/print VectorArray results from the queries created above.

    # Query results
    results = QueryResult(assembly)
    stress_query = FieldQuery(f=Field.Stress) # a field is needed to create a field query
    stress_VectorArray = simulator.sample(stress_query, results)
    

  • This sampling is stored in a VectorArray of size n_tuple by dimension, n_tuple is the number of points sampled in the component/assembly and dimension depends on the query type. For example, FieldQuery(Field.Displacement) has a dimension of 3 for each displacement component, thus VectorArray.get(i, 2) would get the z-displacement of point i.

    # Get the number of tuples (vectors) in the VectorArray (one per point results are sampled at)
    n_tuples = stress_VectorArray.n_tuples()
    
    # Get the dimension of each vector - in this case 6, one value for each stress component
    dimension = stress_VectorArray.dimension()
    
    stress_xx = stress_VectorArray.get(i, 0); # stress_xx at i-th sample point
    stress_yy = stress_VectorArray.get(i, 1); # stress_yy at i-th sample point
    stress_ij = stress_VectorArray.get(i, j); # jth stress component (0, 1, 2, 3, 4, or 5) of the ith tuple
    

• Export results to a file: Once the query is created, results can be used to write a .vtu file which contains all solution fields via results.writeVTK which has two inputs:

  • file name "*.vtu"

  • UnitSystem is an attribute in the *.vtu ( Note that this argument is just metadata in the *.vtu file. The result magnitudes are in the unit of the scenario regardless of the unit system specified. )

  # Write results to a file (unit system is saved as a metadata in the vtu file)
  results = QueryResult(assembly)
  stress_query = FieldQuery(Field.VonMisesStress)
  simulator.sample(stress_query, results)
  
  # Write results to a file (unit system is saved as a metadata in the vtu file)
  results.writeVTK("cantilever_beam_results.vtu", UnitSystem.MeterKilogramSecond)
  

Sampling on Custom Geometries

• Another important functionality is the ability to specify a point set to sample on. This can be done by creating a new sampling_component consisting of vertices at the desired sampling locations. Note that the sampling_component is not to be included in the assembly that is being simulated.

  # Create an empty MeshModel to store vertices for desired sampling locations
  sampling_geometry = MeshModel()
  
  # Add the desired sampling point locations and store their ID for later use
  test_index1 = sampling_geometry.addVertex([6.0, 0.0, 0.0])
  test_index2 = sampling_geometry.addVertex([3.0, 0.0, 0.0])
  
  sampling_component = MaterialDomain(sample_geometry, "material_name", scenario)
  
  # Define the assembly to be simulated (note that the sampling_component is not included)
  component = MaterialDomain(component_geometry, "material_name", scenario)
  assembly = [component]
  simulator = StressSimulator(assembly, scenario)
  simulator.solve()
  
  displacement_query = FieldQuery(f=Field.Displacement)
  sampled_results = QueryResult(sampling_component) # sample only the custom point set
  sampled_VectorArray = simulator.sample(displacement_query, sampled_results)
  
  displacement_index1 = sampled_VectorArray.get(test_index1, 0) # x-displacement at above specified index1
  displacement_index2 = sampled_VectorArray.get(test_index2, 1) # y-displacement at above specified index2
  

LevelOpt

• LevelOpt is another ScenarioDescriptor (LevelOptScenarioDescriptor) for level set structural topology optimization. It supports all of the standard linear elastic boundary conditions and a variety of metadata optimization parameters.

  • LevelOpt scenario is created using LevelOptScenarioDescriptor. It takes the following shared inputs from previous ScenarioDescriptor classes:

    • boundary_conditions note, with LevelOpt, multiple load cases can be input to a single optimization scenario using unique load_case_id. This allows LevelOpt to consider each load individually rather than the net load.

      # multiple load case optimization needs unique load_case_ids
      fixed_boundary1 = FixedBoundaryDescriptor()
      fixed_boundary1.boundary = MeshModel("fixed.ply")
      fixed_boundary1.load_case_id = 0
      
      fixed_boundary2 = FixedBoundaryDescriptor()
      fixed_boundary2.boundary = MeshModel("fixed.ply")
      fixed_boundary2.load_case_id = 1
      
      fixed_boundaries = [fixed_boundary1, fixed_boundary2]
      
      load1 = VectorForceDescriptor()
      load1.boundary = MeshModel("load1.ply")
      load1.direction = [1, 0, 0]
      load1.magnitude = 1000
      load1.load_case_id = 0
      
      load2 = VectorForceDescriptor()
      load2.boundary = MeshModel("load1.ply")
      load2.direction = [0, 1, 0]
      load2.magnitude = 1000
      load2.load_case_id = 1
      
      load_cases = [load1, load2]
      

    • internal_conditions

    • resolution or cell_size

    • units

    • solver_type used in the initial linear elastic simulation with the following options:

      • MKL_PardisoLDLT (default) direct solver

      • AMGCL_amg_rigid_body iterative solver

    • basis_order used in the initial linear elastic simulation

LevelOpt Metadata Settings

• The set of metadata parameters specific to the LevelOptScenarioDescriptor include:

  • vol_frac_cons (volume fraction constraint) sets the target volume of the final design as a fraction of the initial volume.

  • voxelSize, or level set cell size, determines the size of the level set grid cells as a fraction of the FEA grid cell size.

  • move_limit controls the extent of changes per optimization step, as a factor of the voxelSize.

  • opt_max_iter (optimization max iterations) sets the maximum number of iterations for the optimization process. Each iteration refines the design by updating the topology based on the objective function and constraints.

  • fix_thickness specifies the region around boundary conditions that remains unchanged as a factor of level set grid cell size.

  • smooth_iter defines the frequency the geometry is smoothed during the optimization process as a number of iterations.

  • enable_fixed_interfaces allows specifying if the interface between the design domain (optimized) and non-design domain should be fully preserved. True preserves the interface and False allows material to be removed at the interface.

  • num_load_cases is a input required for an optimization scenario with multiple load_case_id. This enables individual load cases to be considered separately during optimization instead of a net load.

  # Setup the LevelOpt optimization scenario descriptor
  leveloptscenario = LevelOptScenarioDescriptor()
  
  # Standard scenario parameters
  leveloptscenario.materials = {"Aluminum": material}
  leveloptscenario.metadata.cell_size = 2.5
  leveloptscenario.metadata.units = UnitSystem.MeterKilogramSecond
  leveloptscenario.metadata.solver_override = "MKL_PardisoLDLT"   # direct LinearElastic solver
  leveloptscenario.metadata.basis_order = 1                   # linear order for LinearElastic solver
  leveloptscenario.boundary_conditions = fixed_boundaries + load_cases    # consist of two sets of load cases
  
  # LevelOpt specific metadata
  leveloptscenario.optimization_metadata.vol_frac_cons = 0.2  # 20% target volume
  leveloptscenario.optimization_metadata.voxelSize = 0.5      # 0.5 (default) level set cell size factor
  leveloptscenario.optimization_metadata.move_limit = 1.0     # 1.0 (default) move limit factor
  leveloptscenario.optimization_metadata.opt_max_iter = 10    # 10 design iterations
  leveloptscenario.optimization_metadata.fix_thickness = 4  # 4 (default) level set cells which are unchanged near BCs
  leveloptscenario.optimization_metadata.smooth_iter = 1      # smoothing performed at every 1 iteration
  leveloptscenario.optimization_metadata.enable_fixed_interfaces = True   # retain full contact between assembly components
  leveloptscenario.optimization_metadata.num_load_cases = 2   # 2 load cases created
  

LevelOpt Optimization and Supported Queries

• To run an optimization scenario a design domain Material Domain must first be defined. For assemblies, additional non-design Material Domain(s) need to be defined. Note, the design domain MeshModel requires a unique instance_id.

  # make sure to include instance_id for design_geometry/design_domain
  design_geometry = MeshModel("design_domain.stl")
  design_geometry.instance_id = "ex_design_domain"   # can be any unique id/name
  design_geometry.refine()
  
  # Create material domains and assembly
  design_domain = MaterialDomain(design_geometry, "Aluminum", leveloptscenario)
  component1 = MaterialDomain(body1, "Aluminum", leveloptscenario)
  component2 = MaterialDomain(body2, "Aluminum", leveloptscenario)
  
  assembly = [design_domain, component1, component2]
  

• The LevelOpt optimization solver then takes the following arguments:

  • design Material Domain

  • assembly or list of all design and non-design Material Domain(s)

  • LevelOptScenarioDescriptor which describes the optimization scenario parameters and inputs.

  • starting design MeshModel, or if no starting design is used, a None argument.

  # Create and run the simulation
  optimizer = LevelOpt(design_domain, assembly, leveloptscenario, None)
  optimizer.optimize()
  
  # Create and run the simulation with optional starting design
  starting_design = MeshModel("starting_design.ply")
  optimizer = LevelOpt(design_domain, assembly, leveloptscenario, starting_design)
  optimizer.optimize()
  

• LevelOpt Queries include support for two GlobalQueryType classes, Compliance and VolumeFraction and two FieldQuery classes, BoundarySensitivity and BoundaryVelocity. A discrete index argument corresponding to the optimization iteration is required to query for any of the previously described quantities.

  # LevelOpt global queries
  compliance_query = GlobalQuery(GlobalQueryType.Compliance, DiscreteIndex(iteration))
  volume_fraction_query = GlobalQuery(GlobalQueryType.VolumeFraction, DiscreteIndex(iteration))
  
  compliance = optimizer.sample(compliance_query)
  volume_fraction = optimizer.sample(volume_fraction_query)
  
  compliance_value = compliance.get(0,0)
  volume_fraction_value = volume_fraction.get(0,0)
  
  # LevelOpt field queries
  boundary_sensitivity_query = FieldQuery(f=Field.BoundarySensitivity, scheme=DiscreteIndex(0))
  boundary_velocity_query = FieldQuery(f=Field.BoundaryVelocity, scheme=DiscreteIndex(0))
  
  r_s = optimizer.sample(boundary_sensitivity_query) # no QueryResult() object, only defined on design boundary
  r_v = optimizer.sample(boundary_velocity_query) # no QueryResult() object, only defined on design boundary
  
  sensitivity_value = r_s.get(0,0) # tuple of dimension 1 for each point
  velocity_value = r_v.get(0,0) # tuple of dimension 1 for each point
  

• Output mesh designs can be stored using .getDesigns() from the LevelOpt scenario optimization and saved optionally as .ply files in a desired directory with .writePLY()

# run LevelOpt Scenario
optimizer = LevelOpt(design_domain, assembly, leveloptscenario, None)
optimizer.optimize()

# get the design iterations from the optimization scenario
designs = optimizer.getDesigns()

# write the design output mesh .ply files to desired directory
for i, design in enumerate(designs):
   ply_filename = os.path.join(ply_dir, f"optimized_design_{i}.ply")
   design.writePLY(ply_filename)

Examples

Complete examples of performing a structural simulation, a modal simulation, or a static thermal simulation are available, as well as the required geometry files: beam geometry, restraint geometry, load geometry.

A complete example of performing a LevelOpt bracket optimization is available along with all geometry files: LevelOpt bracket geometry, restraint, load1, load2, and a reference starting design.