Skip to content
Snippets Groups Projects
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
bonemodel.py 18.87 KiB
# -*- coding: utf-8 -*-
# Bone FEA model (formely: "mandible model")
#
# main Metafor script with default parameters.
# R.Boman

from wrap import *
import boneload
import os

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

    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
    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['narch'] = 1           # number of archives
    
    # gmsh toolbox
    p['use_gmshOld'] = False  # use old gmsh interface

    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['bone'], p['use_gmshOld'])

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

    # create groups of fixed nodes
    for fix in p['fixations']:
        create_group(fix['nodes'], nods_no, nods_pos, domain, groups, fix['name'])

    # create groups of additional extractors
    for extr in p['extractors']:
        create_group(extr['nodes'], nods_no, nods_pos, domain, groups, extr['name'])

    # --------------------------------------------------------------------------

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

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

    # 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':
        # 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':
        pass
    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, p['narch'], 1.0)

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

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

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

    # history curves
    # by default, extract all forces on all groups
    valuesmanager = metafor.getValuesManager()
    idx=1
    valuesmanager.add(idx, MiscValueExtractor(metafor,EXT_T),'time'); idx+=1
    for name, group in groups.items():
        # skip the main mesh and its surface
        if name in ['all','surf']: continue
        # skip all additional extractors
        if name in [e['name'] for e in p['extractors']]: continue
        # muscles and loads should remain
        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

    # add additional extractors

    # conversion dicts
    var2field_db = {
        'dx': Field1D(TX, RE),
        'dy': Field1D(TY, RE),
        'dz': Field1D(TZ, RE),
        'x0': Field1D(TX, AB),
        'y0': Field1D(TY, AB),
        'z0': Field1D(TZ, AB),
        'fx': Field1D(TX, GF1),
        'fy': Field1D(TY, GF1),
        'fz': Field1D(TZ, GF1)
    }

    var2field_if = { 
        'sigvm': IF_EVMS, 
        'sigxx': IF_SIG_XX, 
        'sigyy': IF_SIG_YY, 
        'sigzz': IF_SIG_ZZ,
        'sigxy': IF_SIG_XY, 
        'sigyz': IF_SIG_YZ, 
        'sigxz': IF_SIG_XZ,
        'epsxx': IF_GL_STRAIN_XX,
        'epsyy': IF_GL_STRAIN_YY,
        'epszz': IF_GL_STRAIN_ZZ,
        'epsxy': IF_GL_STRAIN_XY,
        'epsyz': IF_GL_STRAIN_YZ,
        'epsxz': IF_GL_STRAIN_XZ 
    }

    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}'
            target = groups[extr['name']]
            # nodal values
            if var.lower() in var2field_db:
                field = var2field_db[var.lower()]
                valuesmanager.add(idx, DbNodalValueExtractor(target, field), fname)
            elif var.lower() in var2field_if:
                field = var2field_if[var.lower()]
                valuesmanager.add(idx, IFNodalValueExtractor(target, field), fname)
            else:
                raise Exception(f'Unknown variable: {var}')
            # if the target is a single node, display the extractor [TSC-EXT]
            # otherwise, the user should read the ascii file in the workspace.
            if target.getNumberOfMeshPoints() == 1:
                metafor.getTestSuiteChecker().checkExtractor(idx)
            idx+=1

    return metafor

# ------------------------------------------------------------------------------
# utility functions involving Metafor objects


def import_mesh(domain, filename, use_gmshOld=True):
    """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}')
    f = os.path.join(os.path.dirname(__file__), mandible_gmsh)
    if use_gmshOld:
        from toolbox.gmshOld import GmshImport
        mesher = GmshImport(f, domain)
    else:
        from toolbox.gmsh import GmshImport # FIXME: does not work!
        mesher = GmshImport(f, domain)
        # mesher.verb = True # debug
        mesher.writeMshFile = True
    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 (.off, .ply, .stl)
      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)

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

    if gname in groups:
        raise Exception(f'create_group: group "{gname}" already exists!')

    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'\tgroup #{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 mesh file (ply, off or stl)
    """
    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