gmsh + pymapdl + cadquery

lego.py

#!/usr/bin/env python
# coding: utf-8

# In[12]:


import cadquery as cq
import pyvista as pv

from ansys.mapdl.reader  import save_as_archive
from ansys.mapdl.core import launch_mapdl


# In[13]:


#####
# Inputs
######
lbumps = 6     # number of bumps long
wbumps = 2     # number of bumps wide
thin = True    # True for thin, False for thick

#
# Lego Brick Constants-- these make a Lego brick a Lego :)
#
pitch = 8.0
clearance = 0.1
bumpDiam = 4.8
bumpHeight = 1.8
if thin:
    height = 3.2
else:
    height = 9.6

t = (pitch - (2 * clearance) - bumpDiam) / 2.0
postDiam = pitch - t  # works out to 6.5
total_length = lbumps*pitch - 2.0*clearance
total_width = wbumps*pitch - 2.0*clearance

# make the base
s = cq.Workplane("XY").box(total_length, total_width, height)

# shell inwards not outwards
s = s.faces("<Z").shell(-1.0 * t)

# make the bumps on the top
s = (s.faces(">Z").workplane().
    rarray(pitch, pitch, lbumps, wbumps, True).circle(bumpDiam / 2.0)
    .extrude(bumpHeight))

# add posts on the bottom. posts are different diameter depending on geometry
# solid studs for 1 bump, tubes for multiple, none for 1x1
tmp = s.faces("<Z").workplane(invert=True)

if lbumps > 1 and wbumps > 1:
    tmp = (tmp.rarray(pitch, pitch, lbumps - 1, wbumps - 1, center=True).
        circle(postDiam / 2.0).circle(bumpDiam / 2.0).extrude(height - t))
elif lbumps > 1:
    tmp = (tmp.rarray(pitch, pitch, lbumps - 1, 1, center=True).
        circle(t).extrude(height - t))
elif wbumps > 1:
    tmp = (tmp.rarray(pitch, pitch, 1, wbumps - 1, center=True).
        circle(t).extrude(height - t))
else:
    tmp = s


# export cad design to .stl file
filename = 'lego_example.step'
cq.exporters.export(tmp, filename)

###############################################################################
# plot just the tesselated CAD

import pyvista as pv
cq.exporters.export(tmp, './lego_example.stl')
mesh = pv.read('./lego_example.stl')


###############################################################################
import gmsh

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 1)

gmsh.model.add("t20")

# Load a STEP file (using `importShapes' instead of `merge' allows to directly
# retrieve the tags of the highest dimensional imported entities):
v = gmsh.model.occ.importShapes(filename)


# Get the bounding box of the volume:
gmsh.model.occ.synchronize()

# Specify a global mesh size and mesh the partitioned model:
sz = 0.2
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", sz/2)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", sz*2)
# gmsh.option.setNumber('Mesh.SecondOrderIncomplete', 1)  # <-- that's it!

gmsh.model.mesh.generate(3)
# gmsh.model.mesh.setOrder(2)  # quadratic
gmsh_out = "from_gmsh.msh"
gmsh.write(gmsh_out)
gmsh.finalize()

###############################################################################
# output to pyvista
grid = pv.read(gmsh_out)

# grid.plot(color='w', show_edges=True, pbr=True, background='w', metallic=1, split_sharp_edges=True, roughness=0.6)
# grid.extract_cells(range(1, 100000, 10)).plot(color='w', show_edges=True)


# use pyvista to read .stl file and save to .cdb file
# reader = pv.get_reader('./lego_example.stl')
# mesh = reader.read()
grid.points
grid.points /= 1000
cdb_filename = "lego_example.cdb"
save_as_archive(cdb_filename, grid, include_surface_elements=False)
grid.plot()


# In[16]:


mapdl = launch_mapdl()

# load a lego file and mesh it
mapdl.cdread("db", cdb_filename)

mapdl.prep7()
print(mapdl.shpp("SUMM"))

# specify shell thickness
# mapdl.sectype(1, "shell")
# mapdl.secdata(0.01)
# mapdl.emodif("ALL", "SECNUM", 1)

# specify material properties
# using aprox values for AISI 5000 Series Steel
# http://www.matweb.com/search/datasheet.aspx?matguid=89d4b891eece40fbbe6b71f028b64e9e
mapdl.units("SI")  # not necessary, but helpful for book keeping
mapdl.mp("EX", 1, 200e9)  # Elastic moduli in Pa (kg/(m*s**2))
mapdl.mp("DENS", 1, 7800)  # Density in kg/m3
mapdl.mp("NUXY", 1, 0.3)  # Poissons Ratio
mapdl.emodif("ALL", "MAT", 1)

# Run an unconstrained modal analysis
# for the first 20 modes above 1 Hz
mapdl.modal_analysis(nmode=20, freqb=1)


# result = mapdl.result
# result.animate_nodal_displacement(0)

# In[ ]:

# 
# mapdl._result_file
# /tmp/ansys_pbofrfnayt/file.rst

from ansys.dpf import core as dpf

model = dpf.Model(mapdl._result_file)
print(model)

results = model.results
displacements = results.displacement()
fields = displacements.outputs.fields_container()

# Finally, extract the data of the displacement field:
field = fields[0]
field.plot()

# disp = field.data
# disp

###############################################################################

mapdl.exit()


# In[ ]: