-
Boman Romain authoredBoman Romain authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
bonemodel2.py 13.60 KiB
# -*- coding: utf-8 -*-
# Bone FEA model (formely: "mandible model")
#
# main cxxfem script with default parameters.
# R.Boman
import cxxfem as fem
from . import boneload
import os
import gmsh
if not gmsh.isInitialized():
gmsh.initialize()
import numpy as np
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
# material properties
p['Young'] = 17000. # [MPa] elastic modulus - bone: 17-20 GPa
p['Poisson'] = 0.3 # [-] Poisson's ratio
p.update(d)
return p
def solve(p={}):
import time
start_time = time.perf_counter()
p = parms(p)
import json
print(f'parameters={json.dumps(p, indent=4, sort_keys=True)}')
pbl = fem.Problem()
# load main mandible mesh (volume or surface - it will be meshed in 3D)
import_mesh(p['bone'])
# -- define groups
# group (3,1)='all' is the group containing all volume elements
groups = {}
groups['all'] = (3, 1)
nodeTags, coord = gmsh.model.mesh.getNodesForPhysicalGroup(*groups['all'])
print(f"the mandible mesh has {len(nodeTags)} nodes")
# extract external surface of the mesh => groups['surf']
# nods_no: tags of the nodes
# nods_pos: coordinates of these nodes
groups['surf'] = (2, 1)
nods_no, nods_pos = gmsh.model.mesh.getNodesForPhysicalGroup(
*groups['surf'])
print(f"the mandible surface has {len(nods_no)} nodes")
nods_pos = np.reshape(nods_pos, (len(nods_no), 3))
# 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, 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)
# save loads to file
with open(f'loads_{name}.tsv','w') as f:
f.write('node#\tfx\tfy\tfz\n')
for i,load in enumerate(loads):
f.write(f"{int(ntags[i])}\t{load.x[0]}\t{load.x[1]}\t{load.x[2]}\n")
# create groups of fixed nodes
for fix in p['fixations']:
create_group(fix['nodes'], nods_no, nods_pos, groups, fix['name'])
# create groups of additional extractors
for extr in p['extractors']:
create_group(extr['nodes'], nods_no, nods_pos, groups, extr['name'])
# --------------------------------------------------------------------------
# print('== ADDING MEDIUM...')
print('')
do_not_delete = []
# material
do_not_delete.append( fem.Medium(pbl, "all", E=p['Young'], nu=p['Poisson']) )
# add Dirichlet conditions (fixations)
# print('== ADDING DIRICHLETs...')
for fix in p['fixations']:
for d in fix['direction']:
do_not_delete.append(fem.Dirichlet(pbl, fix['name'], dof=d.upper(), val=0.0) )
# add Neumann conditions (loads) from muscle loads
# print('== ADDING NEUMANNs...')
for name, gr in mgroups.items():
# print(f'len(gr.ntags)={len(gr.ntags)}')
for i,load in enumerate(gr.loads):
do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'X', load.x[0]) ) # TODO: regler pbl de cast "gr.ntags[i] => size_t"
do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'Y', load.x[1]) ) # TODO: regler "dof=" qui marche pas
do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'Z', load.x[2]) )
# add Neumann conditions (loads) from external loads
for load in p['loads']:
name, nodes, tris, ntags = \
create_group(load['nodes'], nods_no, nods_pos, groups, load['name'])
for ntag in ntags:
for i in range(3):
do_not_delete.append(fem.Neumann(pbl, int(ntag), 'XYZ'[i], load['values'][i]) )
# python only!
# do_not_delete.append(fem.Dirichlet(pbl, "all", dof='T', val=0.0) )
# Time integration scheme
# print (pbl)
# print('SOLVE.......')
pbl.quadrature = 'Gauss1' # 1 GP/tetra is enough
solver = fem.Solver(pbl)
solver.linear_solver = "pardiso"
solver.solve()
post = fem.Post(solver)
post.build_views()
post.show_views([ 'smooth_stress_tensor' ])
post.set_camera_3D()
post.write()
post.deform()
print('extracting results...')
print(f'\tinternal energy = {solver.int_energy:.3f} N.mm')
# fixations (reaction forces)
for fix in p['fixations']:
grpname = fix['name']
v = np.array(post.probe('force_vector', grpname))
v = np.reshape(v, (v.size//3, 3))
v = np.sum(v, axis=0)
print(f'\t{grpname}_Fx = {v[0]:.3f} N')
print(f'\t{grpname}_Fy = {v[1]:.3f} N')
print(f'\t{grpname}_Fz = {v[2]:.3f} N')
# muscle groups
for name, group in mgroups.items():
v=np.array( [0.0, 0.0, 0.0] )
for tag in group.ntags:
v += np.array(post.probe('force_vector', int(tag)))
print(f'\t{name}_Fx = {v[0]:.3f} N')
print(f'\t{name}_Fy = {v[1]:.3f} N')
print(f'\t{name}_Fz = {v[2]:.3f} N')
# additional extractors
var2field = {
'dx': ('X', 'mm'),
'dy': ('Y', 'mm'),
'dz': ('Z', 'mm'),
'sigxx': ('smooth_stress_xx', 'MPa'),
'sigyy': ('smooth_stress_yy', 'MPa'),
'sigzz': ('smooth_stress_zz', 'MPa'),
'sigxy': ('smooth_stress_xy', 'MPa'),
'sigxz': ('smooth_stress_xz', 'MPa'),
'sigyz': ('smooth_stress_yz', 'MPa'),
'epsxx': ('smooth_strain_xx', ''),
'epsyy': ('smooth_strain_yy', ''),
'epszz': ('smooth_strain_zz', ''),
'epsxy': ('smooth_strain_xy', ''),
'epsxz': ('smooth_strain_xz', ''),
'epsyz': ('smooth_strain_yz', ''),
'sigvm': ('smooth_stress_tensor', 'MPa') # post.probe computes sigvm
}
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}.tsv'
grpname = extr['name']
v = []
units = ""
if var.lower() in var2field:
# extract values from gmsh views
field = var2field[var.lower()][0]
units = var2field[var.lower()][1]
v = np.array(post.probe(field, grpname))
elif var.lower() in ['x0', 'y0', 'z0']:
# extract coordinates of the nodes
units = 'mm'
v = np.array(post.probe('coords', grpname))
v = np.reshape(v, (v.size//3, 3))
v = v[:, ['x0', 'y0', 'z0'].index(var.lower())]
elif var.lower() in ['fx', 'fy', 'fz']:
# extract forces
units = 'N'
v = np.array(post.probe('force_vector', grpname))
v = np.reshape(v, (v.size//3, 3))
v = v[:, ['fx', 'fy', 'fz'].index(var.lower())]
else:
print(f'unknown variable: "{var}"')
if len(v) > 0:
np.savetxt(fname, v, delimiter='\t', header=var.lower())
if len(v) == 1:
print(f'\t{grpname}_{var} = {v[0]:.3f} {units}')
else:
print(f'\t{grpname}_{var} = ({len(v)} values in "{fname}")')
end_time = time.perf_counter()
elapsed_time = end_time - start_time
print(f'Total Execution time: {elapsed_time:.1f}s')
args = fem.parseargs()
if not args.nogui:
post.view()
# ------------------------------------------------------------------------------
# utility functions involing Gmsh objects
def import_mesh(filename):
"""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
print(f'loading file {plyfile} into Gmsh...')
# gmsh.model.add("fossil")
gmsh.merge(plyfile)
gmsh.model.geo.addSurfaceLoop(surfaceTags=[1], tag=1)
gmsh.model.geo.addVolume(shellTags=[1], tag=1)
gmsh.model.geo.addPhysicalGroup(dim=2, tags=[1], tag=1)
gmsh.model.setPhysicalName(dim=2, tag=1, name='surf')
gmsh.model.geo.addPhysicalGroup(dim=3, tags=[1], tag=1)
gmsh.model.setPhysicalName(dim=3, tag=1, name='all')
elif ext in ['.geo', '.msh']:
gmsh.merge(mandible_gmsh)
else:
raise Exception(
f'Unknown extension: {ext}, please use .ply, .stl, .msh or .geo')
# generate mesh
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3)
# gmsh.fltk.run()
def create_group(mesh_file_or_coord, all_no, all_coords, 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
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)
# print(f'ntags={ntags}')
if gname is None:
gname, _ = os.path.splitext(os.path.basename(mesh_file_or_coord))
# create a Gmsh group
# => a discrete entity with "points" elements
print(f'adding gmsh group "{gname}"...')
etag = gmsh.model.addDiscreteEntity(dim=0)
point_type = gmsh.model.mesh.getElementType("Point", 1)
# print(f'point_type={point_type}')
gmsh.model.mesh.addElementsByType(etag, point_type, [], ntags)
ptag = gmsh.model.geo.addPhysicalGroup(dim=0, tags=[etag])
gmsh.model.setPhysicalName(dim=0, tag=ptag, name=gname)
groups[gname] = (0, ptag)
# required for the nodes to be linked to the physical group
gmsh.model.geo.synchronize()
nods, pos = gmsh.model.mesh.getNodesForPhysicalGroup(*groups[gname])
# print(f'nods={nods}')
# print(f'pos={pos}')
print(f'\tgroup {groups[gname]} ({gname}) '
f'created from {mesh_file_or_coord} ({len(nods)} 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