diff --git a/models/bonemodel.py b/models/bonemodel.py
index d2624db217d9e0ab9ce0345e78de724ac934126c..00d1c7acbeed271bc9f1346e22175ce0b007bf4b 100644
--- a/models/bonemodel.py
+++ b/models/bonemodel.py
@@ -18,23 +18,18 @@ def parms(d={}):
"""
p = {}
- # Geomagic/MeshLab input files
+ # folder containing meshes
path = 'dolicorhynchops/10k'
+
# mandible mesh (.ply/.geo/.msh)
- p['bone'] = f'{path}/mandible.ply'
- # one or several vertices (.ply/.off) (used if p['food']=='fixteeth')
- p['contact_pts'] = f'{path}/teeth.off'
- # left temporomandibular joint - 1 vertex (.ply/.off/coordinates)
- p['axis_pt1'] = [0.1458893, -73.13895, 227.3909]
- # right temporomandibular joint - 1 vertex (.ply/.off/coordinates)
- p['axis_pt2'] = f'{path}/RTMJ.off'
+ 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'
+ 'method': 'T' # loading model: 'U', 'T', 'T+N', 'D' or 'N'
},
{
'file': f'{path}/Rmuscle.off',
@@ -43,12 +38,35 @@ def parms(d={}):
'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['food'] = 'fixteeth' # type of food: 'fixteeth' or 'cylinder'
- p['fixations'] = {
- 'contact_pts': ['x','y','z'], # constrained degrees of freedom
- 'axis_pt1': ['x','y','z'],
- 'axis_pt2': ['x','y','z']
- }
+
# rigid cylinder (used if p['food']=='cylinder')
R = 20
p['cylinder'] = {
@@ -90,13 +108,12 @@ def getMetafor(p={}):
# load main mandible mesh (volume or surface - it will be meshed in 3D)
import_mesh(domain, p['bone'], p['use_gmshOld'])
- # define groups
+ # -- 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']
@@ -106,14 +123,10 @@ def getMetafor(p={}):
groups['surf'].addMeshPointsFromObject(groups['all'], BoundarySelector())
nods_no, nods_pos = extract_nodes_from_group(groups['surf'])
- # load other surface meshes => groups['axis_pt1'], ...
- for meshname in ['axis_pt1', 'axis_pt2']:
- create_group(p[meshname], nods_no, nods_pos, domain, groups, meshname)
-
# load muscle groups
mgroups = {} # stores muscle group data and loads
for muscle in p['muscles']:
- # load focal point
+ # 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)
@@ -128,7 +141,11 @@ def getMetafor(p={}):
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'])
+
# --------------------------------------------------------------------------
# material
@@ -152,13 +169,13 @@ def getMetafor(p={}):
fct.setData(0.0, 0.0)
fct.setData(1.0, 1.0)
- txt2key = { 'x': TX, 'y':TY, 'z':TZ }
-
- # axis of rotation
- for gname in ['axis_pt1', 'axis_pt2']:
- for d in p['fixations'][gname]:
- domain.getLoadingSet().define(groups[gname], Field1D(txt2key[d],RE), 0.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)}')
@@ -169,6 +186,16 @@ def getMetafor(p={}):
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':
@@ -193,10 +220,7 @@ def getMetafor(p={}):
ci.addProperty(prp2)
domain.getInteractionSet().add(ci)
elif p['food']=='fixteeth':
- # teeth are fixed along X
- create_group(p['contact_pts'], nods_no, nods_pos, domain, groups, 'contact_pts')
- for d in p['fixations']['contact_pts']:
- domain.getLoadingSet().define(groups['contact_pts'], Field1D(txt2key[d],RE), 0.0)
+ pass
else:
# p['food'] is a filename
ply2mtf(domain, p['food'])
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
index 3e62d531e80a855a030a699c6b44dad57cc16103..1bc8373ee05dd0d72b6641093eb99c7018c4574d 100644
--- a/models/bonemodel2.py
+++ b/models/bonemodel2.py
@@ -18,7 +18,7 @@ def parms(d={}):
"""
p = {}
- # Geomagic/MeshLab input folder
+ # folder containing meshes
path = 'dolicorhynchops/10k'
# mandible mesh (.ply/.geo/.msh)
@@ -29,7 +29,7 @@ def parms(d={}):
'file': f'{path}/Lmuscle.ply',
'force': 100.,
'focalpt': f'{path}/LmuscleF.off',
- 'method': 'T' # loading model: 'U', 'T', 'T+N'
+ 'method': 'T' # loading model: 'U', 'T', 'T+N', 'D' or 'N'
},
{
'file': f'{path}/Rmuscle.off',
@@ -63,7 +63,7 @@ def parms(d={}):
}
]
- p['loads'] = []
+ p['loads'] = [] # additional loads
# material properties
p['Young'] = 17000. # [MPa] elastic modulus - bone: 17-20 GPa
@@ -80,15 +80,15 @@ def solve(p={}):
p = parms(p)
- # import json
- # print(f'parameters={json.dumps(p, indent=4, sort_keys=True)}')
+ 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
+ # -- define groups
# group (3,1)='all' is the group containing all volume elements
groups = {}
@@ -97,19 +97,18 @@ def solve(p={}):
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
+ # 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))
- # print(f'nods_pos={nods_pos}')
# load muscle groups
mgroups = {} # stores muscle group data and loads
for muscle in p['muscles']:
- # load focal point
+ # load focal point / loading direction
if isinstance(muscle['focalpt'], str):
fullpath = os.path.join(
os.path.dirname(__file__), muscle['focalpt'])