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

import cxxfem as fem
from . import boneload
import os
import gmsh
if not gmsh.isInitialized():
    gmsh.initialize()
import numpy as np


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

    # material properties
    p['Young'] = 17000.      # [MPa] elastic modulus - bone: 17-20 GPa
    p['Poisson'] = 0.3       # [-] Poisson's ratio

    p.update(d)
    return p


def solve(p={}):

    import time
    start_time = time.perf_counter()

    p = parms(p)

    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

    # group (3,1)='all' is the group containing all volume elements
    groups = {}
    groups['all'] = (3, 1)
    nodeTags, coord = gmsh.model.mesh.getNodesForPhysicalGroup(*groups['all'])
    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
    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))

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

        # save loads to file
        with open(f'loads_{name}.tsv','w') as f:
            f.write('node#\tfx\tfy\tfz\n')
            for i,load in enumerate(loads):
                f.write(f"{int(ntags[i])}\t{load.x[0]}\t{load.x[1]}\t{load.x[2]}\n")   

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

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

    # --------------------------------------------------------------------------
    # print('== ADDING MEDIUM...')
    print('')

    do_not_delete = []

    # material
    do_not_delete.append( fem.Medium(pbl, "all", E=p['Young'], nu=p['Poisson']) )

    # add Dirichlet conditions (fixations)
    # print('== ADDING DIRICHLETs...')
    for fix in p['fixations']:
        for d in fix['direction']:
            do_not_delete.append(fem.Dirichlet(pbl, fix['name'], dof=d.upper(), val=0.0) )

    # add Neumann conditions (loads) from muscle loads
    # print('== ADDING NEUMANNs...')
    for name, gr in mgroups.items():
        # print(f'len(gr.ntags)={len(gr.ntags)}')
        for i,load in enumerate(gr.loads):
            do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'X', load.x[0]) )  # TODO: regler pbl de cast "gr.ntags[i] => size_t"
            do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'Y', load.x[1]) )  # TODO: regler "dof=" qui marche pas
            do_not_delete.append(fem.Neumann(pbl, int(gr.ntags[i]), 'Z', load.x[2]) )

    # add Neumann conditions (loads) from external loads
    for load in p['loads']:
        name, nodes, tris, ntags = \
            create_group(load['nodes'], nods_no, nods_pos, groups, load['name'])
        for ntag in ntags:
            for i in range(3):
                do_not_delete.append(fem.Neumann(pbl, int(ntag), 'XYZ'[i], load['values'][i]) )


    # python only!
    # do_not_delete.append(fem.Dirichlet(pbl, "all", dof='T', val=0.0) )

    # Time integration scheme    
    # print (pbl)
    # print('SOLVE.......')
    pbl.quadrature = 'Gauss1' # 1 GP/tetra is enough

    solver = fem.Solver(pbl)
    solver.linear_solver = "pardiso"
    solver.solve()

    post = fem.Post(solver)
    post.build_views()
    post.show_views([ 'smooth_stress_tensor' ])
    post.set_camera_3D()
    post.write()
    post.deform()

    print('extracting results...')
    print(f'\tinternal energy = {solver.int_energy:.3f} N.mm')
    # fixations (reaction forces)
    for fix in p['fixations']:
        grpname = fix['name']
        v = np.array(post.probe('force_vector', grpname))
        v = np.reshape(v, (v.size//3, 3))
        v = np.sum(v, axis=0)
        print(f'\t{grpname}_Fx = {v[0]:.3f} N')
        print(f'\t{grpname}_Fy = {v[1]:.3f} N')
        print(f'\t{grpname}_Fz = {v[2]:.3f} N')

    # muscle groups
    for name, group in mgroups.items():
        v=np.array( [0.0, 0.0, 0.0] )
        for tag in group.ntags:
            v += np.array(post.probe('force_vector', int(tag)))
        print(f'\t{name}_Fx = {v[0]:.3f} N')
        print(f'\t{name}_Fy = {v[1]:.3f} N')
        print(f'\t{name}_Fz = {v[2]:.3f} N')

    # additional extractors

    var2field = { 
        'dx': ('X', 'mm'), 
        'dy': ('Y', 'mm'), 
        'dz': ('Z', 'mm'),  
        'sigxx': ('smooth_stress_xx', 'MPa'),  
        'sigyy': ('smooth_stress_yy', 'MPa'), 
        'sigzz': ('smooth_stress_zz', 'MPa'), 
        'sigxy': ('smooth_stress_xy', 'MPa'), 
        'sigxz': ('smooth_stress_xz', 'MPa'), 
        'sigyz': ('smooth_stress_yz', 'MPa'), 
        'epsxx': ('smooth_strain_xx', ''), 
        'epsyy': ('smooth_strain_yy', ''), 
        'epszz': ('smooth_strain_zz', ''), 
        'epsxy': ('smooth_strain_xy', ''), 
        'epsxz': ('smooth_strain_xz', ''), 
        'epsyz': ('smooth_strain_yz', ''),
        'sigvm': ('smooth_stress_tensor', 'MPa')  # post.probe computes sigvm                  
    }

    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}.tsv'
            grpname = extr['name']
            v = []
            units = ""
            if var.lower() in var2field:
                # extract values from gmsh views
                field = var2field[var.lower()][0]
                units = var2field[var.lower()][1]
                v = np.array(post.probe(field, grpname))              
            elif var.lower() in ['x0', 'y0', 'z0']:
                # extract coordinates of the nodes
                units = 'mm'
                v = np.array(post.probe('coords', grpname))
                v = np.reshape(v, (v.size//3, 3))
                v = v[:, ['x0', 'y0', 'z0'].index(var.lower())] 
            elif var.lower() in ['fx', 'fy', 'fz']: 
                # extract forces
                units = 'N'
                v = np.array(post.probe('force_vector', grpname))
                v = np.reshape(v, (v.size//3, 3))
                v = v[:, ['fx', 'fy', 'fz'].index(var.lower())]
            else:
                print(f'unknown variable: "{var}"')

            if len(v) > 0:
                np.savetxt(fname, v, delimiter='\t', header=var.lower())
                if len(v) == 1:
                    print(f'\t{grpname}_{var} = {v[0]:.3f} {units}')
                else:
                    print(f'\t{grpname}_{var} = ({len(v)} values in "{fname}")')                


    end_time = time.perf_counter()
    elapsed_time = end_time - start_time
    print(f'Total Execution time: {elapsed_time:.1f}s')

    args = fem.parseargs()
    if not args.nogui:
        post.view()


# ------------------------------------------------------------------------------
# utility functions involing Gmsh objects


def import_mesh(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
        print(f'loading file {plyfile} into Gmsh...')
        # gmsh.model.add("fossil")
        gmsh.merge(plyfile)
        gmsh.model.geo.addSurfaceLoop(surfaceTags=[1], tag=1)
        gmsh.model.geo.addVolume(shellTags=[1], tag=1)
        gmsh.model.geo.addPhysicalGroup(dim=2, tags=[1], tag=1)
        gmsh.model.setPhysicalName(dim=2, tag=1, name='surf')
        gmsh.model.geo.addPhysicalGroup(dim=3, tags=[1], tag=1)
        gmsh.model.setPhysicalName(dim=3, tag=1, name='all')
    elif ext in ['.geo', '.msh']:
        gmsh.merge(mandible_gmsh)
    else:
        raise Exception(
            f'Unknown extension: {ext}, please use .ply, .stl, .msh or .geo')

    # generate mesh
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(3)
    # gmsh.fltk.run()


def create_group(mesh_file_or_coord, all_no, all_coords, 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'\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)
    # print(f'ntags={ntags}')

    if gname is None:
        gname, _ = os.path.splitext(os.path.basename(mesh_file_or_coord))

    # create a Gmsh group
    #   => a discrete entity with "points" elements
    print(f'adding gmsh group "{gname}"...')
    etag = gmsh.model.addDiscreteEntity(dim=0)
    point_type = gmsh.model.mesh.getElementType("Point", 1)
    # print(f'point_type={point_type}')
    gmsh.model.mesh.addElementsByType(etag, point_type, [], ntags)
    ptag = gmsh.model.geo.addPhysicalGroup(dim=0, tags=[etag])
    gmsh.model.setPhysicalName(dim=0, tag=ptag, name=gname)

    groups[gname] = (0, ptag)

    # required for the nodes to be linked to the physical group
    gmsh.model.geo.synchronize()
    nods, pos = gmsh.model.mesh.getNodesForPhysicalGroup(*groups[gname])
    # print(f'nods={nods}')
    # print(f'pos={pos}')

    print(f'\tgroup {groups[gname]} ({gname}) '
          f'created from {mesh_file_or_coord} ({len(nods)} 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