add first contact model - to be continued

models/contact/cylinder.geo

+// gmsh file 
+// used to create a meshed cylindrical surface
+//  .load it in gmsh
+//  .generate mesh with Mesh/2D
+//  .export mesh in .stl format 
+//  .stl can be read by meshlab, cleaned and converted to ply
+//
+// using the cmd line:
+//      gmsh -2 -format stl -o cylinder.stl cylinder.geo 
+//
+SetFactory("OpenCASCADE");
+
+R=20;                   // <= radius of the cylinder
+Ly=100;                 // <= length of the cylinder
+
+Cylinder(1) = {-R-8, -Ly/2, -150, 0, Ly, 0, R, 2*Pi};
+
+
+Field[1] = MathEval;
+Field[1].F = "4.0";     // <= length of the edges of the triangles 
+Background Field = 1;
+
+

models/dolicorhynchops/10k/LmuscleF.off

 NOFF
 1 0 0
--100.1458893 -173.13895 227.3909 -0.593782 0.162079 0.788133 
+-100.1458893 -173.13895 227.3909

models/dolicorhynchops_10k_contact.py

+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Dolicorhynchops osborni FHSM VP404
+#   10k faces on the mandible surface
+
+
+def parms(d={}):
+    p = {}
+    path = 'dolicorhynchops/10k'
+    p['mandible'] = f'{path}/mandible.ply'
+    p['teeth'] = f'{path}/teeth.off'
+    p['LTMJ'] = f'{path}/LTMJ.off'
+    p['RTMJ'] = f'{path}/RTMJ.off'
+    p['muscles'] = [
+        {
+            'file': f'{path}/Lmuscle.ply',
+            'force': 100.,
+            'focalpt': f'{path}/LmuscleF.off',
+            'method': 'T'                              # 'U', 'T', 'T+N'
+        },
+        {
+            'file': f'{path}/Rmuscle.off',
+            'force': 100.,
+            'focalpt': f'{path}/RmuscleF.off',
+            'method': 'T'
+        }
+    ]
+    p['food'] = 'contact/cylinder.ply'
+    p.update(d)
+    return p
+
+
+def getMetafor(p={}):
+    import mandiblemodel as model
+    return model.getMetafor(parms(p))

models/mandiblemodel.py

@@ -25,6 +25,7 @@ def parms(d={}):
         { '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
 
     # material properties
     p['density'] = 1.850e-9  # [T/mm³] - bone: 1.850 kg/l
@@ -53,6 +54,9 @@ 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
@@ -260,3 +264,36 @@ class MuscleGroup:
         self.tris = tris 
         self.ntags = ntags
         self.loads = loads
+
+
+def ply2mtf(domain, meshfile):
+    """load a surface into Metafor's memory from a ply file
+    """
+    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')