Getting Started with PyDAGMC

Getting Started with PyDAGMC#

Welcome to the pydagmc tutorial! PyDAGMC is a Python module designed to simplify the management and manipulation of hierarchical geometric data structures used in computational modeling and simulation.

This guide will walk you through the essentials:

  1. Installation & Setup: Get up and running in minutes.

  2. Loading a Model: Import your DAGMC geometry files.

  3. Exploring the Model: Navigate volumes, surfaces, and groups with ease.

  4. Accessing Entity Properties: Retrieve IDs, names, materials, and geometric data.

  5. Modifying Geometry: Alter entity properties and relationships on the fly.

  6. Creating New Geometry: Programmatically build models from scratch.

  7. Working with Triangles: Access the underlying mesh connectivity and coordinates.

  8. Exporting Data: Save your work for visualization or further analysis.

  9. Deleting Entities: Removing parts of your model.

Ready to dive in? Let’s get started!

1. Installation & Setup#

PyDAGMC depends on MOAB. Please ensure MOAB is installed in your Python environment. You can find MOAB installation instructions here.

This tutorial also uses numpy for some data handling.

[1]:

import pydagmc
import numpy as np
from pathlib import Path
import urllib.request
import os

# Helper function to download files if they don't exist
def download_file(url, filename_str):
    filename = Path(filename_str)
    if not filename.exists():
        print(f"Downloading {filename.name} from {url}...")
        try:
            urllib.request.urlretrieve(url, filename)
            print(f"{filename.name} downloaded successfully.")
        except Exception as e:
            print(f"Error downloading {filename.name}: {e}")
            raise
    else:
        print(f"{filename.name} already exists.")
    return filename

# Example model files
FUEL_PIN_URL = 'https://tinyurl.com/y3ugwz6w' # fuel_pin.h5m
FUEL_PIN_FILE = download_file(FUEL_PIN_URL, "fuel_pin.h5m")

CUBE_STL_URL = 'https://raw.githubusercontent.com/svalinn/pydagmc/master/test/cube.stl' # A simple cube STL
CUBE_STL_FILE = download_file(CUBE_STL_URL, "cube.stl")

Downloading fuel_pin.h5m from https://tinyurl.com/y3ugwz6w...
fuel_pin.h5m downloaded successfully.
Downloading cube.stl from https://raw.githubusercontent.com/svalinn/pydagmc/master/test/cube.stl...
cube.stl downloaded successfully.

2. Loading a Model#

The primary interface to your geometry is the pydagmc.Model class. You can create a Model instance by providing the path to a DAGMC file (e.g., .h5m) or by passing an existing pymoab.core.Core object.

[2]:

# Load the example fuel pin model
try:
    model = pydagmc.Model(str(FUEL_PIN_FILE))
    print("Model loaded successfully.")
except ModuleNotFoundError as e:
    print(f"Failed to load model: MOAB might not be installed or found. {e}")
except RuntimeError as e:
    print(f"Runtime error loading model (file might be corrupted or MOAB issue): {e}")
except Exception as e:
    print(f"An unexpected error occurred during model loading: {e}")

# Get a summary of the model's contents
print(repr(model))

Model loaded successfully.
Model: 4 Volumes, 21 Surfaces, 5 Groups

The repr(model) output (e.g., Model: 4 Volumes, 21 Surfaces, 5 Groups) gives a quick overview of the main geometric entities in the loaded file.

3. Exploring the Model#

Once a model is loaded, you can access its constituent Volumes, Surfaces, and Groups through various properties of the Model object.

3.1 Volumes#

Volumes represent distinct 3D regions in your geometry, often associated with materials.

[3]:

# List of all Volume objects
all_volumes = model.volumes
print(f"Number of volumes: {len(all_volumes)}")

# Dictionary mapping volume ID (int) to Volume object
volumes_by_id = model.volumes_by_id
if 1 in volumes_by_id:
    vol1 = volumes_by_id[1]
    print(f"\nDetails for Volume ID {vol1.id}:")
    print(f"  Object: {vol1}")
    print(f"  Material: {vol1.material}")
    print(f"  Number of surfaces: {len(vol1.surfaces)}")

# Dictionary mapping material name (str) to a list of Volume objects
volumes_by_material = model.volumes_by_material
print("\nVolumes grouped by material:")
for material_name, vols_list in volumes_by_material.items():
    print(f"  Material '{material_name}': IDs {[v.id for v in vols_list]}")

# Find volumes by a specific material name (case-sensitive)
fuel_volumes = model.find_volumes_by_material('fuel')
print(f"\nVolumes with material 'fuel': IDs {[v.id for v in fuel_volumes]}")

# List of Volume objects that do not have an assigned material group
volumes_without_material = model.volumes_without_material
print(f"\nNumber of volumes without an assigned material: {len(volumes_without_material)}")

Number of volumes: 4

Details for Volume ID 1:
  Object: Volume 1, 4092 triangles
  Material: fuel
  Number of surfaces: 3

Volumes grouped by material:
  Material 'fuel': IDs [np.int32(1), np.int32(2)]
  Material '41': IDs [np.int32(3)]
  Material 'Graveyard': IDs [np.int32(6)]

Volumes with material 'fuel': IDs [np.int32(1), np.int32(2)]

Number of volumes without an assigned material: 0

3.2 Surfaces#

Surfaces are 2D entities that typically form the boundaries of volumes. Each surface has a “sense” indicating which volumes lie on its forward and reverse sides (relative to the surface normal).

[4]:

# List of all Surface objects
all_surfaces = model.surfaces
print(f"Number of surfaces: {len(all_surfaces)}")

# Dictionary mapping surface ID (int) to Surface object
surfaces_by_id = model.surfaces_by_id
if 1 in surfaces_by_id:
    surf1 = surfaces_by_id[1]
    print(f"\nDetails for Surface ID {surf1.id}:")
    print(f"  Object: {surf1}")
    print(f"  Number of triangles: {surf1.num_triangles}")

    fwd_vol = surf1.forward_volume
    rev_vol = surf1.reverse_volume
    print(f"  Forward volume ID: {fwd_vol.id if fwd_vol else 'None'}")
    print(f"  Reverse volume ID: {rev_vol.id if rev_vol else 'None'}")

    # List of Volume objects that this surface bounds
    parent_volumes = surf1.volumes
    print(f"  Parent volume IDs: {[v.id for v in parent_volumes]}")

Number of surfaces: 21

Details for Surface ID 1:
  Object: Surface 1, 2048 triangles
  Number of triangles: 2048
  Forward volume ID: 1
  Reverse volume ID: 2
  Parent volume IDs: [np.int32(1), np.int32(2)]

3.3 Groups#

Groups are named collections of volumes and/or surfaces. They are commonly used to assign properties like materials (e.g., a group named mat:fuel containing all fuel volumes) or boundary conditions.

[5]:

# List of all Group objects
all_groups = model.groups
print(f"Number of groups: {len(all_groups)}")

# Dictionary mapping group name (str) to Group object
groups_by_name = model.groups_by_name
if 'mat:fuel' in groups_by_name:
    fuel_group = groups_by_name['mat:fuel']
    print(f"\nDetails for Group '{fuel_group.name}':")
    # print(repr(fuel_group)) # For a more verbose output
    print(f"  ID: {fuel_group.id}")
    print(f"  Contained volume IDs: {fuel_group.volume_ids}")
    print(f"  Contained surface IDs: {fuel_group.surface_ids}")

    # Check if a specific volume is part of this group
    if 1 in model.volumes_by_id:
        vol1_check = model.volumes_by_id[1]
        print(f"  Is Volume {vol1_check.id} in group '{fuel_group.name}'? {vol1_check in fuel_group}")

# List of all unique group names
print(f"\nAll group names in model: {list(model.group_names)}")

Number of groups: 5

Details for Group 'mat:fuel':
  ID: 2
  Contained volume IDs: [1 2]
  Contained surface IDs: []
  Is Volume 1 in group 'mat:fuel'? True

All group names in model: ['picked', 'mat:fuel', 'mat:Graveyard', 'mat:41', 'temp:300']

4. Accessing Entity Properties#

Volume, Surface, and Group objects are all derived from a common base class GeometrySet. They share fundamental properties like ID, category, and geometric dimension, alongside their specific attributes.

[6]:

if 1 in model.volumes_by_id:
    vol1 = model.volumes_by_id[1]
    print(f"Common Properties for Volume ID {vol1.id}:")
    print(f"  Category Tag: '{vol1.category}'")
    print(f"  Geometric Dimension Tag: {vol1.geom_dimension}")
    print(f"Specific Properties for Volume ID {vol1.id}:")
    print(f"  Calculated geometric volume: {vol1.volume:.2f}")
    print(f"  Total triangles in bounding surfaces: {vol1.num_triangles}")

if 1 in model.surfaces_by_id:
    surf1 = model.surfaces_by_id[1]
    print(f"\nCommon Properties for Surface ID {surf1.id}:")
    print(f"  Category Tag: '{surf1.category}'")
    print(f"  Geometric Dimension Tag: {surf1.geom_dimension}")
    print(f"Specific Properties for Surface ID {surf1.id}:")
    print(f"  Calculated geometric area: {surf1.area:.2f}")
    print(f"  Number of triangles: {surf1.num_triangles}")

if 'mat:fuel' in model.groups_by_name:
    fuel_group = model.groups_by_name['mat:fuel']
    print(f"\nCommon Properties for Group '{fuel_group.name}':")
    print(f"  Name: {fuel_group.name}") # Specific to Group
    print(f"  ID: {fuel_group.id}")
    print(f"  Category Tag: '{fuel_group.category}'")
    print(f"  Geometric Dimension Tag: {fuel_group.geom_dimension}")

Common Properties for Volume ID 1:
  Category Tag: 'Volume'
  Geometric Dimension Tag: 3
Specific Properties for Volume ID 1:
  Calculated geometric volume: 6157.48
  Total triangles in bounding surfaces: 4092

Common Properties for Surface ID 1:
  Category Tag: 'Surface'
  Geometric Dimension Tag: 2
Specific Properties for Surface ID 1:
  Calculated geometric area: 1759.29
  Number of triangles: 2048

Common Properties for Group 'mat:fuel':
  Name: mat:fuel
  ID: 2
  Category Tag: 'Group'
  Geometric Dimension Tag: 4

Note on ``category`` and ``geom_dimension`` tags: pydagmc strives for consistency. If it loads an entity that, for instance, has a geometric dimension appropriate for a surface but lacks the ‘Surface’ category tag, pydagmc will automatically assign the correct category tag when you access it as a pydagmc.Surface object (often with a warning). This helps ensure the model data is well-formed.

5. Modifying Geometry#

pydagmc enables you to alter various aspects of your geometric model programmatically.

5.1 Changing Entity IDs#

You can change the ID of a Volume, Surface, or Group. pydagmc ensures that new IDs are unique within their respective entity types to prevent conflicts. Assigning None to an ID will cause pydagmc to automatically assign the next available highest ID.

[7]:

if 2 in model.volumes_by_id:
    vol_to_modify = model.volumes_by_id[2]
    original_id = vol_to_modify.id
    new_id_candidate = 999
    print(f"Volume ID {original_id} initial ID: {vol_to_modify.id}")

    # Change to a new, unused ID
    if new_id_candidate not in model.volumes_by_id:
        vol_to_modify.id = new_id_candidate
        print(f"ID changed to: {vol_to_modify.id}")
        assert original_id not in model.volumes_by_id # Old ID is now free
        assert new_id_candidate in model.volumes_by_id # New ID is registered
    else:
        print(f"Could not change to {new_id_candidate} as it's already in use.")

    # Attempting to assign an already used ID raises a ValueError
    if 1 in model.volumes_by_id: # Assuming Volume 1 exists
        try:
            print(f"\nAttempting to set ID to 1 (which is used by another volume)...")
            vol_to_modify.id = 1
        except ValueError as e:
            print(f"  Error: {e}")

    # Assign None to auto-assign the next available ID
    current_max_id = max(model.used_ids[pydagmc.Volume], default=0)
    vol_to_modify.id = None
    print(f"\nID after assigning None: {vol_to_modify.id} (expected {current_max_id + 1})")
    assert vol_to_modify.id == current_max_id + 1

    # Revert ID for consistency in subsequent tutorial cells
    # This logic finds a safe ID to revert to, preferring original_id if free.
    if original_id not in model.volumes_by_id:
         vol_to_modify.id = original_id
    else:
        fallback_id = 1000
        while fallback_id in model.volumes_by_id: fallback_id +=1
        vol_to_modify.id = fallback_id
    print(f"ID reset for consistency: {vol_to_modify.id}")

Volume ID 2 initial ID: 2
ID changed to: 999

Attempting to set ID to 1 (which is used by another volume)...
  Error: Volume ID 1 is already in use in this model.

ID after assigning None: 1000 (expected 1000)
ID reset for consistency: 2

5.2 Changing a Volume’s Material#

Modifying a Volume’s material property will:

  1. Update the material name associated with the volume.

  2. Remove the volume from its previous material group (e.g., mat:old_material).

  3. If a group for the new material (e.g., mat:new_material) doesn’t exist, it will be created.

  4. Add the volume to the new material group.

[8]:

if 1 in model.volumes_by_id:
    vol1 = model.volumes_by_id[1]
    original_material = vol1.material
    print(f"Volume {vol1.id} initial material: '{original_material}'")

    new_material_name = "graphite"
    vol1.material = new_material_name
    print(f"Volume {vol1.id} new material: '{vol1.material}'")

    # Verify changes
    assert vol1.material == new_material_name
    new_group_name = f"mat:{new_material_name}"
    assert new_group_name in model.groups_by_name
    assert vol1 in model.groups_by_name[new_group_name].volumes
    print(f"Volume {vol1.id} is now in group '{new_group_name}'.")

    old_group_name = f"mat:{original_material}"
    if old_group_name in model.groups_by_name:
        assert vol1 not in model.groups_by_name[old_group_name].volumes
        print(f"Volume {vol1.id} removed from group '{old_group_name}'.")

    # Revert material for tutorial consistency
    vol1.material = original_material
    print(f"\nVolume {vol1.id} material reverted to '{vol1.material}'.")

Volume 1 initial material: 'fuel'
Volume 1 new material: 'graphite'
Volume 1 is now in group 'mat:graphite'.
Volume 1 removed from group 'mat:fuel'.

Volume 1 material reverted to 'fuel'.

5.3 Changing Surface Senses#

You can change which volumes a Surface separates by setting its forward_volume, reverse_volume, or both via the senses property. This also updates the underlying parent-child relationships in MOAB.

[9]:

if 1 in model.surfaces_by_id and 1 in model.volumes_by_id and 3 in model.volumes_by_id:
    surf_to_modify = model.surfaces_by_id[1]
    vol_sense_a = model.volumes_by_id[1]
    vol_sense_b = model.volumes_by_id[3] # A different volume

    original_fwd = surf_to_modify.forward_volume
    original_rev = surf_to_modify.reverse_volume
    print(f"Surface {surf_to_modify.id} original senses: Fwd={original_fwd.id if original_fwd else 'N'}, Rev={original_rev.id if original_rev else 'N'}")

    # Change forward volume
    surf_to_modify.forward_volume = vol_sense_b
    print(f"After changing fwd_vol: Fwd={surf_to_modify.forward_volume.id}, Rev={surf_to_modify.reverse_volume.id if surf_to_modify.reverse_volume else 'N'}")

    # Change reverse volume
    surf_to_modify.reverse_volume = vol_sense_a
    print(f"After changing rev_vol: Fwd={surf_to_modify.forward_volume.id}, Rev={surf_to_modify.reverse_volume.id}")

    # Revert to original senses for tutorial consistency
    surf_to_modify.senses = [original_fwd, original_rev]
    print(f"Senses reverted: Fwd={surf_to_modify.forward_volume.id if surf_to_modify.forward_volume else 'N'}, Rev={surf_to_modify.reverse_volume.id if surf_to_modify.reverse_volume else 'N'}")

Surface 1 original senses: Fwd=1, Rev=2
After changing fwd_vol: Fwd=3, Rev=2
After changing rev_vol: Fwd=3, Rev=1
Senses reverted: Fwd=1, Rev=2

5.4 Modifying Groups#

You can rename Group objects and manage the entities (volumes or surfaces) they contain.

Group Merging: If two distinct group entities in the underlying MOAB data happen to have the same name, accessing the model.groups_by_name property will cause pydagmc to automatically merge them. The entities from one group will be moved into the other, and the redundant group entity will be deleted. This ensures that model.groups_by_name always provides a unique Group object for each name.

[10]:

if 'mat:fuel' in model.groups_by_name and 3 in model.volumes_by_id:
    target_group = model.groups_by_name['mat:fuel']
    vol_for_group_ops = model.volumes_by_id[3] # Typically in 'mat:41'

    original_group_name = target_group.name
    print(f"Working with group: '{original_group_name}', ID: {target_group.id}")
    print(f"Initial volumes: {target_group.volume_ids}")

    # Change group name
    candidate_new_name = "core_ fissile_material"
    if candidate_new_name not in model.group_names:
        target_group.name = candidate_new_name
        print(f"Group name changed to: '{target_group.name}'")
        assert candidate_new_name in model.groups_by_name
        assert original_group_name not in model.groups_by_name
    else:
        print(f"Skipping name change to '{candidate_new_name}' as it's already in use.")

    # Add a volume to the group (use current name via target_group.name)
    active_target_group = model.groups_by_name[target_group.name] # Re-fetch by current name
    if vol_for_group_ops not in active_target_group:
        active_target_group.add_set(vol_for_group_ops)
        print(f"Added Volume {vol_for_group_ops.id} to group '{active_target_group.name}'. Volumes: {active_target_group.volume_ids}")

    # Remove the volume from the group
    active_target_group.remove_set(vol_for_group_ops)
    print(f"Removed Volume {vol_for_group_ops.id} from group '{active_target_group.name}'. Volumes: {active_target_group.volume_ids}")

    # Revert name for tutorial consistency
    if active_target_group.name != original_group_name:
         active_target_group.name = original_group_name
         print(f"Group name reverted to: '{active_target_group.name}'")
    final_target_group = model.groups_by_name[original_group_name] # ensure we have the correct object reference

    original_fuel_vols = set(final_target_group.volume_ids)
    print(f"Group '{final_target_group.name}' (handle {final_target_group.handle}) initially has volumes: {original_fuel_vols}")

    # Create a temporary group
    temp_group = model.create_group(name="temp_for_merge", group_id=888)
    vol_to_add_to_temp = None
    if 4 in model.volumes_by_id: # Assuming vol ID 4 exists and is not in 'mat:fuel'
        vol_to_add_to_temp = model.volumes_by_id[4]
        if vol_to_add_to_temp.id not in original_fuel_vols:
            temp_group.add_set(vol_to_add_to_temp)
            print(f"Created group '{temp_group.name}' (ID {temp_group.id}) with volume {vol_to_add_to_temp.id}.")
        else:
            vol_to_add_to_temp = None # Avoid confusion if already present
            print(f"Volume {model.volumes_by_id[4].id} already in {final_target_group.name}, merge won't show addition.")

    print("Accessing `model.groups_by_name` to trigger merge...")
    final_groups_map = model.groups_by_name
    merged_group = final_groups_map[final_target_group.name]

    print(f"After merge, group '{merged_group.name}' (handle {merged_group.handle}) exists.")
    assert merged_group.handle == final_target_group.handle # Original group should be preserved

    final_merged_vols = set(merged_group.volume_ids)
    print(f"Its volumes are now: {final_merged_vols}")
    if vol_to_add_to_temp:
        assert vol_to_add_to_temp.id in final_merged_vols
        print(f"Volume {vol_to_add_to_temp.id} (from temp group) is now in the merged group.")
        # Clean up by removing the added volume
        merged_group.remove_set(vol_to_add_to_temp)
        print(f"Cleaned up: Removed {vol_to_add_to_temp.id} from '{merged_group.name}'.")

Working with group: 'mat:fuel', ID: 2
Initial volumes: [1 2]
Group name changed to: 'core_ fissile_material'
Added Volume 3 to group 'core_ fissile_material'. Volumes: [1 2 3]
Removed Volume 3 from group 'core_ fissile_material'. Volumes: [1 2]
Group name reverted to: 'mat:fuel'
Group 'mat:fuel' (handle 12682136550675316816) initially has volumes: {np.int32(1), np.int32(2)}
Accessing `model.groups_by_name` to trigger merge...
After merge, group 'mat:fuel' (handle 12682136550675316816) exists.
Its volumes are now: {np.int32(1), np.int32(2)}

6. Creating New Geometry#

pydagmc allows you to build new geometric models from scratch or add to existing ones.

[11]:

# For clarity, let's work with a new, empty model for creation tasks
creative_model = pydagmc.Model()
print(f"New empty model for creation: {creative_model}")

# Create a Volume
created_vol_id = 101
my_volume = creative_model.create_volume(global_id=created_vol_id)
print(f"Created Volume: ID {my_volume.id}, Material: {my_volume.material}")
assert my_volume in creative_model.volumes_without_material # New volumes have no material

# Create another Volume (ID will be auto-assigned)
another_my_volume = creative_model.create_volume()
print(f"Created another Volume: ID {another_my_volume.id}")
assert another_my_volume.id == created_vol_id + 1

# Create a Surface by loading geometry from an STL file
created_surf_id = 201
my_surface = creative_model.create_surface(global_id=created_surf_id, filename=CUBE_STL_FILE)
print(f"Created Surface from STL: ID {my_surface.id}, Triangles: {my_surface.num_triangles}")
assert my_surface.num_triangles == 12 # The example cube.stl has 12 triangles

# Create an empty Surface
empty_surface = creative_model.create_surface(global_id=202)
print(f"Created empty Surface: ID {empty_surface.id}, Triangles: {empty_surface.num_triangles}")

# Create a Group
my_group_name = "custom_assembly"
my_group_id = 301
my_group = creative_model.create_group(name=my_group_name, group_id=my_group_id)
print(f"Created Group: Name '{my_group.name}', ID {my_group.id}")

# Establish relationships
my_surface.senses = [my_volume, another_my_volume] # Surface separates the two volumes
print(f"Surface {my_surface.id} senses set: Fwd={my_surface.forward_volume.id}, Rev={my_surface.reverse_volume.id}")

my_group.add_set(my_volume)
my_group.add_set(my_surface)
print(f"Group '{my_group.name}' contents: Vols={my_group.volume_ids}, Surfs={my_group.surface_ids}")

# Assign material to one of the new volumes
my_volume.material = "custom_material"
print(f"Volume {my_volume.id} material set to '{my_volume.material}'")
assert "custom_material" in creative_model.volumes_by_material

# Batch-add groups using `model.add_groups()`
batch_group_data = {
    ("boundary_elements", 401): [my_surface], # Can use object or ID
    ("other_material", 402): [another_my_volume.id]
}
creative_model.add_groups(batch_group_data)
print("Batch-added groups:")
print(f"  Group 'boundary_elements' surfaces: {creative_model.groups_by_name['boundary_elements'].surface_ids}")
print(f"  Group 'other_material' volumes: {creative_model.groups_by_name['other_material'].volume_ids}")

New empty model for creation: Model: 0 Volumes, 0 Surfaces, 0 Groups
Created Volume: ID 101, Material: None
Created another Volume: ID 102
Created Surface from STL: ID 201, Triangles: 12
Created empty Surface: ID 202, Triangles: 0
Created Group: Name 'custom_assembly', ID 301
Surface 201 senses set: Fwd=101, Rev=102
Group 'custom_assembly' contents: Vols=[101], Surfs=[201]
Volume 101 material set to 'custom_material'
Batch-added groups:
  Group 'boundary_elements' surfaces: [201]
  Group 'other_material' volumes: [102]

7. Working with Triangles#

Geometric surfaces in DAGMC are represented by triangular facets. pydagmc provides access to the connectivity (which vertices form each triangle) and coordinates of these triangles. This data can be retrieved for individual surfaces, or aggregated for volumes or groups.

[12]:

# Using the 'creative_model' and the cube surface created earlier (ID was 201)
if created_surf_id in creative_model.surfaces_by_id:
    cube_surface_from_creative_model = creative_model.surfaces_by_id[created_surf_id]
    print(f"Triangle data for Surface {cube_surface_from_creative_model.id} (the cube):")

    # Raw MOAB entity handles for each triangle in the surface
    triangle_handles = cube_surface_from_creative_model.triangle_handles
    print(f"  Number of triangle handles: {len(triangle_handles)}")

    # Connectivity: NumPy array (N_triangles, 3) of MOAB vertex entity handles for each triangle
    vertex_handles_per_triangle = cube_surface_from_creative_model.triangle_conn
    print(f"  Shape of vertex handles per triangle: {vertex_handles_per_triangle.shape}")

    # Coordinates: NumPy array (N_triangles * 3, 3) of XYZ coordinates for all triangle vertices (expanded list)
    expanded_triangle_vertex_coords = cube_surface_from_creative_model.triangle_coords
    print(f"  Shape of expanded triangle vertex coordinates: {expanded_triangle_vertex_coords.shape}")

    # Get connectivity and coordinates optimized for space or direct indexing
    # `compress=True`: Returns unique vertex coordinates and connectivity as indices into that unique array.
    conn_indices_unique, unique_coords = cube_surface_from_creative_model.get_triangle_conn_and_coords(compress=True)
    print(f"\nCompressed: unique_coords shape {unique_coords.shape}, conn_indices (to unique_coords) shape {conn_indices_unique.shape}")
    if conn_indices_unique.size > 0:
        first_tri_coords_via_unique = unique_coords[conn_indices_unique[0]]
        print(f"Coordinates of 1st triangle (from unique):\n{first_tri_coords_via_unique}")

    # `compress=False`: Returns expanded coordinates (each triangle's vertices listed independently)
    # and connectivity as simple indices (0,1,2 for tri1; 3,4,5 for tri2 etc.) into that expanded array.
    conn_indices_expanded, expanded_coords_again = cube_surface_from_creative_model.get_triangle_conn_and_coords(compress=False)
    print(f"Expanded: expanded_coords shape {expanded_coords_again.shape}, conn_indices (to expanded_coords) shape {conn_indices_expanded.shape}")
    if conn_indices_expanded.size > 0:
        first_tri_coords_via_expanded = expanded_coords_again[conn_indices_expanded[0]]
        print(f"Coordinates of 1st triangle (from expanded):\n{first_tri_coords_via_expanded}")

    # Mapping: Get a dictionary from triangle EntityHandle to its vertex indices (into a coordinate array)
    triangle_to_coord_indices_map, map_coords_array = cube_surface_from_creative_model.get_triangle_coordinate_mapping(compress=True)
    if len(triangle_to_coord_indices_map) > 0:
        first_tri_handle_example = triangle_handles[0]
        indices_for_handle = triangle_to_coord_indices_map[first_tri_handle_example]
        coords_for_handle = map_coords_array[indices_for_handle]
        print(f"\nCoordinates of 1st triangle (handle {first_tri_handle_example}) via mapping:\n{coords_for_handle}")

Triangle data for Surface 201 (the cube):
  Number of triangle handles: 12
  Shape of vertex handles per triangle: (12, 3)
  Shape of expanded triangle vertex coordinates: (36, 3)

Compressed: unique_coords shape (8, 3), conn_indices (to unique_coords) shape (12, 3)
Coordinates of 1st triangle (from unique):
[[-35.  60.  20.]
 [-55.  60.  20.]
 [-35.  40.  20.]]
Expanded: expanded_coords shape (36, 3), conn_indices (to expanded_coords) shape (12, 3)
Coordinates of 1st triangle (from expanded):
[[-35.  60.  20.]
 [-55.  60.  20.]
 [-35.  40.  20.]]

Coordinates of 1st triangle (handle 2305843009213693953) via mapping:
[[-35.  60.  20.]
 [-55.  60.  20.]
 [-35.  40.  20.]]

8. Exporting Data#

You can save your entire pydagmc.Model to an H5M file or export specific entities (like a Surface, Volume, or Group) to a VTK file for visualization in tools like ParaView or VisIt.

[13]:

# Using the 'creative_model' from section 6
output_h5m_path = Path("pydagmc_tutorial_creative_model.h5m")
creative_model.write_file(output_h5m_path)
print(f"Creative model saved to: {output_h5m_path}")

# Export the cube surface (ID 201 in creative_model) to VTK
if created_surf_id in creative_model.surfaces_by_id:
    cube_surf_to_export = creative_model.surfaces_by_id[created_surf_id]
    output_vtk_path_str = "cube_surface_export_from_creative_model" # .vtk suffix added automatically if missing
    cube_surf_to_export.to_vtk(output_vtk_path_str)
    print(f"Surface {cube_surf_to_export.id} exported to: {output_vtk_path_str}.vtk")

    # Export one of the volumes (ID 101 in creative_model) to VTK
    if created_vol_id in creative_model.volumes_by_id:
      vol_to_export = creative_model.volumes_by_id[created_vol_id]
      output_vol_vtk_str = "created_volume_export"
      vol_to_export.to_vtk(output_vol_vtk_str)
      print(f"Volume {vol_to_export.id} exported to: {output_vol_vtk_str}.vtk")

Creative model saved to: pydagmc_tutorial_creative_model.h5m
Surface 201 exported to: cube_surface_export_from_creative_model.vtk
Volume 101 exported to: created_volume_export.vtk

9. Deleting Entities#

Entities can be removed from the model using their delete() method. This action removes the entity from the underlying MOAB database and updates pydagmc’s internal tracking. The Python object representing the deleted entity will no longer be valid for model operations.

[14]:

# Using the 'creative_model'
# Let's delete 'another_my_volume' (ID was likely 102 in creative_model)
vol_id_to_delete_from_creative = another_my_volume.id # from section 6

if vol_id_to_delete_from_creative in creative_model.volumes_by_id:
    vol_object_to_delete = creative_model.volumes_by_id[vol_id_to_delete_from_creative]
    print(f"Preparing to delete Volume ID {vol_object_to_delete.id} from creative_model.")
    assert vol_id_to_delete_from_creative in creative_model.used_ids[pydagmc.Volume]

    vol_object_to_delete.delete()
    print(f"Volume with original ID {vol_id_to_delete_from_creative} deleted.")

    assert vol_id_to_delete_from_creative not in creative_model.volumes_by_id
    assert vol_id_to_delete_from_creative not in creative_model.used_ids[pydagmc.Volume]

    # After deletion, the object's .handle and .model attributes are None.
    # Accessing properties that rely on the model will raise an error.
    print(f"Deleted object's handle: {vol_object_to_delete.handle}") # Expected: None
    try:
        _ = vol_object_to_delete.category
    except AttributeError as e:
        print(f"Error accessing property of deleted volume (expected): {e}")
else:
    print(f"Volume with ID {vol_id_to_delete_from_creative} not found in creative_model for deletion.")

Preparing to delete Volume ID 102 from creative_model.
Volume with original ID 102 deleted.
Deleted object's handle: None
Error accessing property of deleted volume (expected): 'NoneType' object has no attribute 'category_tag'

10. A Look Under the Hood: How Key Operations Work#

This section provides a brief, high-level overview of the internal PyMOAB operations that pydagmc uses to perform common tasks. This is not information you need to use the library, but it may be helpful for understanding its behavior or for debugging.

  • Loading a Model (``Model(filename)``):

    • Initializes pymoab.core.Core().

    • Calls mb.load_file(filename).

    • Pre-populates used_ids by querying existing entities.

  • Accessing Entities (e.g., ``model.volumes_by_id``):

    • Uses mb.get_entities_by_type_and_tag() with the root set and the CATEGORY_TAG to find all entity sets of a specific type (e.g., “Volume”).

    • For each handle found, a PyDAGMC object (e.g., Volume(model, handle)) is instantiated.

    • These are often collected into dictionaries keyed by ID for quick access.

  • Triangle Data (``.triangle_coords``, ``.triangle_conn``):

    • Surface._get_triangle_sets() returns itself.

    • Volume._get_triangle_sets() iterates through its child surfaces and collects their handles.

    • Group._get_triangle_sets() iterates through its child surfaces and volumes (recursively calling their _get_triangle_sets()).

    • GeometrySet.triangle_handles then uses mb.get_entities_by_type(handle, types.MBTRI) on these collected surface/set handles.

    • mb.get_connectivity() and mb.get_coords() are then used on the triangle handles.

  • Material Assignment (``Volume.material``):

    • The Volume.groups property finds all parent groups of the volume.

    • It then iterates through these groups, checking if their group.name starts with "mat:". The first one found determines the material.

    • Setting Volume.material = "new_mat" will:

      • Remove the volume from its current material group (if any).

      • Find or create a Group named "mat:new_mat".

      • Add the volume to this new/existing material group.

  • Entity Creation (e.g., ``Model.create_volume()``):

    • Calls mb.create_meshset() to get a new entity set handle.

    • Wraps this handle in a temporary GeometrySet to set the default CATEGORY_TAG and GEOM_DIMENSION_TAG appropriate for the entity type (e.g., “Volume”, 3).

    • Then, the specific class (e.g., Volume) is instantiated with this configured handle.

    • An ID is assigned, either user-provided or auto-generated to be unique within the model for that entity type.

  • Volume Calculation (``Volume.volume``): The volume of a Volume object is calculated using the formula derived from the divergence theorem:

    \[V = \frac{1}{3} \sum_{i \in \text{surfaces}} \text{sign}_i \sum_{j \in \text{triangles in surface } i} (\mathbf{a}_j \times \mathbf{b}_j) \cdot \mathbf{c}_j / 2\]

    This simplifies to:

    \[V = \frac{1}{6} \sum_{s \in \text{Surfaces}} \text{sense}(s,V) \sum_{t \in \text{Triangles in } s} \det(\mathbf{v}_{t1}, \mathbf{v}_{t2}, \mathbf{v}_{t3})\]

    where \(\text{sense}(s,V)\) is +1 if \(V\) is the forward volume of surface \(s\), and -1 if it’s the reverse volume. \(\mathbf{v}_{t1}, \mathbf{v}_{t2}, \mathbf{v}_{t3}\) are the vertex coordinates of triangle \(t\).

    PyDAGMC iterates through each Surface of the Volume. For each triangle in the surface, it computes \(\sum (\mathbf{r}_0 \cdot (\mathbf{r}_1 - \mathbf{r}_0) \times (\mathbf{r}_2 - \mathbf{r}_0))\). The sum over all triangles is then scaled by \(1/6\). The sign depends on whether the current Volume is the forward or reverse volume of the Surface.

11. Conclusion#

This tutorial has introduced the core functionalities of pydagmc for interacting with DAGMC geometric models. You’ve learned how to:

  • Load and get an overview of models.

  • Access and query volumes, surfaces, and groups and their properties.

  • Modify existing geometry by changing IDs, materials, surface senses, and group memberships.

  • Create new geometric entities and structure them into a model.

  • Retrieve detailed triangular mesh data.

  • Save models and export parts for visualization.

  • Remove entities from a model.

pydagmc offers a convenient Python API over PyMOAB for DAGMC-specific operations, facilitating scripting for geometry preprocessing, analysis, and modification tasks. For further details, consult the pydagmc source code and the PyMOAB documentation.

[15]:

# Clean up downloaded and created files at the end of the session
print("Cleaning up tutorial files...")
if 'FUEL_PIN_FILE' in locals() and FUEL_PIN_FILE.exists(): FUEL_PIN_FILE.unlink()
if 'CUBE_STL_FILE' in locals() and CUBE_STL_FILE.exists(): CUBE_STL_FILE.unlink()
if 'output_h5m_path' in locals() and output_h5m_path.exists(): output_h5m_path.unlink()
if 'output_vtk_path_str' in locals() and Path(f"{output_vtk_path_str}.vtk").exists(): Path(f"{output_vtk_path_str}.vtk").unlink()
if 'output_vol_vtk_str' in locals() and Path(f"{output_vol_vtk_str}.vtk").exists(): Path(f"{output_vol_vtk_str}.vtk").unlink()
print("Cleanup attempts complete.")

Cleaning up tutorial files...
Cleanup attempts complete.
Contents