add external nodal forces as well as external pressures (new "muscle" forces with "D" or "N" code)

cxxfem/src/femSolver.cpp

@@ -67,8 +67,8 @@ void Solver::solve()
 
     create_partition();
     fill_partition();
-    fill_rhs();
     fill_lhs();
+    fill_rhs();
     apply_bc();
     solve_static();
 
@@ -206,11 +206,11 @@ void Solver::fill_partition()
     std::cout << "\telasped = " << timers["fill_partition"].read() << "s\n";
 }
 
-/// fill right-hand side of the system of equations
+/// fill left-hand side of the system of equations
 
-void Solver::fill_rhs()
+void Solver::fill_lhs()
 {
-    timers["fill_rhs"].start();
+    timers["fill_lhs"].start();
 
     std::cout << "assembly of LHS...\n";
 
@@ -527,15 +527,15 @@ void Solver::fill_rhs()
         for (size_t j = 0; j < Kijv[i].size(); ++j)
             delete Kijv[i][j];
 
-    timers["fill_rhs"].stop();
-    std::cout << "\telasped = " << timers["fill_rhs"].read() << "s\n";
+    timers["fill_lhs"].stop();
+    std::cout << "\telasped = " << timers["fill_lhs"].read() << "s\n";
 }
 
 /// compute right-hand side 'f' (sources)
 
-void Solver::fill_lhs()
+void Solver::fill_rhs()
 {
-    timers["fill_lhs"].start();
+    timers["fill_rhs"].start();
 
     std::cout << "assembly of RHS...\n";
 
@@ -627,8 +627,8 @@ void Solver::fill_lhs()
         }
     }
 
-    timers["fill_lhs"].stop();
-    std::cout << "\telasped = " << timers["fill_lhs"].read() << "s\n";
+    timers["fill_rhs"].stop();
+    std::cout << "\telasped = " << timers["fill_rhs"].read() << "s\n";
 
 }
 

models/boneload.py

@@ -461,24 +461,37 @@ def create_loads(nodes, tris, total_force, target, method='T'):
         "adds the normal traction (i.e., pressure) that
         muscle fibers impose on the skull as they wrap around
         it."
+    method = 'D' (Directional-Traction Model)
+        all the nodal forces are applied along a prescribed direction 
+        given by "target" see as a vector
+    method = 'N' (Normal-Traction Model)
+        all the nodal forces are applied along the normal to each triangle
     """
 
     # options
     if method == 'U':
-        project = False
         methodtxt = 'Ad Hoc Uniform Traction Model'
+        project = False
         normal_comp = None
     elif method == 'T':
+        methodtxt = 'Tangential-Traction Model'
         project = True
         normal_comp = False
-        methodtxt = 'Tangential-Traction Model'
     elif method == 'T+N':
+        methodtxt = 'Tangential-Plus-Normal-Traction Model'
         project = True
         normal_comp = True
-        methodtxt = 'Tangential-Plus-Normal-Traction Model'
+    elif method == 'D':
+        methodtxt = 'Directional-Traction Model'
+        project = False
+        normal_comp = None
+    elif method == 'N':
+        methodtxt = 'Normal-Traction Model'
+        project = False
+        normal_comp = None
     else:
         raise Exception(
-            "unknown method: choose from 'uniform', 'tangential' or 'tangential+normal'")
+            f"unknown method ({method}): choose from 'U', 'T', 'T+N', 'D' or 'N'")    
 
     if normal_comp:
         muscle = SurfMesh(nodes, tris, vertices=False)
@@ -492,11 +505,17 @@ def create_loads(nodes, tris, total_force, target, method='T'):
     for tri in tris:
         p1, p2, p3 = [Pt(nodes[tri[i]]) for i in range(3)]
         centre = (p1+p2+p3)/3       # barycentre
-        e1 = p2-p1                  # edge 1
-        e2 = p3-p1                  # edge 2
-        direct = targetP-centre
-        direct.normalize()          # direction of the traction (normalised)
+        e1 = p2 - p1                # edge 1
+        e2 = p3 - p1                # edge 2
         n = e1.cross(e2)            # normal vector
+        if method == 'D':
+            direct = targetP        # force direction is the same for all triangles
+        elif method == 'N':
+            direct = n              # force direction is normal to the triangle
+        else:
+            direct = targetP-centre # force direction varies from 1 triangle to another
+        direct.normalize()          # direction of the traction (normalised)
+
         area2 = abs(n)              # 2 x area
         area = area2/2              # triangle area
         total_area += area

models/bonemodel2.py

@@ -63,6 +63,8 @@ def parms(d={}):
         }
     ]
 
+    p['loads'] = []
+
     # material properties
     p['Young'] = 17000.      # [MPa] elastic modulus - bone: 17-20 GPa
     p['Poisson'] = 0.3       # [-] Poisson's ratio
@@ -134,7 +136,6 @@ def solve(p={}):
     for fix in p['fixations']:
         create_group(fix['nodes'], nods_no, nods_pos, groups, fix['name'])
 
-
     # --------------------------------------------------------------------------
     # print('== ADDING MEDIUM...')
     print('')
@@ -159,6 +160,15 @@ 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]) )
 
+    # 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) )
 

models/dolicorhynchops/dolicorhynchops_10k_extloads.py

+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+# Dolicorhynchops osborni FHSM VP404
+#   10k faces on the mandible surface
+
+
+def parms(d={}):
+    p = {}
+    import os
+    path = os.path.join(os.path.dirname(__file__),'10k')
+    p['bone'] = f'{path}/mandible.stl'
+    p['muscles'] = [
+        {
+            'file': f'{path}/Lmuscle.stl',
+            'force': 100.,
+            # f'{path}/LmuscleF.off',
+            'focalpt': [-100.1458893, -173.13895, 227.3909],
+            'method': 'U'
+        },
+        {
+            'file': f'{path}/Rmuscle.off',
+            'force': 100.,
+            'focalpt': f'{path}/RmuscleF.off',
+            'method': 'T'
+        }
+    ]
+    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'] = [
+        {
+            '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.]
+        }      
+    ]
+
+
+
+    # material properties
+    p['density'] = 1.850e-9  # [T/mm³]
+    p['Young'] = 20000.      # [MPa]
+    p['Poisson'] = 0.3       # [-]
+
+    p.update(d)
+    return p
+
+
+def getMetafor(p={}):
+    import bonemodel as model
+    return model.getMetafor(parms(p))
+
+
+if __name__ == "__main__":
+    import models.bonemodel2 as model
+    model.solve(parms())