Mesh integration #10

eitan-weinstein · 2023-07-27T16:22:08Z

This PR integrates tetrahedral meshing into the magnet coil creation, so that it can be directly incorporated into the parametric stellarator, if desired. ExampleScript.py was also modified to allow for the meshing.

…roid script

gonuke

A few initial thoughts... the integration with the @connoramoreno code is good. For the mesh analysis code... I've suggested some things to break up the code.

gonuke · 2023-07-27T23:00:10Z

parametric_stellarator.py

+        h5m_path = f'{cwd}/coil_mesh.h5m'
+        h5m_conversion = subprocess.run(f'mbconvert {exo_path} {h5m_path}', shell = True)


I think we can do this with pymoab commands by just reading and writing a file in different formats, possibly in a new MOAB instance but maybe just into a new meshset? It will be lighter weight than launching a subprocess.

I looked in the MOAB bitbucket for a pymoab version of mbconvert, but could not find anything. Is there a specific function that you have in mind to replace the subprocess command?

Something like:

myFileset = pymoab.EntitySet() myMoab.read_file('file.exo', fileset = myFileset) myMoab.write_file('file.h5m', fileset = myFileset)

gonuke · 2023-07-27T23:03:08Z

coils_tet_mesh.py

+    cubit.cmd(f'import step "{coils_path}" heal')
+    cubit.cmd('volume  size auto factor 5')
+
+    # For-loop iterates over all 48 coils included in the coils.step file


Comment still refers to 48 coils

gonuke · 2023-07-27T23:03:45Z

coils_tet_mesh.py

+
+        # Create the mesh
+        cubit.cmd(f'volume {coil_num} scheme tetmesh')
+        cubit.cmd(f'volume {coil_num} scheme tetmesh')


Repeated line... does it do anything?

gonuke · 2023-07-27T23:05:44Z

coils_tet_mesh.py

+        cubit.cmd(f'mesh volume {coil_num}')
+
+        # Set log file path and initiate logging
+        log_split = log_general_path.split('.')


Do this once before the loop, perhaps name the components for clarity something like log_prefix and log_ext?

gonuke · 2023-07-27T23:06:36Z

coils_tet_mesh.py

+        # Mesh quality metrics
+        quality_metric_list = ['shape', 'aspect ratio', 'condition no.', 'distortion', 'element volume',
+                                'equiangle skew', 'equivolume skew','inradius', 'jacobian', 
+                                'normalized inradius', 'relative size', 'scaled jacobian', 'shape and size']


Define this once before the loop

gonuke · 2023-07-27T23:09:20Z

coils_tet_mesh.py

+        # Read the text file
+        with open(log_path, 'r') as file:
+            text = file.read()
+
+        # Cut out each of the string 'tables' 
+        start_table = 'Function Name'
+        end_table = 'Finished Command:'
+        ranges = []
+        start_index = 0
+
+        while True:
+            # Find the starting index after the previous occurrence
+            start_index = text.find(start_table, start_index)
+
+            # Break the loop if the starting cue is not found
+            if start_index == -1:
+                break
+
+            # Find the ending index after the starting index
+            end_index = text.find(end_table, start_index + 1)
+
+            # Break the loop if the ending cue is not found
+            if end_index == -1:
+                break
+
+            # Extract the range between the starting and ending index
+            range_text = text[start_index : end_index + 1]
+            ranges.append(range_text)
+
+            # Update the starting index for the next iteration
+            start_index = end_index + 1


Maybe this can be a function to reduce cognitive load

Arguably, the file reading should perhaps be outside the function for separation of concerns

gonuke · 2023-07-27T23:12:32Z

coils_tet_mesh.py

+            range_split = qual_metric[qual_metric.index('-\n'):].split(' ')
+
+            row_extraction = []
+            for i in range_split:
+                if i != '' and "(" not in i:
+                    row_extraction.append(i)
+            offset = len(row_extraction) - 4
+            row_extraction[0] = [''.join(row_extraction[:offset])]
+
+            for j in np.arange(1,5):
+                row_extraction[j] = row_extraction[offset + (j-1) : offset + j]
+            row_extraction = row_extraction[:5]
+            row_extraction[0][0] = row_extraction[0][0].split('-\n')[1]
+
+            data_lists = [function_names, averages, std_devs, minima, maxima]
+            indices = np.arange(5)
+            for index, data_list in zip(indices, data_lists):
+                if index == 0:
+                    data_list.append(row_extraction[index][0])
+                else:
+                    data_list.append(float(row_extraction[index][0]))


Maybe another function to reduce cognitive load

We'll also want to review whether there are better ways to do this parsing....

gonuke · 2023-07-27T23:14:09Z

coils_tet_mesh.py

+        # Recreate the geometry and mesh with the best_analytic
+        cubit.cmd('reset')
+        cubit.cmd(f'import step "{coils_path}" heal')
+        cubit.cmd('volume  size auto factor 5')
+        cubit.cmd(f'volume {coil_num} scheme tetmesh')
+        cubit.cmd(f'volume {coil_num} scheme tetmesh')
+        cubit.cmd(f'mesh volume {coil_num}')
+        cubit.cmd(f'mesh volume {coil_num}')
+        cubit.cmd(f'quality volume {coil_num} {best_analytic} global draw mesh list')


I'm not sure this is doing anything meaningful...

eitan-weinstein · 2023-07-28T16:51:35Z

Thanks for the comments, Paul. I'll get on them now and see if I can implement these improvements.

…larator" This reverts commit a67859b.

eitan-weinstein · 2023-08-07T21:28:05Z

Sorry about the whirlwind of commits in the past few minutes. I had made an accidental commit when trying to update the coils_tet_mesh.py script and then made a few compounding mistakes on GitHub trying to fix it all. It looks like it should be updated properly now. My newest change incorporates the PyMOAB commands in the mesh file conversion, as well as the export of a .csv file containing the mesh quality analytics for each coil.

gonuke · 2023-08-08T01:20:13Z

somewhere you have lost the integration of meshing within the main package

eitan-weinstein · 2023-08-08T15:48:58Z

Alright, the integration has been re-incorporated and updated to include the PyMOAB commands.

gonuke

Looks great! The integrated meshing is great, still some more refining on the stand-alone analysis.

gonuke · 2023-08-08T21:05:20Z

parametric_stellarator.py

+            # Get current working directory
+            cwd = os.getcwd()      
+            # Exodus export
+            exo_path = f'{cwd}/coil_mesh.exo'


Define the base_name of the file

base_name = 'coil_mesh' exo_path = f'{file_path}/{base_name}.exo'

and use that in both places. This will facilitate changing it in the future.

gonuke · 2023-08-08T21:06:47Z

parametric_stellarator.py

+        # Conditionally export tetrahedral meshing
+        if magnets['meshing']:
+            # Get current working directory
+            cwd = os.getcwd()      


Suggested change

cwd = os.getcwd()

file_path = os.getcwd()

gonuke · 2023-08-08T21:07:37Z

parametric_stellarator.py

+            # Create a new meshset for each coil iteration
+            coil_meshset = mb.create_meshset()
+            # Add the current coil's mesh to the meshset
+            mb.add_entities(coil_meshset, [exodus_set])
+            # Set .h5m path
+            h5m_path = f'{cwd}/coil_mesh.h5m'
+            # Write the current coil's meshset to an individual .h5m file
+            mb.write_file(h5m_path, [coil_meshset])


Suggested change

# Create a new meshset for each coil iteration

coil_meshset = mb.create_meshset()

# Add the current coil's mesh to the meshset

mb.add_entities(coil_meshset, [exodus_set])

# Set .h5m path

h5m_path = f'{cwd}/coil_mesh.h5m'

# Write the current coil's meshset to an individual .h5m file

mb.write_file(h5m_path, [coil_meshset])

# Set .h5m path

h5m_path = f'{file_path}/{base_name}.h5m'

# Write the current coil's meshset to an individual .h5m file

mb.write_file(h5m_path, [exodus_set])

coils_tet_mesh.py

gonuke · 2023-08-08T21:15:45Z

coils_tet_mesh.py

+    log_split = log_general_path.split('.')
+    log_prefix = log_split[0]
+    log_ext = log_split[1]


might be easier with the pathlib package/module??

gonuke · 2023-08-08T21:20:17Z

coils_tet_mesh.py

+        # Override the current mesh analytic with the best quality metric
+        cubit.cmd(f'quality volume {coil_num} {best_analytic} global draw mesh list')


May not be necessary?? We have this in the log file and/or data frame already

gonuke · 2023-08-08T21:21:35Z

coils_tet_mesh.py

+    # For-loop iterates over all of the coils included in the coils.step file
+    for coil_num in np.arange(1,total_coil_count + 1):
+
+        # Create the mesh
+        cubit.cmd(f'volume {coil_num} scheme tetmesh')
+        cubit.cmd(f'mesh volume {coil_num}')
+
+        # Set log file path and initiate logging
+        log_path = f"{log_prefix}_{coil_num}.{log_ext}"
+        cubit.cmd(f'logging on file "{log_path}"')
+
+        # Run mesh quality metrics for each relevant metric
+        for metric in quality_metric_list:
+            qual_command = f'quality volume {coil_num} {metric} global draw mesh list'
+            mesh_analytics = cubit.cmd(qual_command)


Separation of concerns & cognitive load:
Make this a function that creates and does the quality measures on each mesh

gonuke · 2023-08-08T21:22:34Z

coils_tet_mesh.py

+        # MESH QUALITY ANALYSIS
+
+        # Read the log file
+        with open(log_path, 'r') as file:
+            text = file.read()
+
+        # Extract the data from the log file using the extract_tabular_data function
+        function_names, averages, std_devs, minima, maxima = extract_tabular_data(text)
+
+        # Separate multi-word function names to be able to match-case with Cubit names
+        function_names = split_function_names(function_names)
+
+        # Compile the extracted tabular data into a Pandas data frame
+        quality_metrics_df = pd.DataFrame({'Function Name': function_names,
+                                           'Average': averages,
+                                           'Standard Deviation': std_devs,
+                                           'Minimum': minima,
+                                           'Maximum': maxima})
+
+        # Remove rows where the standard deviation equals 0
+        # Found this for a few cases and it appeared too implausible to keep
+        quality_metrics_df = quality_metrics_df.loc[quality_metrics_df['Standard Deviation'] != 0]


Separation of concerns & cognitive load:
Make this a method that builds & cleans the dataframe

gonuke · 2023-08-08T21:23:32Z

coils_tet_mesh.py

+    # MESH EXPORT
+
+    # Export the mesh as an EXODUS file
+    cubit.cmd(f'export mesh "{exo_path}"')
+
+    # Initialize the MOAB core instance
+    mb = core.Core()
+
+    # Load the EXODUS file
+    exodus_set = mb.create_meshset()
+    mb.load_file(exo_path, exodus_set)
+
+    # Create a new meshset for each coil iteration
+    coil_meshset = mb.create_meshset()
+
+    # Add the current coil's mesh to the meshset
+    mb.add_entities(coil_meshset, [exodus_set])
+
+    # Write the current coil's meshset to an individual .h5m file
+    mb.write_file(h5m_path, [coil_meshset])


Method to export the mesh - make it optional based on the user input

Also update as in the parametric stellarator code

Consider a method that can be reused here and in parametric_stellarator for mesh export as H5M.

gonuke · 2023-08-08T21:24:19Z

coils_tet_mesh.py

+    # CSV EXPORT
+
+    # Rename 'index' column to 'Coil Number'
+    top_results.reset_index(inplace = True)
+    top_results.rename(columns = {'index' : 'Coil Number'}, inplace = True)
+
+    # Fix indexing for coil numbers
+    top_results = top_results.transpose()
+    top_results.iloc[0] = np.arange(1, total_coil_count + 1)
+    top_results = top_results.transpose()
+    print(top_results.head())
+
+    # Export the top_results data frame as a .csv
+    top_results.to_csv(csv_path)


Function to write CSV - optional?
Separation of concerns: separate analysis from I/O

gonuke

Some quick feedback - very much in the fine tuning style aspects

gonuke · 2023-08-09T18:58:49Z

parametric_stellarator.py

+        # Conditionally export tetrahedral meshing
+        if magnets['meshing']:
+            # Assign desired directory
+            file_path = os.getcwd()
+            # Define the base name of the file
+            base_name = 'coil_mesh'
+            # Exodus export
+            exo_path = f'{file_path}/{base_name}.exo'
+            cubit.cmd(f'export mesh "{exo_path}"')
+            # Initialize the MOAB core instance
+            mb = core.Core()
+            # Load the EXODUS file
+            exodus_set = mb.create_meshset()
+            mb.load_file(exo_path, exodus_set)
+            # Set .h5m path
+            h5m_path = f'{file_path}/{base_name}.h5m'
+            # Write the current coil's meshset to an individual .h5m file
+            mb.write_file(h5m_path, [exodus_set])


A couple of fine tuning thoughts related to best practice, readability and cognitive load:

Consider whether all the comments are necessary; good choice of variable names and code structure can make comments redundant and interrupt the flow reading a snippet of code. Comments can often be reserved for explain it more complicated algorithmic choices.

If the two filename definitions are adjacent, it can be more readable by their association to each other as well as removing them from the "action" flow of the read/write steps

gonuke

A few more ways to streamline and simplify some of the algorithms.

I think we perhaps should have sat down to discuss the cleanest way to extract data from the logs - it seems far more complicated than it probably needs to be.

There may be some value in separating the coils_tet_mesh.py script from the rest of this PR.

gonuke · 2023-08-14T17:22:13Z

coils_tet_mesh.py

+# Create an argument parser
+parser = argparse.ArgumentParser(description="Magnet Coils Tetrahedral Mesh Script")
+
+# Define command-line arguments
+parser.add_argument("coils_path", type=str, help="Path to the coils .step file")
+parser.add_argument("total_coil_count", type=int, help="Number of coils")
+parser.add_argument("--cwd", action="store_true", help="Export outputs to current working directory")
+parser.add_argument("--mesh-export", action="store_true", help="Export meshes as .exo and .h5m files")
+parser.add_argument("--log-path", type=str, help="Path to .log file prefix")
+parser.add_argument("--exo-path", type=str, help="Path to .exo file")
+parser.add_argument("--h5m-path", type=str, help="Path to .h5m file")
+parser.add_argument("--csv-export", action="store_true", help="Export data as a .csv file")
+parser.add_argument("--csv-path", type=str, help="Path to .csv file")
+
+# Parse the command-line arguments
+args = parser.parse_args()
+


This needs to be buried in a method as well.

One goal is to allow this to be imported as a module into another python script. Ideally, no code should run at the time of import - or only code that is necessary at the time of import. As it is, this parser code will run when the module is imported and may cause issues.

That method can then either be called from within the conditional below, or from within coils_tet_mesh directly

gonuke · 2023-08-14T17:31:44Z

coils_tet_mesh.py

+# Define command-line arguments
+parser.add_argument("coils_path", type=str, help="Path to the coils .step file")
+parser.add_argument("total_coil_count", type=int, help="Number of coils")
+parser.add_argument("--cwd", action="store_true", help="Export outputs to current working directory")


For most tools, this is the default, and then the other override it. You've created an unusual hybrid where the other options define the output locations and this option overrides those. For the sake of moving forward, I'd recommend removing many of these options - we can add them back in as needed.

gonuke · 2023-08-14T19:59:55Z

coils_tet_mesh.py

+        for char in name:
+            # Check if the character is an uppercase letter
+            if char.isupper() and current_word.isupper() and len(current_word) > 1:
+                # If the character is uppercase and the current_word is also uppercase
+                # and not a single character, consider it a separate word and append
+                # the current_word to the split_words list.
+                split_words.append(current_word)
+                current_word = char
+            elif char.isupper() and current_word and not current_word.isupper():
+                # If the character is uppercase and the current_word exists but is not
+                # uppercase, consider it a separate word and append the current_word
+                # to the split_words list.
+                split_words.append(current_word)
+                current_word = char
+            else:
+                # If the character is not uppercase or doesn't meet the above conditions,
+                # append it to the current_word to build the word.
+                current_word += char
+
+        # Append the last word to the split_words list
+        split_words.append(current_word)
+


A slightly simpler approach might be to insert a special character (underscore?) before each upper case and then use the python split()

Suggested change

for char in name:

# Check if the character is an uppercase letter

if char.isupper() and current_word.isupper() and len(current_word) > 1:

# If the character is uppercase and the current_word is also uppercase

# and not a single character, consider it a separate word and append

# the current_word to the split_words list.

split_words.append(current_word)

current_word = char

elif char.isupper() and current_word and not current_word.isupper():

# If the character is uppercase and the current_word exists but is not

# uppercase, consider it a separate word and append the current_word

# to the split_words list.

split_words.append(current_word)

current_word = char

else:

# If the character is not uppercase or doesn't meet the above conditions,

# append it to the current_word to build the word.

current_word += char

# Append the last word to the split_words list

split_words.append(current_word)

name_split = ''

delim = '_'

for char in name:

# insert a delim character before each upper case letter to facilitate splitting

if char.isupper():

name_split += delim

name_split += char

split_words = name_split.split(delim)

gonuke · 2023-08-15T13:36:01Z

coils_tet_mesh.py

+        indices = np.arange(5)
+        for index, data_list in zip(indices, data_lists):
+            if index == 0:
+                data_list.append(row_extraction[index][0])
+            else:
+                data_list.append(float(row_extraction[index][0]))


Suggested change

indices = np.arange(5)

for index, data_list in zip(indices, data_lists):

if index == 0:

data_list.append(row_extraction[index][0])

else:

data_list.append(float(row_extraction[index][0]))

for index, data_list in enumerate(data_lists):

if index == 0:

data_list.append(row_extraction[index][0])

else:

data_list.append(float(row_extraction[index][0]))

eitan-weinstein added 3 commits July 17, 2023 13:11

Adding toroid and coil tet mesh scripts

aea3373

Incorporating tet meshing for the coils as a function and removing to…

046a1d5

…roid script

Integrating tet meshing into magnet coils and parametric stellarator

a67859b

gonuke requested changes Jul 27, 2023

View reviewed changes

eitan-weinstein added 4 commits August 7, 2023 16:02

drop

2af2411

Incorporating PyMOAB Functionality and .csv Export

087b458

Mistake in last commit. PyMOAB and .csv

88dcfdd

Revert "Integrating tet meshing into magnet coils and parametric stel…

f00a903

…larator" This reverts commit a67859b.

eitan-weinstein force-pushed the mesh_integration branch from c4757b4 to f00a903 Compare August 7, 2023 21:23

eitan-weinstein added 3 commits August 8, 2023 10:34

Putting back integration into ExampleScript.py

048d516

Putting back integration into magnet_coils.py

d607d6a

Putting back integration into parametric_stellarator.py

3bb1afc

gonuke requested changes Aug 8, 2023

View reviewed changes

eitan-weinstein added 2 commits August 9, 2023 10:12

Tweaking file path and simplifying PyMOAB meshsets

694a58e

Missing apostrophe in base_name string

37f6b5a

gonuke requested changes Aug 9, 2023

View reviewed changes

eitan-weinstein and others added 3 commits August 9, 2023 16:53

Separated out each function and enabled command-line execution

59bb7f3

Improving readability of meshing export

cf75ab9

Including argparser

45d802b

gonuke requested changes Aug 15, 2023

View reviewed changes

gonuke merged commit 45d802b into svalinn:main Oct 2, 2023

		h5m_path = f'{cwd}/coil_mesh.h5m'
		h5m_conversion = subprocess.run(f'mbconvert {exo_path} {h5m_path}', shell = True)

		# Override the current mesh analytic with the best quality metric
		cubit.cmd(f'quality volume {coil_num} {best_analytic} global draw mesh list')

Mesh integration #10

Mesh integration #10

Conversation

eitan-weinstein commented Jul 27, 2023 • edited Loading

gonuke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eitan-weinstein commented Jul 28, 2023

eitan-weinstein commented Aug 7, 2023

gonuke commented Aug 8, 2023

eitan-weinstein commented Aug 8, 2023

gonuke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gonuke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gonuke left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eitan-weinstein commented Jul 27, 2023 •

edited

Loading