Preprocessing with AxisVM

We are going to read a tetrahedral mesh from a .vtk file, send its surface into AxisVM to select points of applications for the applications of natural and essential boundary conditions, identify the selected items, build a finite element model, do a simple analysis and save the results in a simple, but powerful format.

The mesh for this tutorial is among the examples of AxisVM:

from axisvm import examples
import pyvista as pv
import numpy as np

from polymesh import PolyData, CartesianFrame, PointData
from polymesh.cells import TET4 as CellData
from polymesh.tri.triutils import edges_tri
from polymesh.topo import unique_topo_data, detach

vtkpath = examples.download_stand_vtk()
pd = PolyData.read(vtkpath)
c = pd.coords() / 1000
t = pd.topology()[:, :4]
c, t = detach(c, t)
csys = CartesianFrame(dim=3)

# pointdata and celldata
pd = PointData(coords=c, frame=csys)
cd = CellData(topo=t, frames=csys)

mesh = PolyData(pd, cd, frame=csys)
surface = mesh.surface()
surface.nummrg()

PolyData({})

surface.plot(notebook=False)

We can access the local and global indices of the surface mesh. This migh be handy if we want to convert data between the two sets.

len(surface.pd), surface.pd.id, surface.pd.gid

(2869,
 array([   0,    1,    2, ..., 2866, 2867, 2868]),
 <Array [0, 1, 2, 3, ... 2885, 2886, 2887, 2888] type='2869 * int32'>)

from polymesh.topo import rewire, detach

coords, topo = surface.coords(), surface.topology()
edges, edgeIDs = unique_topo_data(edges_tri(topo))

# in AxisVM, numbering starts with 1
edges += 1
edgeIDs += 1

from axisvm.com.client import start_AxisVM
axvm = start_AxisVM(visible=True, daemon=True)

For the example at hand, we build the AxisVM model with an arbitrary material model.

modelId = axvm.Models.New()
axm = axvm.Models.Item[modelId]
axm.Settings.EditingTolerance = -1
wdir = ""

from axisvm.com.tlb import RPoint3d
foo = lambda x : RPoint3d(x=x[0], y=x[1], z=x[2])
axm.BeginUpdate()
axm.Nodes.BulkAdd(list(map(foo, coords)))
axm.EndUpdate()

from axisvm.com.tlb import lgtStraightLine, RLineData

def gen_line(edge):
    return RLineData(
        NodeId1 = edge[0],
        NodeId2 = edge[1],
        GeomType = lgtStraightLine
    )

axm.BeginUpdate()
axm.Lines.BulkAdd(list(map(gen_line, edges)))
axm.EndUpdate()

from axisvm.com.tlb import vTop
axm.View = vTop
axm.FitInView()

from axisvm.com.tlb import ndcEuroCode
axm.Settings.NationalDesignCode = ndcEuroCode
matId = axm.Materials.AddFromCatalog(ndcEuroCode, "S 235")

from axisvm.com.tlb import RSurfaceAttr, lnlTensionAndCompression, \
    RResistancesXYZ, schLinear, stShell, RElasticFoundationXYZ, \
    RNonLinearityXYZ, RSurface

SurfaceAttr = RSurfaceAttr(
    Thickness=1,
    SurfaceType=stShell,
    RefZId=0,
    RefXId=0,
    MaterialId=matId,
    ElasticFoundation=RElasticFoundationXYZ(0, 0, 0),
    NonLinearity=RNonLinearityXYZ(lnlTensionAndCompression,
                                  lnlTensionAndCompression,
                                  lnlTensionAndCompression),
    Resistance=RResistancesXYZ(0, 0, 0),
    Charactersitics=schLinear)

def gen_surface(edges):
    return RSurface(
        N=3,
        LineIndex1 = edges[0],
        LineIndex2 = edges[1],
        LineIndex3 = edges[2],
        Attr = SurfaceAttr,
        DomainIndex = 0
    )
axm.BeginUpdate()
axm.Surfaces.BulkAdd(list(map(gen_surface, edgeIDs)))
axm.EndUpdate()

0

We can add data to the AxisVM Model. The following cell adds the global ids of the nodes to the model. The attached data is then stored with the model file in binary form.

attachments = axm.Nodes.Attachments
gid = surface.pd.gid
nN = axm.Nodes.Count
for i in range(nN):
    attachments.AddData('gid', i+1, [gid[i]])

Use AxisVM to select points of applications for nodal loads and supports. You will be prompted to select node in AxisVM using the GUI. Each selection goes into a custom part folder.

from axisvm.com.tlb import RPartItem, pitNode
axvm.BringToFront()
CustomParts = axm.CustomParts
CustomPartsFolder = CustomParts.RootFolder
i = CustomPartsFolder.AddSubFolder('SigmaEpsilon')
Folder = CustomPartsFolder.SubFolder[i]
iebc = Folder.AddPart('EBC', [])[-1]
inbc = Folder.AddPart('NBC', [])[-1]

nodes_f = axm.Nodes.select_IDs(msg='Select nodes where nodal loads are to be imposed!')
parts = list(map(lambda id : RPartItem(ItemType=pitNode, Id=id), nodes_f))
Folder.AddPartItemsToPart(inbc, parts)
axm.SelectAll(False)
axm.Refresh()

0

nodes_u = axm.Nodes.select_IDs(msg='Select nodes where displacement penalties are to be enforced!')
parts = list(map(lambda id : RPartItem(ItemType=pitNode, Id=id), nodes_u))
Folder.AddPartItemsToPart(iebc, parts)
axm.SelectAll(False)
axm.Refresh()

0

Read the selected nodes from their custom part folders.

nbc_items = Folder.GetPart(inbc)[0]
nbc_ids = list(map(lambda i : i.Id, filter(lambda i : i.ItemType == 0, nbc_items)))
ebc_items = Folder.GetPart(iebc)[0]
ebc_ids = list(map(lambda i : i.Id, filter(lambda i : i.ItemType == 0, ebc_items)))

Build the finite element model and perform a linear elastic analysis assuming small strains and displacements.

from sigmaepsilon.solid import Structure, PointData, FemMesh
from polymesh.space import StandardFrame
from sigmaepsilon.solid.fem.cells import TET4 as CellData
from neumann.array import repeat

# the global frame of the workspace
GlobalFrame = StandardFrame(dim=3)

# coordinates and topology
coords = mesh.coords()
topo = mesh.topology()

# essential boundary conditions
ebc_gids = [gid[id-1] for id in ebc_ids]
fixity = np.zeros((coords.shape[0], 6), dtype=bool)
fixity[ebc_gids, :3] = True
fixity[:, 3:] = True

# natural boundary conditions
F = 10
nbc_gids = [gid[id-1] for id in nbc_ids]
loads = np.zeros((coords.shape[0], 6))
loads[nbc_gids, 2] = -F

# pointdata
pd = PointData(coords=coords, frame=GlobalFrame,
               loads=loads, fixity=fixity)

# celldata
frames = repeat(GlobalFrame.show(), topo.shape[0])
cd = CellData(topo=topo, frames=frames)

# define a stiffness matrix for Hooke's model
E = 12000.
nu = 0.2
Hooke = np.array([
    [1, nu, nu, 0, 0, 0],
    [nu, 1, nu, 0, 0, 0],
    [nu, nu, 1, 0, 0, 0],
    [0., 0, 0, (1-nu)/2, 0, 0],
    [0., 0, 0, 0, (1-nu)/2, 0],
    [0., 0, 0, 0, 0, (1-nu)/2]]) * (E / (1-nu**2))

mesh = FemMesh(pd, cd, model=Hooke, frame=GlobalFrame)

structure = Structure(mesh=mesh)

structure.linear_static_analysis()
dofsol = structure.nodal_dof_solution()

Plot the model, coloured with vertical displacements:

structure.mesh.plot(notebook=True, scalars=dofsol[:, 2], backend='k3d')

f:\GitHub\sigmaepsilon\.sigeps\lib\site-packages\traittypes\traittypes.py:97: UserWarning: Given trait value dtype "float32" does not match required type "float32". A coerced copy has been created.
  warnings.warn(