change metafor model

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

models/bonemodel2.py

@@ -18,7 +18,7 @@ def parms(d={}):
     """
     p = {}
 
-    # Geomagic/MeshLab input folder
+    # folder containing meshes
     path = 'dolicorhynchops/10k'
 
     # mandible mesh (.ply/.geo/.msh)
@@ -29,7 +29,7 @@ def parms(d={}):
             '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',
@@ -63,7 +63,7 @@ def parms(d={}):
         }
     ]
 
-    p['loads'] = []
+    p['loads'] = [] # additional loads
 
     # material properties
     p['Young'] = 17000.      # [MPa] elastic modulus - bone: 17-20 GPa
@@ -80,15 +80,15 @@ def solve(p={}):
 
     p = parms(p)
 
-    # import json
-    # print(f'parameters={json.dumps(p, indent=4, sort_keys=True)}')
+    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
+    # -- define groups
 
     # group (3,1)='all' is the group containing all volume elements
     groups = {}
@@ -97,19 +97,18 @@ def solve(p={}):
     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
+    #   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))
-    # print(f'nods_pos={nods_pos}')
 
     # 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'])