Skip to content
Snippets Groups Projects
Commit c38d462a authored by Boman Romain's avatar Boman Romain
Browse files

add rigid cylinder

parent c04b93ec
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
models/dolicorhynchops_10k_contact.py
+ 13
− 2
View file @ c38d462a
......@@ -8,7 +8,7 @@ def parms(d={}):
p = {}
path = 'dolicorhynchops/10k'
p['mandible'] = f'{path}/mandible.ply'
p['teeth'] = f'{path}/teeth.off'
# p['teeth'] = f'{path}/teeth.off'
p['LTMJ'] = f'{path}/LTMJ.off'
p['RTMJ'] = f'{path}/RTMJ.off'
p['muscles'] = [
......@@ -25,7 +25,18 @@ def parms(d={}):
'method': 'T'
}
]
p['food'] = 'contact/cylinder.ply'
p['food'] = 'cylinder'
R = 10
p['cylinder'] = {
'origin': [-R-8, 0, -150],
'radius': R,
'width': 100,
'friction': 0.1,
'penalty': 1e2
}
p['tolNR'] = 1e-4 # [-] equilibrium tolerance
p['dt0'] = 2e-2 # [s] time step size
p.update(d)
return p
......
models/mandiblemodel.py
+ 105
− 13
View file @ c38d462a
......@@ -18,14 +18,24 @@ def parms(d={}):
# Geomagic/MeshLab input files
path = 'dolicorhynchops/10k'
p['mandible'] = f'{path}/mandible.ply' # mandible mesh (.ply/.geo/.msh)
p['teeth'] = f'{path}/teeth.off' # one or several vertices (.ply/.off)
p['teeth'] = f'{path}/teeth.off' # one or several vertices (.ply/.off) (used if p['food']=='fixteeth')
p['LTMJ'] = f'{path}/LTMJ.off' # left temporomandibular joint - 1 vertex (.ply/.off)
p['RTMJ'] = f'{path}/RTMJ.off' # right temporomandibular joint - 1 vertex (.ply/.off)
p['muscles'] = [
{ 'file': f'{path}/Lmuscle.ply', 'force': 100., 'focalpt':f'{path}/LmuscleF.off', 'method':'T'},
{ 'file': f'{path}/Rmuscle.off', 'force': 100., 'focalpt':f'{path}/RmuscleF.off', 'method':'T'}
]
p['food'] = None
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
......@@ -34,6 +44,7 @@ def parms(d={}):
# numerical parameters
p['tolNR'] = 1e-6 # [-] equilibrium tolerance
p['dt0'] = 1.0 # [s] time step size
p.update(d)
return p
......@@ -54,9 +65,6 @@ def getMetafor(p={}):
# load main mandible mesh (volume or surface - it will be meshed in 3D)
import_mesh(domain, p['mandible'])
if p['food']:
ply2mtf(domain, p['food'])
# define groups
# group #1 is the group containing all volume elements
......@@ -64,7 +72,7 @@ def getMetafor(p={}):
groups = {}
groups['all'] = groupset(1)
print("the whole mesh has %d nodes" % groups['all'].getNumberOfMeshPoints())
print("the mandible mesh has %d nodes" % groups['all'].getNumberOfMeshPoints())
# extract external surface of the mesh => groups['surf']
# nods_no: tags of the nodes
......@@ -74,7 +82,7 @@ def getMetafor(p={}):
nods_no, nods_pos = extract_nodes_from_group(groups['surf'])
# load other surface meshes => groups['LTMJ'], ...
for meshname in ['LTMJ', 'RTMJ', 'teeth']:
for meshname in ['LTMJ', 'RTMJ']:
create_group(p[meshname], nods_no, nods_pos, domain, groups, meshname)
# load muscle groups
......@@ -115,14 +123,11 @@ def getMetafor(p={}):
fct.setData(0.0, 0.0)
fct.setData(1.0, 1.0)
# axis of rotation
for gname in ['LTMJ', 'RTMJ']:
for d in [ TX, TY, TZ ]:
domain.getLoadingSet().define(groups[gname], Field1D(d,RE), 0.0)
for gname in ['teeth']:
for d in [ TX ]:
domain.getLoadingSet().define(groups[gname], Field1D(d,RE), 0.0)
mshpts = geometry.getMesh().getPointSet()
for name, gr in mgroups.items():
for i,load in enumerate(gr.loads):
......@@ -131,13 +136,45 @@ 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)
# 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':
# teeth are fixed along X
create_group(p['teeth'], nods_no, nods_pos, domain, groups, 'teeth')
for d in [ TX ]:
domain.getLoadingSet().define(groups['teeth'], Field1D(d,RE), 0.0)
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, 1.0)
tsm.setInitialTime(0.0, p['dt0'])
tsm.setNextTime(1.0, 1, 1.0)
mim = metafor.getMechanicalIterationManager()
......@@ -156,7 +193,13 @@ def getMetafor(p={}):
if name in ['all','surf']: continue
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)
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
......@@ -297,3 +340,52 @@ def ply2mtf(domain, meshfile):
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
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment