mandiblemodel.py

# -*- coding: utf-8 -*-
# mandible model

from wrap import *
import boneload
import os

# enable full parallelism
StrVectorBase.useTBB()
StrMatrixBase.useTBB()
ContactInteraction.useTBB()

def parms(d={}):
    """default parameters of the model
    """
    p = {}

    # 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) (used if p['food']=='fixteeth')
    p['LTMJ'] = [0.1458893, -73.13895, 227.3909] # left temporomandibular joint - 1 vertex (.ply/.off/coordinates)
    p['RTMJ'] = f'{path}/RTMJ.off'      # right temporomandibular joint - 1 vertex (.ply/.off/coordinates) 
    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'] = 'fixteeth'  # type of food: 'fixteeth' or 'cylinder'
    p['fixations'] = {  
        'teeth': ['x','y','z'],  # constrained degrees of freedom
        'LTMJ': ['x','y','z'],
        'RTMJ': ['x','y','z']
    }
    # 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.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['mandible'])

    # 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 other surface meshes  => groups['LTMJ'], ...
    for meshname in ['LTMJ', 'RTMJ']:
        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
        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)
        
    # ------------------------------------------------------------------

    # 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)

    txt2key = { 'x': TX, 'y':TY, 'z':TZ }

    # axis of rotation
    for gname in ['LTMJ', 'RTMJ']:
        for d in p['fixations'][gname]:
            domain.getLoadingSet().define(groups[gname], Field1D(txt2key[d],RE), 0.0)

    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)

    # 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 p['fixations']['teeth']:
            domain.getLoadingSet().define(groups['teeth'], Field1D(txt2key[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, p['dt0'])
    tsm.setNextTime(1.0, 1, 1.0)

    mim = metafor.getMechanicalIterationManager()
    mim.setResidualTolerance(p['tolNR'])

    # parallel solver
    solman = metafor.getSolverManager()
    solver = DSSolver()
    solman.setSolver(solver)

    # history curves
    valuesmanager = metafor.getValuesManager()
    idx=1
    valuesmanager.add(idx, MiscValueExtractor(metafor,EXT_T),'time'); idx+=1
    for name, group in groups.items():
        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)
            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
    
    return metafor

# ------------------------------------------------------------------------------
# utility functions involing Metafor objects


def import_mesh(domain, 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
        # 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}')
    # from toolbox.gmsh import GmshImport # FIXME: does not work!
    from toolbox.gmshOld import GmshImport
    f = os.path.join(os.path.dirname(__file__), mandible_gmsh)
    mesher = GmshImport(f, domain)
    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
      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'nodes={nodes}')
        print(f'mesh_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))
    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'group #{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 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')


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