-
Boman Romain authoredBoman Romain authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
bonemodel.py 18.87 KiB
# -*- coding: utf-8 -*-
# Bone FEA model (formely: "mandible model")
#
# main Metafor script with default parameters.
# R.Boman
from wrap import *
import boneload
import os
def parms(d={}):
"""default parameters of the model
"""
p = {}
# folder containing meshes
path = 'dolicorhynchops/10k'
# mandible mesh (.ply/.geo/.msh)
p['bone'] = f'{path}/mandible.ply'
p['muscles'] = [
{
'file': f'{path}/Lmuscle.ply',
'force': 100.,
'focalpt': f'{path}/LmuscleF.off',
'method': 'T' # loading model: 'U', 'T', 'T+N', 'D' or 'N'
},
{
'file': f'{path}/Rmuscle.off',
'force': 100.,
'focalpt': f'{path}/RmuscleF.off',
'method':'T'
}
]
# list of fixed nodes
# with "name" a description string used for output purposes,
# "nodes" a list of nodes (coordinates) or a mesh file,
# "direction" a list of constrained degrees of freedom
p['fixations'] = [
{
'name': 'contact_pts',
'nodes': [[-10.20362, -17.46838, -229.9061],
[-11.92466, 26.3042, -229.5354]],
'direction': 'x'
},
{
'name': 'axis_pt1',
'nodes': f'{path}/LTMJ.off',
'direction': ['x', 'y', 'z']
},
{
'name': 'axis_pt2',
'nodes': [-8.716309, 79.13171, 233.8385],
'direction': ['x', 'z']
}
]
p['loads'] = [] # additional loads
p['extractors'] = [] # additional extractors
p['food'] = 'fixteeth' # type of food: 'fixteeth' or 'cylinder'
# rigid cylinder (used if p['food']=='cylinder')
R = 20
p['cylinder'] = {
'origin': [-R-8, 0, -150],
'radius': R,
'width': 100,
'friction': 0.1,
'penalty': 1e3
}
# material properties
p['density'] = 1.850e-9 # [T/mm³] - bone: 1.850 kg/l
p['Young'] = 17000. # [MPa] elastic modulus - bone: 17-20 GPa
p['Poisson'] = 0.3 # [-] Poisson's ratio
# numerical parameters
p['tolNR'] = 1e-6 # [-] equilibrium tolerance
p['dt0'] = 1.0 # [s] time step size
p['narch'] = 1 # number of archives
# gmsh toolbox
p['use_gmshOld'] = False # use old gmsh interface
p.update(d)
return p
def getMetafor(p={}):
p = parms(p)
import json
print(f'parameters={json.dumps(p, indent=4, sort_keys=True)}')
metafor = Metafor()
domain = metafor.getDomain()
geometry = domain.getGeometry()
geometry.setDim3D()
# load main mandible mesh (volume or surface - it will be meshed in 3D)
import_mesh(domain, p['bone'], p['use_gmshOld'])
# -- define groups
# group #1 is the group containing all volume elements
groupset = geometry.getGroupSet()
groups = {}
groups['all'] = groupset(1)
print("the mandible mesh has %d nodes" % groups['all'].getNumberOfMeshPoints())
# extract external surface of the mesh => groups['surf']
# nods_no: tags of the nodes
# nods_pos: coordinates of these nodes
groups['surf'] = groupset.add(Group(groupset.getLargestNo()+1))
groups['surf'].addMeshPointsFromObject(groups['all'], BoundarySelector())
nods_no, nods_pos = extract_nodes_from_group(groups['surf'])
# load muscle groups
mgroups = {} # stores muscle group data and loads
for muscle in p['muscles']:
# load focal point / loading direction
if isinstance(muscle['focalpt'], str):
fullpath = os.path.join(os.path.dirname(__file__), muscle['focalpt'])
focalnodes, _ = boneload.load_msh(fullpath)
else: # coordinates in array or tuple
focalnodes = [ muscle['focalpt'] ]
# load surface mesh => groups[name (from filename)]
name, nodes, tris, ntags = \
create_group(muscle['file'], nods_no, nods_pos, domain, groups)
# create loads on the surface
loads = boneload.create_loads(nodes, tris, \
total_force=muscle['force'], target=focalnodes[0], method=muscle['method'])
# store data in a structure
mgroups[name] = MuscleGroup(nodes, tris, ntags, loads)
# create groups of fixed nodes
for fix in p['fixations']:
create_group(fix['nodes'], nods_no, nods_pos, domain, groups, fix['name'])
# create groups of additional extractors
for extr in p['extractors']:
create_group(extr['nodes'], nods_no, nods_pos, domain, groups, extr['name'])
# --------------------------------------------------------------------------
# material
materials = domain.getMaterialSet()
materials.define (1, ElastHypoMaterial)
materials(1).put(MASS_DENSITY, p['density'])
materials(1).put(ELASTIC_MODULUS, p['Young'])
materials(1).put(POISSON_RATIO, p['Poisson'])
prp1 = ElementProperties(TetraVolume3DElement)
prp1.put(MATERIAL, 1)
prp1.put (CAUCHYMECHVOLINTMETH, VES_CMVIM_STD)
app = FieldApplicator(1)
app.push(groups['all'])
app.addProperty(prp1)
domain.getInteractionSet().add(app)
# boundary conditions
fct = PieceWiseLinearFunction()
fct.setData(0.0, 0.0)
fct.setData(1.0, 1.0)
# add Dirichlet conditions (fixations)
txt2key = { 'x': TX, 'y': TY, 'z': TZ }
for fix in p['fixations']:
for d in fix['direction']:
domain.getLoadingSet().define(groups[fix['name']], Field1D(txt2key[d],RE), 0.0)
# add Neumann conditions (loads) from muscle loads
mshpts = geometry.getMesh().getPointSet()
for name, gr in mgroups.items():
# print(f'len(gr.ntags)={len(gr.ntags)}')
for i,load in enumerate(gr.loads):
# print(f'i={i}')
target = mshpts(gr.ntags[i])
domain.getLoadingSet().define(target, Field1D(TX,GF1), load.x[0], fct)
domain.getLoadingSet().define(target, Field1D(TY,GF1), load.x[1], fct)
domain.getLoadingSet().define(target, Field1D(TZ,GF1), load.x[2], fct)
# add Neumann conditions (loads) from external loads
dirs = [TX,TY,TZ]
for load in p['loads']:
name, nodes, tris, ntags = \
create_group(load['nodes'], nods_no, nods_pos, domain, groups, load['name'])
for ntag in ntags:
for i in range(3):
target = mshpts(ntag)
domain.getLoadingSet().define(target, Field1D(dirs[i],GF1), load['values'][i], fct)
# type of food
if p['food']=='cylinder':
# contact with a rigid cylinder
cylinder = Cylinder(domain, center=p['cylinder']['origin'], \
width=p['cylinder']['width'], \
radius=p['cylinder']['radius'])
mater2 = materials.define(2, CoulombContactMaterial)
mater2.put(COEF_FROT_DYN, p['cylinder']['friction'])
mater2.put(COEF_FROT_STA, p['cylinder']['friction'])
mater2.put(PEN_TANGENT, p['cylinder']['penalty'])
mater2.put(PEN_NORMALE, p['cylinder']['penalty'])
mater2.put(PROF_CONT, p['cylinder']['radius']/2.)
prp2 = ElementProperties(Contact3DElement)
prp2.put(MATERIAL, 2)
ci = RdContactInteraction(2)
ci.push(groups['surf'])
ci.setTool(cylinder.skin)
ci.addProperty(prp2)
domain.getInteractionSet().add(ci)
elif p['food']=='fixteeth':
pass
else:
# p['food'] is a filename
ply2mtf(domain, p['food'])
# Time integration scheme
ti = AlphaGeneralizedTimeIntegration(metafor)
metafor.setTimeIntegration(ti)
tsm = metafor.getTimeStepManager()
tsm.setInitialTime(0.0, p['dt0'])
tsm.setNextTime(1.0, p['narch'], 1.0)
mim = metafor.getMechanicalIterationManager()
mim.setResidualTolerance(p['tolNR'])
# parallel solver
solman = metafor.getSolverManager()
solver = DSSolver()
solman.setSolver(solver)
# enable full parallelism
StrVectorBase.useTBB()
StrMatrixBase.useTBB()
ContactInteraction.useTBB()
# history curves
# by default, extract all forces on all groups
valuesmanager = metafor.getValuesManager()
idx=1
valuesmanager.add(idx, MiscValueExtractor(metafor,EXT_T),'time'); idx+=1
for name, group in groups.items():
# skip the main mesh and its surface
if name in ['all','surf']: continue
# skip all additional extractors
if name in [e['name'] for e in p['extractors']]: continue
# muscles and loads should remain
for key, c in zip([TX,TY,TZ],['x','y','z']):
fname = f'{name}_F{c}'
valuesmanager.add(idx, DbNodalValueExtractor(group, Field1D(key,GF1)), SumOperator(), fname)
metafor.getTestSuiteChecker().checkExtractor(idx)
idx+=1
if p['food'] == 'cylinder':
for key, c in zip([TX,TY,TZ],['x','y','z']):
fname = f'cylinder_F{c}'
valuesmanager.add(idx, InteractionValueExtractor(ci, key), fname)
metafor.getTestSuiteChecker().checkExtractor(idx)
idx+=1
# add additional extractors
# conversion dicts
var2field_db = {
'dx': Field1D(TX, RE),
'dy': Field1D(TY, RE),
'dz': Field1D(TZ, RE),
'x0': Field1D(TX, AB),
'y0': Field1D(TY, AB),
'z0': Field1D(TZ, AB),
'fx': Field1D(TX, GF1),
'fy': Field1D(TY, GF1),
'fz': Field1D(TZ, GF1)
}
var2field_if = {
'sigvm': IF_EVMS,
'sigxx': IF_SIG_XX,
'sigyy': IF_SIG_YY,
'sigzz': IF_SIG_ZZ,
'sigxy': IF_SIG_XY,
'sigyz': IF_SIG_YZ,
'sigxz': IF_SIG_XZ,
'epsxx': IF_GL_STRAIN_XX,
'epsyy': IF_GL_STRAIN_YY,
'epszz': IF_GL_STRAIN_ZZ,
'epsxy': IF_GL_STRAIN_XY,
'epsyz': IF_GL_STRAIN_YZ,
'epsxz': IF_GL_STRAIN_XZ
}
for extr in p['extractors']:
# if only one variable is given, convert it to a list
if isinstance(extr['variables'], str):
extr['variables'] = [extr['variables']]
for var in extr['variables']:
fname = f'{extr["name"]}_{var}'
target = groups[extr['name']]
# nodal values
if var.lower() in var2field_db:
field = var2field_db[var.lower()]
valuesmanager.add(idx, DbNodalValueExtractor(target, field), fname)
elif var.lower() in var2field_if:
field = var2field_if[var.lower()]
valuesmanager.add(idx, IFNodalValueExtractor(target, field), fname)
else:
raise Exception(f'Unknown variable: {var}')
# if the target is a single node, display the extractor [TSC-EXT]
# otherwise, the user should read the ascii file in the workspace.
if target.getNumberOfMeshPoints() == 1:
metafor.getTestSuiteChecker().checkExtractor(idx)
idx+=1
return metafor
# ------------------------------------------------------------------------------
# utility functions involving Metafor objects
def import_mesh(domain, filename, use_gmshOld=True):
"""import the volume mesh
input file type can be:
- .msh: the mesh is simply loaded
- .geo: the mesh is generated
- .ply/.stl: a .geo is written in the workspace and loaded into gmsh
"""
# import gmsh file (accept .geo, .msh, .ply or .stl)
mandible_gmsh = filename
# make the path absolute
if not os.path.isabs(mandible_gmsh):
mandible_gmsh = os.path.join(os.path.dirname(__file__), mandible_gmsh)
# check whether file exists
if not os.path.isfile(mandible_gmsh):
raise Exception(f'{mandible_gmsh} not found!')
# extract extension
noext,ext = os.path.splitext(mandible_gmsh)
if ext in ['.ply', '.stl']:
plyfile = mandible_gmsh
# create 2 workspace files
basename = os.path.join(os.getcwd(), os.path.basename(noext))
mshfile = f"{basename}.msh"
geofile = f"{basename}.geo"
# if msh exists, we can load it directly if the ply file is older than
# the msh file. Otherwise, we create a geo file
if os.path.isfile(mshfile) and \
os.path.getmtime(mshfile) > os.path.getmtime(plyfile):
mandible_gmsh=mshfile
else:
print(f'creating {geofile}')
f = open(geofile, "w")
f.write(f'Merge "{plyfile}";\n')
f.write('Surface Loop(1) = {1};\n')
f.write('Volume(1) = {1};\n')
f.write('Physical Volume(1) = { 1 };\n')
f.write('Mesh.MshFileVersion = 2.2;\n')
f.close()
mandible_gmsh=geofile
elif ext in ['.geo', '.msh']:
pass
else:
raise Exception(f'Unknown extension: {ext}, please use .ply, .stl, .msh or .geo')
print(f'importing {mandible_gmsh}')
f = os.path.join(os.path.dirname(__file__), mandible_gmsh)
if use_gmshOld:
from toolbox.gmshOld import GmshImport
mesher = GmshImport(f, domain)
else:
from toolbox.gmsh import GmshImport # FIXME: does not work!
mesher = GmshImport(f, domain)
# mesher.verb = True # debug
mesher.writeMshFile = True
mesher.execute()
def extract_nodes_from_group(group):
"""input: a Group object of Metafor (a group of mesh points)
output: 2 arrays containing the node numbers and the node coordinates
"""
nods_no = []
nods_pos = []
for i in range(group.getNumberOfMeshPoints()):
nods_no.append(group.getMeshPoint(i).getNo())
pos = group.getMeshPoint(i).getPos(Configuration().currentConf)
nods_pos.append([ pos.get1(), pos.get2(), pos.get3()])
# print(f'nods_no={nods_no}')
return nods_no, nods_pos
def create_group(mesh_file_or_coord, all_no, all_coords, domain, groups, gname=None):
"""loads the mesh file "mesh_file_or_coord" and identify vertices among vertices
stored in all_coords.
Then, the routine creates a group named "gname" (built from mesh_file_or_coord if None)
mesh_file_or_coord can be either:
a filename (.off, .ply, .stl)
an array of points: [ [x1,y1,z1], [x2,y2,z2], ... ]
one single point: [x1,y1,z1]
"""
if isinstance(mesh_file_or_coord, str):
# load mesh file
print(f'\n* create_group {mesh_file_or_coord}...')
fullpath = os.path.join(os.path.dirname(__file__), mesh_file_or_coord)
nodes, tris = boneload.load_msh(fullpath)
elif type(mesh_file_or_coord) in [list, tuple]:
print(f'\n* create_group {gname}...')
if type(mesh_file_or_coord[0]) in [list, tuple]:
# array of points
nodes = mesh_file_or_coord
else:
# one single point
nodes = [ mesh_file_or_coord ]
print(f'\tnodes={nodes}')
print(f'\tmesh_file_or_coord={mesh_file_or_coord}')
tris = []
else:
raise Exception(f'Data is not understood: {mesh_file_or_coord}')
# identify nodes of the mesh among the given vertices of the main mesh
# ntags = boneload.identify_nodes_brutal(nodes, all_no, all_coords)
ntags = boneload.identify_nodes_SaP(nodes, all_no, all_coords)
# create a metafor group
groupset = domain.getGeometry().getGroupSet()
if gname is None:
gname, _ = os.path.splitext(os.path.basename(mesh_file_or_coord))
group = groupset.add(Group(groupset.getLargestNo()+1))
if gname in groups:
raise Exception(f'create_group: group "{gname}" already exists!')
groups[gname] = group
# add all the identified vertices to the group
nodeset = domain.getGeometry().getMesh().getPointSet()
for tag in ntags:
group.addMeshPoint(nodeset(tag))
print(f'\tgroup #{group.getNo()} ({gname}) ' \
f'created from {mesh_file_or_coord} ({group.getNumberOfMeshPoints()} nodes)')
return gname, nodes, tris, ntags
class MuscleGroup:
"""convenient structure for storing muscle data and loads
"""
def __init__(self, nodes, tris, ntags, loads):
self.nodes = nodes
self.tris = tris
self.ntags = ntags
self.loads = loads
def ply2mtf(domain, meshfile):
"""load a surface into Metafor's memory from a mesh file (ply, off or stl)
"""
fullpath = os.path.join(os.path.dirname(__file__), meshfile)
nodes, tris = boneload.load_msh(fullpath)
mesh = domain.getGeometry().getMesh()
meshPtSet = mesh.getPointSet()
# create group
groupset = domain.getGeometry().getGroupSet()
grp = groupset.add(Group(groupset.getLargestNo()+1))
# create mesh points
idx0 = meshPtSet.getLargestNo()+1
for i, n in enumerate(nodes):
pt = meshPtSet.define(idx0+i, n[0], n[1], n[2])
grp.addMeshPoint(pt)
skinset = mesh.getSkinSet()
skin = skinset.add(Skin(skinset.getLargestNo()+1))
print(f'ply2mtf: creating Skin #{skin.getNo()}')
sideset = mesh.getSideSet()
# create triangles
idx1 = sideset.getLargestNo()+1
for i,tnods in enumerate(tris):
mesh.define(idx1+i, CELL_TRIANGLE, grp, [ meshPtSet(idx0+n) for n in tnods])
skin.push(sideset(idx1+i))
print(f'ply2mtf: creating Group #{grp.getNo()} meshed with {grp.getNumberOfMeshPoints()} nodes')
class Cylinder:
def __init__(self, domain, center, width, radius):
geometry = domain.getGeometry()
ox, oy, oz = center
pointset = geometry.getPointSet()
idx = pointset.getLargestNo()
p1 = pointset.define(idx+1, ox, oy-width/2, oz-radius)
p2 = pointset.define(idx+2, ox+radius, oy-width/2, oz)
p3 = pointset.define(idx+3, ox, oy-width/2, oz+radius)
p4 = pointset.define(idx+4, ox-radius, oy-width/2, oz)
p11 = pointset.define(idx+11, ox, oy+width/2, oz-radius)
p12 = pointset.define(idx+12, ox+radius, oy+width/2, oz)
p13 = pointset.define(idx+13, ox, oy+width/2, oz+radius)
p14 = pointset.define(idx+14, ox-radius, oy+width/2, oz)
curveset = geometry.getCurveSet()
idx = curveset.getLargestNo()
l1 = curveset.add(Arc(idx+1, p1, p2, p3))
l2 = curveset.add(Arc(idx+2, p3, p4, p1))
l11 = curveset.add(Arc(idx+11, p11, p12, p13))
l12 = curveset.add(Arc(idx+12, p13, p14, p11))
l21 = curveset.add(Line(idx+21, p1, p11))
l22 = curveset.add(Line(idx+22, p3, p13))
wireset = geometry.getWireSet()
idx = wireset.getLargestNo()
w1 = wireset.add(Wire(idx+1, [l1, l22, l11, l21]))
w2 = wireset.add(Wire(idx+2, [l2, l21, l12, l22]))
srfset = geometry.getSurfaceSet()
idx = srfset.getLargestNo()
srf1 = srfset.add(Ruled(idx+1, l1, l11))
srf2 = srfset.add(Ruled(idx+2, l12, l2))
sideset = geometry.getSideSet()
idx = sideset.getLargestNo()
s1 = sideset.add(Side(idx+1, [w1])); s1.setSurface(srf1)
s2 = sideset.add(Side(idx+2, [w2])); s2.setSurface(srf2)
skinset = geometry.getSkinSet()
idx = skinset.getLargestNo()
sk1 = skinset.add(Skin(idx+1, [s1, s2]))
self.skin = sk1