Python API Scripts

JohnM

I'm starting this as a "fresh" Forum topic in the hopes it can become a repository for good API scripts/tricks. I'll add a few of our own here, hopefully others will follow!

JohnM

Python script to create center node of circle between 3 points.
I just adapted this from old FORTRAN (no jokes), and it "seems" to work. I can't tell you what it is doing with more than 3 nodes selected - I tested a few examples and it wasn't clear. But if 3 nodes are selected, it seems to work.

If anyone improves on this (it needs error checking) please re-post.


# Get the coordinates of the points that define the circle from 3 selected nodes.
# Given three points in space (A,B,C)perimeter_node_ids = mw.selected_nodes()
perimeter_node_ids = mw.selected_nodes()
point_1 = mw.node(perimeter_node_ids[0])
point_2 = mw.node(perimeter_node_ids[1])
point_3 = mw.node(perimeter_node_ids[2])

# ... do some math but not this ...
Xa=point_1.x
Ya=point_1.y
Za=point_1.z

Xb=point_2.x
Yb=point_2.y
Zb=point_2.z

Xc=point_3.x
Yc=point_3.y
Zc=point_3.z

# Lengths of AB, AC, AC
AB = (((Xa-Xb)**2) + ((Ya-Yb)**2) + ((Za-Zb)**2))**0.5
BC = (((Xb-Xc)**2) + ((Yb-Yc)**2) + ((Zb-Zc)**2))**0.5
AC = (((Xa-Xc)**2) + ((Ya-Yc)**2) + ((Za-Zc)**2))**0.5
# Direction cosines of AB(ABi,ABj,ABk)
ABi = (Xb-Xa) / AB
ABj = (Yb-Ya) / AB
ABk = (Zb-Za) / AB
# Direction cosines of AC(ACi,ACj,ACk)
ACi = (Xc-Xa) / AC
ACj = (Yc-Ya) / AC
ACk = (Zc-Za) / AC
# Cosine of angle BAC
cosBAC = (AB**2 + AC**2-BC**2) / (2 * AB * AC)
AD = cosBAC * AC
CD = (AC**2-AD**2)**0.5
# Position of point D, which is C projected normally onto AB
Xd = Xa + (AD * ABi)
Yd = Ya + (AD * ABj)
Zd = Za + (AD * ABk)
# Direction cosines of CD(Cdi,CDj,CDk)
CDi = (Xc-Xd) / CD
CDj = (Yc-Yd) / CD
CDk = (Zc-Zd) / CD
# Direction cosines of normal to AB and CD
# to be used for rotations of circle center
Ni = (ABk * CDj)-(ABj * CDk)
Nj = (ABi * CDk)-(ABk * CDi)
Nk = (ABj * CDi)-(ABi * CDj)
# # Diameter of circumscribed circle of a triangle is equal to the
# the length of any side divided by sine of the opposite angle.
# This is done in a coordinate system where X is colinear with AB, Y is // to CD,
# and Z is the normal (N) to X and Y, and the origin is point A
# R = D / 2
sinBAC = (1-cosBAC**2)**0.5
R = (BC / sinBAC) / 2
# center of circumscribed circle is point E
X2e = AB / 2
Y2e = (R**2-X2e**2)**0.5
Z2e = 0
# Transform matrix
# Rotations Translations
#
# ABi , ABj , ABk Xa
# CDi , CDj , CDk Ya
# Ni , Nj , Nk Za
#
# Position of circle center in absolute axis system
X_center = Xa + (X2e * ABi) + (Y2e * CDi) + (Z2e * Ni)
Y_center = Ya + (X2e * ABj) + (Y2e * CDj) + (Z2e * Nj)
Z_center = Za + (X2e * ABk) + (Y2e * CDk) + (Z2e * Nk)
#
center = Vector(X_center,Y_center,Z_center)

# Create the center node
mw.new_node(center)

JohnM

Modify selected nodes to the same constant X value
(for y and z edit the "x" lines below

def mod_node(nodes, value):
for point in nodes:
mw.set_node_x(point, value)


active_nodes = mw.selected_nodes() # Get selected nodes
if len(active_nodes) > 0:
value = mw.input("Input node x value(m): ")
mod_node(active_nodes, float(value))
else:
mw.message("No nodes selected.")

JohnM

Solve then Save

mw.solve()

mw.save()

JohnM

Transfer displacements from solution (one of Victor's examples)

# This example script replicates the functionality of
#
# Mesh tools -> Transfer displacements from solution
#
# but is limited to solutions without time steps or modes and lacks error
# checking.

assert mw.version() == 13

for node_id in solution.all_nodes():

# Read displacement from solution
displacement_x = solution.node_value("displx", node_id)
displacement_y = solution.node_value("disply", node_id)
displacement_z = solution.node_value("displz", node_id)

# Calculate new node position
new_x = mw.node(node_id).x + displacement_x
new_y = mw.node(node_id).y + displacement_y
new_z = mw.node(node_id).z + displacement_z

# Change node position
mw.set_node_x(node_id, new_x)
mw.set_node_y(node_id, new_y)
mw.set_node_z(node_id, new_z)

mw.message("Finished")

CornelisP

I found myself frequently needing to get the forces on some nodes/surface that I hadn't defined before solving the model, so I put together a quick script to do that. One thing I wasn't able to do is check the user's current time step so I can pull just that force. Instead, it spits out the force for every time step and the user has to figure out which one they want. If anyone knows a way to do that, it'd be appreciated.

Sum forces over a given set of nodes


# This script sums the force in a given direction on a user-selected set of nodes options = ["x", "y", "z"] force_dir = mw.input("Choose direction (X, Y, or Z)") force_dir = force_dir.lower() if force_dir == "": mw.message("No direction given!") elif force_dir not in options: mw.message("X, Y, or Z are the only valid directions") else: nodes = solution.selected_nodes() if not nodes: mw.message("No nodes selected! Make sure you didn't select faces or elements instead.") else: force = 0 if solution.time_steps(): steps = solution.time_steps() mess = "" for step in steps: for node in nodes: force = force + solution.node_value("externalforce" + force_dir, node, step) mess = mess + "At " + str(solution.time(step)) + " sec, the sum of forces in " + force_dir +\ " = " + str(force) + " N\n" force = 0 mess = mess[:-1] # Remove extra newline mw.message(mess) else: for node in nodes: force = force + solution.node_value("externalforce" + force_dir, node) mw.message("Sum of forces in " + force_dir + " = " + str(force) + " N")

Hengre

Calling an external program via phyton script

The following lines shows how to call an external program (EXE) from within Mecway using a phyton script.
The BS.EXE is located in a chosen directory,
a input file say "a.inp" (here some mesh data created with Mecway) can belocated in a different directory. The file name and the parameter (20.) is passed to the Executable.

import subprocess
subprocess.call("C:/temp/BIOT/bin/bs.exe C:/temp/BIOT/a.inp 20.")

Following Victors recomment I used in the path spezification the forward slash "/", but
double slash works too. e.g.
subprocess.call("C:\\temp\\BIOT\\bin\\bs.exe C:\\temp\\BIOT\\a.inp 20.")

Regards

Victor

Skew tool. Simple script to deform a mesh like this:

struthon

Hi.
There is "Formula" features that make possible to add new results for Solution. Using this metod has some limitation so it is not possible to make advanced result processing. If it is possible to do it using API script feature? I trying to make a script that will calculate required reinforcement for shell elements and display it as colour maps in the same way as other results are displayed. If the answer for above is "No" the next question is if it is possible to inject additional results directly to .liml file? Thanks!

JohnM

Attached below is get_some_stresses.py , a nice script written by @CornelisP that outputs 6 component stresses to a text file for further postprocessing. Also attached is text_format_solution.liml, which is a sample model with post results written in non-binary format. With a little
scripting you can output your required data, create the results you want, then reformat as a Mecway post output using the sample results as a template.

After you have made something amazing, by all means share

Victor

@struthon. No, the API doesn't allow that but will do in the next release (v15). Yes, you can add data to the .liml file. It accepts both a binary format and a text format which you can work out by looking at files saved by older versions up to v3 https://mecway.com/download/oldversions/mecway30.msi . Or you can put the data in a CCX .frd or Gmsh .msh file which Mecway can open.

struthon

Thanks for those help, looks it will helps me make first tests. Will share my experience in case of success. Do we already know when the v15 is expected?

JohnM

Instantly like this one by Victor that outputs element volumes to a file:


import os
filename = os.path.dirname(mw.file_name()) + "\\volumes.txt"

with open(filename, 'w') as file:
for element_id in solution.all_elements():
file.write(str(solution.volume(element_id)) + "\n")

JohnM

EDITED BY VICTOR. This code is now obsolete. See
https://mecway.com/forum/discussion/comment/6503/#Comment_6503 for a simpler way using solution.interpolate().


Results extraction at a given coordinate without node.
from Victor
12:53AM edited 12:55AM
You can do the interpolation with a script but it's a bit of work. Here's an example for hex8 elements. For other shapes, you would need to use the appropriate shape functions.

To use the script, specify an element and a point inside the element using its isoparametric coordinates r, s, t. If you run it as is with the example model, it uses (0,-1,0) which is in the center of the face that has 222 Pa stress.

def shape_function_hex8(element_node, r, s, t):

  # r,s,t are each in the range [-1,1] with 0,0,0 in the middle of the
  # element. They are in the same directions as the R,S,T used in
  # Mesh Tools -> Refine -> Custom.

  if element_node == 1:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * -1) * (1 + t * -1)
  if element_node == 2:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * -1) * (1 + t * -1)
  if element_node == 3:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * 1) * (1 + t * -1)
  if element_node == 4:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * 1) * (1 + t * -1)
  if element_node == 5:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * -1) * (1 + t * 1)
  if element_node == 6:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * -1) * (1 + t * 1)
  if element_node == 7:
    return 1.0 / 8 * (1 + r * 1) * (1 + s * 1) * (1 + t * 1)
  if element_node == 8:
    return 1.0 / 8 * (1 + r * -1) * (1 + s * 1) * (1 + t * 1)


def interpolate(element_id, variable, r, s, t):
  # Only for hex8 elements.

  sum = 0
  nodes = solution.nodes(element_id)
  for element_node in range(1, len(nodes) + 1):
    shape_function = shape_function_hex8(element_node, r, s, t)
    node_id = nodes[element_node - 1]
    value = solution.node_value(variable, node_id)
    sum += shape_function * value

  return sum


assert mw.version() >= 15, "Wrong Mecway version"

# Enter parameters here
value = interpolate(element_id = 1,
                    variable = "stressxx",
                    r = 0,
                    s = -1,
                    t = 0)

mw.message(str(value))

https://mecway.com/forum/discussion/1121/results-extraction-at-a-given-coordinate-without-node#latest

Victor

Better use the "code" formatting for posting code on the forum otherwise the indenting and possibly some symbols get broken. @JohnM, I've edited your post to apply code formatting.

Sergio

Based on the code of @JohnM:

Info window with diameter, radius and center of a circle by 3 points (just info, no node creation)

@Victor, is there a way to know which length units the user is using? I have hardcoded the results in mm, but would be nice to inform in user units.

# Get the coordinates of the points that define the circle from 3 selected nodes.
# Given three points in space (A,B,C)perimeter_node_ids = mw.selected_nodes()
perimeter_node_ids = mw.selected_nodes()
point_1 = mw.node(perimeter_node_ids[0])
point_2 = mw.node(perimeter_node_ids[1])
point_3 = mw.node(perimeter_node_ids[2])

# ... do some math but not this ...
Xa=point_1.x
Ya=point_1.y
Za=point_1.z

Xb=point_2.x
Yb=point_2.y
Zb=point_2.z

Xc=point_3.x
Yc=point_3.y
Zc=point_3.z

# Lengths of AB, AC, AC
AB = (((Xa-Xb)**2) + ((Ya-Yb)**2) + ((Za-Zb)**2))**0.5
BC = (((Xb-Xc)**2) + ((Yb-Yc)**2) + ((Zb-Zc)**2))**0.5
AC = (((Xa-Xc)**2) + ((Ya-Yc)**2) + ((Za-Zc)**2))**0.5
# Direction cosines of AB(ABi,ABj,ABk)
ABi = (Xb-Xa) / AB
ABj = (Yb-Ya) / AB
ABk = (Zb-Za) / AB
# Direction cosines of AC(ACi,ACj,ACk)
ACi = (Xc-Xa) / AC
ACj = (Yc-Ya) / AC
ACk = (Zc-Za) / AC
# Cosine of angle BAC
cosBAC = (AB**2 + AC**2-BC**2) / (2 * AB * AC)
AD = cosBAC * AC
CD = (AC**2-AD**2)**0.5
# Position of point D, which is C projected normally onto AB
Xd = Xa + (AD * ABi)
Yd = Ya + (AD * ABj)
Zd = Za + (AD * ABk)
# Direction cosines of CD(Cdi,CDj,CDk)
CDi = (Xc-Xd) / CD
CDj = (Yc-Yd) / CD
CDk = (Zc-Zd) / CD
# Direction cosines of normal to AB and CD
# to be used for rotations of circle center
Ni = (ABk * CDj)-(ABj * CDk)
Nj = (ABi * CDk)-(ABk * CDi)
Nk = (ABj * CDi)-(ABi * CDj)
# # Diameter of circumscribed circle of a triangle is equal to the
# the length of any side divided by sine of the opposite angle.
# This is done in a coordinate system where X is colinear with AB, Y is // to CD,
# and Z is the normal (N) to X and Y, and the origin is point A
# R = D / 2
sinBAC = (1-cosBAC**2)**0.5
R = (BC / sinBAC) / 2
# center of circumscribed circle is point E
X2e = AB / 2
Y2e = (R**2-X2e**2)**0.5
Z2e = 0
# Transform matrix
# Rotations Translations
#
# ABi , ABj , ABk Xa
# CDi , CDj , CDk Ya
# Ni , Nj , Nk Za
#
# Position of circle center in absolute axis system
X_center = Xa + (X2e * ABi) + (Y2e * CDi) + (Z2e * Ni)
Y_center = Ya + (X2e * ABj) + (Y2e * CDj) + (Z2e * Nj)
Z_center = Za + (X2e * ABk) + (Y2e * CDk) + (Z2e * Nk)
#
center = Vector(X_center,Y_center,Z_center)

# Create the center node
# mw.new_node(center)

X_center_mm = str(X_center*1000)
Y_center_mm = str(Y_center*1000)
Z_center_mm = str(Z_center*1000)
R_mm = str(R*1000)
D_mm = str(R*2*1000)

STR_center = str("Coordinates of center" + "\n\n" + "X = " + X_center_mm + " mm" + "\n" + "Y = " + Y_center_mm + " mm" + "\n" + "Z = " +  Z_center_mm + " mm" + "\n\n" + "Radius =" + R_mm + " mm" + "\n" + "Diameter =" + D_mm + " mm")
mw.message(STR_center)

Victor

Thanks @Sergio. No way to find the current length unit, sorry.

cwharpe

Spider Beams tool. Draws Line2 elements from a center Hub (origin) to a selected cloud of nodes.

Transition easily from line elements to shells. Save your node budget for areas of interest. Hyper-stiff Spider Beams can approximate behavior of an RBE2 element if set up with care. In my tests I increased Modulus of Spiders between 10e3 to 10e6 greater than the part, chose Spider widths that would fill the volume spanned, and came within 3% (form, displacement, vM) comparing actual Nastran RBE2 elements against Mecway Internal solver.

Same behavior is seen with a meshed hyper-stiff diaphram, but I think Spiders are quicker to set up, and maybe a fun little script to run.

cwharpe

Inspired by others to work on scripting skills. Made a couple modifications on (two) clever scripts previously posted by others above:

1) AnchorSkew (Thanks @Victor --Skew tool). Noticed the skewed mesh would wind up scaled off the display the further it originated from (0,0,0). User now pins the mesh down by effectively picking the Anchor Node last. Also revised the Axis selection method and added three more skewing directions.

2) Align_Nodes (Thanks @JohnM --Modify selected nodes to the same constant X value). Included user selection to translate nodes to either X,Y,or Z constant value planes. Also included option to pre-select endpoints of any line in 3D space, then snap all nodes in the selection cloud to the line through perpendicular translation. It's a cool result, but I'm still wondering how I'll use it. (Maybe more useful to snap a cloud of nodes to any plane in 3D space by selecting (3) nodes? Hmmmm..... @Victor -- what about a 3-pt. pick option in Project-to-Surface?)

I'll listen to any scripting suggestions offered. Much to learn.

P.S. Where did all the good spaghetti coding GOTO when I was away?

Sergio

Well well well, @cwharpe, could it be that I herd that you are open to script suggestion??? :-) I have some in mind that could be very usefull, but no code skills at my side. What about:

1) A script to convert all the element named selections (groups) to component? When we mesh in Salome, part meshes are translated and imported as groups of elements and not Mecway components. Maybe in Salome the groups could be renamed to a specific name patron to make easier to find it in the script
2) A script to read the actual liml file, look for the begining of each component definition, and count the elements/nodes and write to an external file (csv or something easiliy imported in Excel), this is for making reports. Something like component_1, xxx (nodes), xxx (elements). Point extra if the material (with properties) can be included

Regards!

cwharpe

@Sergio:

Ha! I was thinking "script suggestions" for better coding technique. I can see how it might sound I was offering to write more code.

Off the cuff, for wish 1) & 2) I'd say the cleanest file read would be from the Mecway generated CCX_ Input_Deck. Looking for keywords like *ELSET,ELSET=(name), *NSET,NSET=(name), *SURFACE,NAME=(name), etc. Read in element numbers for, say, "SADDLE_PLATE" Named Range or something until the next asterisk (*) is encountered, and so on. The old Fortran Rip-and-Strip data mining techniques come to mind. @Victor's API functions now allow for component name creation, plus element, material property assignments to components. From the short-end of the telescope it doesn't look too far away ???

The CCX_Input_Deck_Mesh_Only is an easier read-in if all you want are nodes & elements.

Post a short Salome file that can be read into Mecway. Maybe the forum at large might have a few more insights.

Victor

@cwharpe, did you know that you can project a set of nodes onto a plane parallel to the 3 global planes by using Mesh tools -> Node coordinates? If you leave some components blank, they're not altered.

As for 3-point plane, I think it's almost as easy to create a temporary tri3 element then project onto surface onto that?

@sergio, the API doesn't yet have function to enumerate the named selections or components, so both your suggestions are difficult and probably require reading the file yourself like cwharpe suggests.

cwharpe

@Victor:
Thanks for the Node coordinates tip. Question about Project to Surface in new thread this date.

Victor

Convert element selections to components

# Convert element selections to components. # Caution - can go wrong if there are node selections in the model. # Copyright and related rights waived via CC0. assert mw.version() >= 17, "Wrong Mecway version" for name in mw.named_selections(): contents = mw.named_selection(name) # Check if it possibly contains element IDs. # Can have false positives for node selections. is_elements = True number_of_elements = len(mw.all_elements()) for thingId in contents: if isinstance(thingId, FaceId): is_elements = False elif thingId > number_of_elements: is_elements = False # Create the component if is_elements: component_name = mw.new_component(name) for element_id in contents: mw.set_element_component(element_id, component_name)

Sergio

Hi @Victor, I have modified this last script to create a new material and assing to every component, and it works very well, the problem is that when I try to assing a Young Modulus of 205000000000 Pa (205000 Mpa) to the material (using mw.set_material_property("STEEL", "youngsmodulus", 205000000000)), an error appear:

Error from Microsoft.VisualBasic:
System.InvalidCastException
Conversion from type "BigInteger' to type 'Double' is not valid

With a value of 205000 it works well, but in the GUI this is 205000 Pa and not MPa

Victor

A workaround is to put a decimal point in it like 205000000000.0

Sergio

It works! Now I can import a big UNV assembly mesh into an empty Mecway file, and get all groups of elements converted to components. I have tunned the script to create and assign a basic steel material and assign to each component, and then delete all element selections and empty components. Even all the components have now a fancy name with leading zeros (so ordering is...beautifull)

Victor

Results extraction at a given coordinate without node

# Interpoloates stress XX at the point in global coordinates (X, Y, Z) = (0.5, 0, 0.5)
assert mw.version() >= 18, "Wrong Mecway version"
value = solution.interpolate("stressxx", Vector(0.5,0,0.5))
mw.message(str(value))

If you run it as is with the example model, it uses (0.5, 0, 0.5) which is in the center of the face that has 222 Pa stress.

cwharpe

PURSE a set of edge nodes together with Line elements

Purpose: From a selection of perimeter nodes, script connects line elements (Purse) that can be Auto-meshed, which helps patch regular shaped voids in shell meshes. For high-density meshes surrounding a void I wanted to automate a sometimes intensive mouse-clicking effort in forming the lines. Concept derived from Purse-Seine fishing nets.



Beginning from a Hub (last node in selection) script performs a brute force distance calc. to all neighbors, then draws line to that nearest node as a newly created [PURSE] component. Makes that node the new Hub, drops out the used up nodes, then progressively updates the shrinking list until the curve is finally closed (or left open) by user.

PURSE works admirably in most cases, but has trouble with narrow corners or slender cracks because the nearest distance might bridge across said void instead of follow the edge. To be fair, there IS an existing routine discussed in the Forum/Manual that works for Any void shape, loosely call it The Edge-face Extrude/Merge-Back(onto) method, which results in extruded shells collapsing into edge Line elements. But I've been finding PURSE to be a touch easier to employ for most of the voids I encounter (holes), and I submit it as a noteworthy alternative.

Usage notes:

• Clever modelers might have saved the void edge faces as a Named Range before refining model, which auto updates all the voids edge nodes.
  
• Hide-all-but-this in Component tree will highlight just the [PURSE] beams for inspection.
  
• UNDO is your friend.
  
• CTRL Re-Pick the last node in your selection makes it the starting Hub node and possibly a different curve shape for irregularly spaced node selections.
  

Pursed sides then Auto-meshed.

disla

Wow!!. It worked nice and soft. Thanks for sharing.

cwharpe

@disla -- you consistently demonstrate some very graphics. Glad it worked.