tangent+normal loads OK

models/boneload.py

@@ -191,7 +191,91 @@ def scale_loads(loads, total_force):
     return loads
 
 
-def create_uniform_loads(nodes, tris, total_force, target=None):
+def compute_normal_force_factor(muscle, centre, nunit, targetP, debug=False, surface=None, muscleF=None):
+
+    if debug:
+        print(f'centre={centre}')
+        print(f'nunit={nunit}')
+        plane = Plane(centre, centre+nunit, targetP)  # only for display
+    
+    nclip = -(targetP - (centre+(targetP*nunit)*nunit)).normalized()
+    planepart = ClipPolyData(muscle.polydata, centre, nclip)
+    planecut = PlaneCut(planepart.clip.GetOutput(), centre, nunit.cross(targetP-centre))
+    closest = ClosestPart(planecut.cutter.GetOutputPort(), centre)
+    pts, ids = get_segments(closest.cfilter.GetOutput())
+    
+    # compute "s(r)"
+    s = compute_length(pts, ids)
+    if debug:  
+        print(f's(r)={s}')
+
+    # front part of the cut
+    planepart2 = ClipPolyData(muscle.polydata, centre, -nclip)
+    planecut2 = PlaneCut(planepart2.clip.GetOutput(), centre, nunit.cross(targetP-centre))
+    closest2 = ClosestPart(planecut2.cutter.GetOutputPort(), centre)
+    pts2, ids2 = get_segments(closest2.cfilter.GetOutput())
+
+    sorted_v = sort_vertices(pts, ids, centre)
+    sorted_v2 = sort_vertices(pts2, ids2, centre)
+
+    fpath = build_clean_path(pts, sorted_v)
+    fpath2 = build_clean_path(pts2, sorted_v2)
+
+    curv, polyx, polyy = compute_curvature(fpath, fpath2, nunit)
+    if debug:
+        print(f'curv={curv}')
+
+    # display things
+
+    # project vertices of fpath onto the plane for display
+    if debug:
+        p0, xaxis, yaxis = build_local_axes(fpath, nunit)
+        xpath = [ (p-p0)*xaxis for p in reversed(fpath2) ] + [ (p-p0)*xaxis for p in fpath ]
+        ypath = [ (p-p0)*yaxis for p in reversed(fpath2) ] + [ (p-p0)*yaxis for p in fpath ]
+
+        import matplotlib.pyplot as plt
+        import numpy as np        
+        plt.plot(xpath, ypath, 'o-', label='path')
+        plt.plot(polyx, polyy, 'r-', label='polynomial')
+
+        plt.xlabel('x')
+        plt.ylabel('y')
+        plt.title('muscle fiber shape')
+        plt.grid(True)
+        plt.legend()
+        plt.show()
+
+        # 3D debugging view
+        view = View()
+
+        view.addActors(surface.actor) 
+        view.addActors(muscle.actor)  
+        view.addActors(muscleF.actor)
+        view.addActors(plane.actor)  
+        view.addActors(planepart.actor) 
+        view.addActors(planecut.actor) 
+        view.addActors(closest.actor) 
+        view.addActors(closest2.actor) 
+
+        surface.actor.GetProperty().SetOpacity(0.1)
+
+        muscle.mapper.SetResolveCoincidentTopologyToPolygonOffset()
+        muscle.actor.GetProperty().SetColor( colors.GetColor3d('red') )
+        muscle.actor.GetProperty().SetOpacity(0.5)
+        muscle.actor.GetProperty().EdgeVisibilityOn()
+        muscle.actor.GetProperty().SetPointSize(1.0)
+
+        muscleF.actor.GetProperty().SetPointSize(15.0)
+        muscleF.actor.GetProperty().SetColor( colors.GetColor3d('green') )
+
+        closest2.actor.GetProperty().SetColor( colors.GetColor3d('pink') )
+
+        view.interact()
+
+    return s, curv
+
+
+def create_loads(nodes, tris, total_force, target=None, normal_comp=False):
     """ Ad Hoc Uniform Traction Model (Grosse et al. 2007)
 
     "In this model, we applied uniform traction to the surfaces
@@ -205,6 +289,9 @@ def create_uniform_loads(nodes, tris, total_force, target=None):
     note: if target==None: the traction is normal to the facet
     """
     project=True
+    if normal_comp:
+        muscle = SurfMesh(nodes, tris, vertices=False)
+
     loads = [Pt([0,0,0]) for n in nodes ]
 
     if target:
@@ -221,19 +308,28 @@ def create_uniform_loads(nodes, tris, total_force, target=None):
             direct = targetP-centre
         else:
             direct = n 
-        direct.normalize()
+        direct.normalize()          # direction of the traction (normalised)
         n = e1.cross(e2)            # normal vector
         area2 = abs(n)              # 2 x area
         area = area2/2              # triangle area
         total_area += area
         nunit = n/area2             # unit normal vector
         
-        force = (area/3)*direct      # compute nodal force
-
+        # we use a unit traction which will be scaled later
+        traction = direct
         if project:
-            ps = force*nunit
+            ps = traction*nunit
             if ps<0.0: # no line of sight => projection onto the element plane
-                force = force - ps*nunit
+                traction = traction - ps*nunit
+                # tangential traction should be normalised if we want to follow
+                # Grosse implementation
+                traction.normalize() # could fail if direct and nunit are aligned in opposite directions
+
+                # we can add normal component here.
+                s, curv = compute_normal_force_factor(muscle, centre, nunit, targetP)
+                traction = traction - (s*curv)*nunit
+
+        force = (area/3)*traction      # compute nodal force
 
         for i in range(3):          # 3 nodes per triangle
             loads[tri[i]]+=force    # assembly
@@ -465,6 +561,7 @@ def sort_segments(pts, ids):
 
     return sorted_ids
 
+
 def sort_vertices(pts, ids, centre):
     """build a list of vertex indices of an open piecewise linear path.
     the points are sorted by increasing distance to "centre" (Pt object)
@@ -528,6 +625,7 @@ def build_clean_path(pts, sorted_v):
             break
     return fpath
 
+
 def build_local_axes(fpath, nunit):
     p0 = fpath[0]     # origin
     p1 = fpath[1]     # next point
@@ -535,6 +633,7 @@ def build_local_axes(fpath, nunit):
     yaxis = nunit 
     return p0, xaxis, yaxis   
 
+
 def compute_curvature(fpath, fpath2, nunit):
 
     # fit a polynomial
@@ -561,10 +660,6 @@ def compute_curvature(fpath, fpath2, nunit):
     return curv, polyx, polyy
 
 
-
-
-
-
 # ------------------------------------------------------------------------------
 # Visualisation
 

models/mandiblemodel.py

@@ -81,7 +81,7 @@ def getMetafor(p={}):
         name, nodes, tris, ntags = \
             create_group(muscle['file'], nods_no, nods_pos, domain, groups)
         # create loads on the surface
-        loads = boneload.create_uniform_loads(nodes, tris, \
+        loads = boneload.create_loads(nodes, tris, \
             total_force=muscle['force'], target=focalnodes[0])
         # store data in a structure
         mgroups[name] = MuscleGroup(nodes, tris, ntags, loads)

models/view_cylinder.py

@@ -13,7 +13,7 @@ if __name__=="__main__":
     # create vector field
     focalnodes, _ = load_msh('loadings/muscleF.off')
     nodes, tris = load_msh('loadings/muscle.ply')
-    loads = create_uniform_loads(nodes, tris, \
+    loads = create_loads(nodes, tris, \
         total_force=10., target=focalnodes[0])
 
     muscle.add_vectors(loads)

models/view_cylinder2.py

@@ -17,13 +17,12 @@ if __name__=="__main__":
     # create vector field
     focalnodes, _ = load_msh('loadings/muscleF.off')
     nodes, tris = load_msh('loadings/muscle.ply')
-    # loads = create_uniform_loads(nodes, tris, \
+    # loads = create_loads(nodes, tris, \
     #     total_force=10., target=focalnodes[0])
 
     # muscle.add_vectors(loads)
     # forces = Arrows(muscle.polydata)
-    targetP = Pt(focalnodes[0]  )
-    total_area=0.
+    targetP = Pt( focalnodes[0] )
 
     # notri = 140
     # notri = 1
@@ -38,72 +37,28 @@ if __name__=="__main__":
         direct = targetP-centre
         direct.normalize()
         n = e1.cross(e2)            # normal vector
-        area2 = abs(n)
-        area = area2/2              # triangle area
-        total_area+=area
-        nunit = n/area2             # unit normal vector
+        nunit = n.normalized()      # unit normal vector
         
-        force = (area/3)*direct      # compute nodal force
+        compute_normal_force_factor(muscle, centre, nunit, targetP, \
+            debug=True, surface=surface, muscleF=muscleF)
 
+    # compute all forces
 
-        plane = Plane(centre, centre+nunit, targetP) 
-        nclip = -(targetP - (centre+(targetP*nunit)*nunit)).normalized()
-        planepart = ClipPolyData(muscle.polydata, centre, nclip)
-        planecut = PlaneCut(planepart.clip.GetOutput(), centre, nunit.cross(targetP-centre))
-        closest = ClosestPart(planecut.cutter.GetOutputPort(), centre)
-        pts, ids = get_segments(closest.cfilter.GetOutput())
-        
-        # compute "s(r)"
-        s = compute_length(pts, ids)  
-        print(f's(r)={s}')
-
-        # front part of the cut
-        planepart2 = ClipPolyData(muscle.polydata, centre, -nclip)
-        planecut2 = PlaneCut(planepart2.clip.GetOutput(), centre, nunit.cross(targetP-centre))
-        closest2 = ClosestPart(planecut2.cutter.GetOutputPort(), centre)
-        pts2, ids2 = get_segments(closest2.cfilter.GetOutput())
-
-        sorted_v = sort_vertices(pts, ids, centre)
-        sorted_v2 = sort_vertices(pts2, ids2, centre)
-
-        fpath = build_clean_path(pts, sorted_v)
-        fpath2 = build_clean_path(pts2, sorted_v2)
+    loads = create_loads(nodes, tris, \
+                total_force=10., target=focalnodes[0], normal_comp=True)
 
-        curv, polyx, polyy = compute_curvature(fpath, fpath2, nunit)
-        print(f'curv={curv}')
+    muscle.add_vectors(loads)
+    forces = Arrows(muscle.polydata)
 
-        # display things
 
-        # project vertices of fpath onto the plane for display
-        p0, xaxis, yaxis = build_local_axes(fpath, nunit)
-        xpath = [ (p-p0)*xaxis for p in reversed(fpath2) ] + [ (p-p0)*xaxis for p in fpath ]
-        ypath = [ (p-p0)*yaxis for p in reversed(fpath2) ] + [ (p-p0)*yaxis for p in fpath ]
 
-        import matplotlib.pyplot as plt
-        import numpy as np        
-        plt.plot(xpath, ypath, 'o-', label='path')
-        plt.plot(polyx, polyy, 'r-', label='polynomial')
-
-        plt.xlabel('x')
-        plt.ylabel('y')
-        plt.title('muscle fiber shape')
-        plt.grid(True)
-        plt.legend()
-        plt.show()
-
-    # -- view
 
     view = View()
 
     view.addActors(surface.actor) 
     view.addActors(muscle.actor)  
-    view.addActors(muscleF.actor)
-    view.addActors(plane.actor)  
-    view.addActors(planepart.actor) 
-    view.addActors(planecut.actor) 
-    view.addActors(closest.actor) 
-    view.addActors(closest2.actor) 
- 
+    view.addActors(muscleF.actor) 
+    view.addActors(forces.actor) 
 
     surface.actor.GetProperty().SetOpacity(0.1)
     
@@ -116,7 +71,4 @@ if __name__=="__main__":
     muscleF.actor.GetProperty().SetPointSize(15.0)
     muscleF.actor.GetProperty().SetColor( colors.GetColor3d('green') )
 
-    closest2.actor.GetProperty().SetColor( colors.GetColor3d('pink') )
-
-    view.interact()
-
+    view.interact()