diff --git a/cxxfem/src/femSolver.cpp b/cxxfem/src/femSolver.cpp
index 57f33c6e0363a5eb197dc1c84b2992aef9d0ec8d..e2a92ad51be6f047046897d5548e724e7b2858f9 100644
--- a/cxxfem/src/femSolver.cpp
+++ b/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";
}
diff --git a/models/boneload.py b/models/boneload.py
index 103b77324d5cd2cfe839d4c6c89d5cfcf87f733a..35c1387a30d5884fd528d6b275dae0e351e2103a 100644
--- a/models/boneload.py
+++ b/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
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
index 3f2634be0d67b11a173a42597f906144ab59d5d0..3e62d531e80a855a030a699c6b44dad57cc16103 100644
--- a/models/bonemodel2.py
+++ b/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) )
diff --git a/models/dolicorhynchops/dolicorhynchops_10k_extloads.py b/models/dolicorhynchops/dolicorhynchops_10k_extloads.py
new file mode 100644
index 0000000000000000000000000000000000000000..4812b3936c8c239d9be8f9eec61c4a4333f750bf
--- /dev/null
+++ b/models/dolicorhynchops/dolicorhynchops_10k_extloads.py
@@ -0,0 +1,77 @@
+#! /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())