new loads

️ More versatile fixation system

It is now possible to define more complex sets of fixations.

Previously, 3 sets of nodes were hard-coded in the input file as "contact_pts", "axis_pt1", "axis_pt2".

Example:

    p['contact_pts'] = [[-10.20362, -17.46838, -229.9061],
                        [-11.92466, 26.3042, -229.5354]]
    p['axis_pt1'] = f'{path}/LTMJ.off'
    p['axis_pt2'] = [-8.716309, 79.13171, 233.8385]
    p['fixations'] = {
        'contact_pts': ['x'],
        'axis_pt1': ['x', 'y', 'z'],
        'axis_pt2': ['x', 'z']
    }

In this new version, the user can define any number of fixation sets and name them as he/she wants. It is thus possible to define much more complex senarios than before.

The previous commands should be gathered in a single list of fixations:

    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']
        }
    ]

This list of fixations can be of any size.

️ In consequence, the old syntax of input files is not compatible with this new version of fossils! (but translation is rather straightforward).

️ Additional loads

Additional loads can be added with a syntax similar to the new syntax for fixations:

    p['loads'] = [
        {
            'name': 'force1',
            'nodes': [-10.20362, -17.46838, -229.9061],
            'values': [0., 0., -100.]
        },
        {
            'name': 'force2',
            'nodes': [-11.92466, 26.3042, -229.5354],
            'values': [0., 0., 100.]
        }      
    ]

with:

• name: a user-defined string, mainly used to name the group of nodes,
• nodes: as for fixations, it can be a point defined by its 3 coordinates between square brackets ([x, y, z]), or a list of points defined in the same manner (e.g. [[x1, y1, z1], [x2, y2, z2], ...]) or a mesh file containing one or more vertices.
• values: the force vector applied to each node.

A new model named models/dolicorhynchops/dolicorhynchops_10k_extloads.py has been created to test this new feature.

️ New load distribution models on surface meshes

In addition to the "Ad Hoc Uniform Traction Model" ('U'), the "Tangential-Traction Model" (T) and the "Tangential-Plus-Normal-Traction Model" (T+N), the user can now use 2 new load distributions:

Directional-Traction Model (D)

The forces are distributed on the nodes of a given surface (file key). They are all pointing in the focal_pt direction seen as a vector. They are scaled so that the total resulting force has the magnitude prescribed by the user (with the key force)

Example:

    p['muscles'] = [
        {
            'method': 'D',
            'file': f'{path}/Lmuscle.stl',  # surface mesh
            'force': 100.,                  # total force
            'focalpt': [-1, 1, 1],          # <= used as a direction
        }
    ]

Note: the keyword focal_pt can be misleading here. The vector given as argument is a direction. However we have kept the string focal_pt as a keyword to have the same keywords for all the methods.

:dart:Normal-Traction Model (N)

In this new method, the nodal forces are applied in the perpendicular direction to each facet of the surface mesh. The focalpt keyword is ignored. Forces are scaled to have the prescribed resulting magnitude (the sum of all the forces is force).

Example:

    p['muscles'] = [
        {
            'method': 'N',
            'file': f'{path}/Lmuscle.stl',  # surface mesh
            'force': 100.,                  # total force
            'focalpt': [0, 0, 0],           # <= ignored
        }
    ]

models/bonemodel.py

@@ -18,23 +18,18 @@ def parms(d={}):
     """
     p = {}
 
-    # Geomagic/MeshLab input files
+    # folder containing meshes
     path = 'dolicorhynchops/10k'
+
     # mandible mesh (.ply/.geo/.msh)
-    p['bone'] = f'{path}/mandible.ply' 
-    # one or several vertices (.ply/.off) (used if p['food']=='fixteeth')
-    p['contact_pts'] = f'{path}/teeth.off'
-    # left temporomandibular joint - 1 vertex (.ply/.off/coordinates)
-    p['axis_pt1'] = [0.1458893, -73.13895, 227.3909]
-    # right temporomandibular joint - 1 vertex (.ply/.off/coordinates) 
-    p['axis_pt2'] = f'{path}/RTMJ.off'
+    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'
+            'method': 'T' # loading model: 'U', 'T', 'T+N', 'D' or 'N'
         },
         { 
             'file': f'{path}/Rmuscle.off', 
@@ -43,12 +38,35 @@ def parms(d={}):
             '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['food'] = 'fixteeth'  # type of food: 'fixteeth' or 'cylinder'
-    p['fixations'] = {  
-        'contact_pts': ['x','y','z'],  # constrained degrees of freedom
-        'axis_pt1': ['x','y','z'],
-        'axis_pt2': ['x','y','z']
-    }
+    
     # rigid cylinder (used if p['food']=='cylinder') 
     R = 20
     p['cylinder'] = {
@@ -90,13 +108,12 @@ def getMetafor(p={}):
     # load main mandible mesh (volume or surface - it will be meshed in 3D)
     import_mesh(domain, p['bone'], p['use_gmshOld'])
 
-    # define groups
+    # -- 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']
@@ -106,14 +123,10 @@ def getMetafor(p={}):
     groups['surf'].addMeshPointsFromObject(groups['all'], BoundarySelector())
     nods_no, nods_pos = extract_nodes_from_group(groups['surf'])
 
-    # load other surface meshes  => groups['axis_pt1'], ...
-    for meshname in ['axis_pt1', 'axis_pt2']:
-        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
+        # 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)
@@ -128,7 +141,11 @@ def getMetafor(p={}):
             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'])
+
     # --------------------------------------------------------------------------
 
     # material 
@@ -152,13 +169,13 @@ def getMetafor(p={}):
     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 ['axis_pt1', 'axis_pt2']:
-        for d in p['fixations'][gname]:
-            domain.getLoadingSet().define(groups[gname], Field1D(txt2key[d],RE), 0.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)}')
@@ -169,6 +186,16 @@ 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)
 
+    # 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':
@@ -193,10 +220,7 @@ def getMetafor(p={}):
         ci.addProperty(prp2)
         domain.getInteractionSet().add(ci)
     elif p['food']=='fixteeth':
-        # teeth are fixed along X
-        create_group(p['contact_pts'], nods_no, nods_pos, domain, groups, 'contact_pts')
-        for d in p['fixations']['contact_pts']:
-            domain.getLoadingSet().define(groups['contact_pts'], Field1D(txt2key[d],RE), 0.0)
+        pass
     else:
         # p['food'] is a filename
         ply2mtf(domain, p['food'])