add new way of adding fixations

models/bonemodel2.py

@@ -18,16 +18,11 @@ def parms(d={}):
     """
     p = {}
 
-    # Geomagic/MeshLab input files
+    # Geomagic/MeshLab input folder
     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['muscles'] = [
         {
@@ -43,12 +38,30 @@ def parms(d={}):
             'method': 'T'
         }
     ]
-    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']
-    }
+
+    # 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']
+        }
+    ]
 
     # material properties
     p['Young'] = 17000.      # [MPa] elastic modulus - bone: 17-20 GPa
@@ -91,10 +104,6 @@ def solve(p={}):
     nods_pos = np.reshape(nods_pos, (len(nods_no), 3))
     # print(f'nods_pos={nods_pos}')
 
-    # load other surface meshes  => groups['axis_pt1'], ...
-    for meshname in ['axis_pt1', 'axis_pt2']:
-        create_group(p[meshname], nods_no, nods_pos, groups, meshname)
-
     # load muscle groups
     mgroups = {}                    # stores muscle group data and loads
     for muscle in p['muscles']:
@@ -121,10 +130,9 @@ def solve(p={}):
             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")   
 
-    if p['food'] == 'fixteeth':
-        # teeth are fixed along chosen directions
-        create_group(p['contact_pts'], nods_no,
-                     nods_pos, groups, 'contact_pts')
+    # create groups of fixed nodes
+    for fix in p['fixations']:
+        create_group(fix['nodes'], nods_no, nods_pos, groups, fix['name'])
 
 
     # --------------------------------------------------------------------------
@@ -136,14 +144,14 @@ def solve(p={}):
     # material
     do_not_delete.append( fem.Medium(pbl, "all", E=p['Young'], nu=p['Poisson']) )
 
+    # add Dirichlet conditions (fixations)
     # print('== ADDING DIRICHLETs...')
-    # axis of rotation
-    for gname in ['axis_pt1', 'axis_pt2']:
-        for d in p['fixations'][gname]:
-            do_not_delete.append(fem.Dirichlet(pbl, gname, dof=d.upper(), val=0.0) )
+    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...')
-    # 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):
@@ -151,20 +159,9 @@ def solve(p={}):
             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]) )
 
-    # type of food
-
-    if p['food'] == 'fixteeth':
-        # teeth are fixed along chosen directions
-        for d in p['fixations']['contact_pts']:
-            do_not_delete.append(fem.Dirichlet(pbl, 'contact_pts', dof=d.upper(), val=0.0) )
-    else:
-        raise Exception('food = "fixteeth" only supported')
-
-
     # python only!
     # do_not_delete.append(fem.Dirichlet(pbl, "all", dof='T', val=0.0) )
 
-    
     # Time integration scheme    
     # print (pbl)
     # print('SOLVE.......')

models/dolicorhynchops/dolicorhynchops_10k.py

@@ -9,10 +9,6 @@ def parms(d={}):
     import os
     path = os.path.join(os.path.dirname(__file__),'10k')
     p['bone'] = f'{path}/mandible.stl'
-    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['muscles'] = [
         {
             'file': f'{path}/Lmuscle.stl',
@@ -28,11 +24,24 @@ def parms(d={}):
             'method': 'T'
         }
     ]
-    p['fixations'] = {
-        'contact_pts': ['x'],
-        'axis_pt1': ['x', 'y', 'z'],
-        'axis_pt2': ['x', 'z']
-    }
+    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']
+        }
+    ]
 
     # material properties
     p['density'] = 1.850e-9  # [T/mm³]