From b979bb3ba686270a95d79ffe41862fe2c1d2d140 Mon Sep 17 00:00:00 2001
From: Romain Boman <r.boman@uliege.be>
Date: Wed, 22 May 2024 16:13:10 +0200
Subject: [PATCH] add basic extractors to cxxfem (to be continued)
---
models/bonemodel2.py | 47 ++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 47 insertions(+)
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
index 417ae52..47af45d 100644
--- a/models/bonemodel2.py
+++ b/models/bonemodel2.py
@@ -64,6 +64,7 @@ def parms(d={}):
]
p['loads'] = [] # additional loads
+ p['extractors'] = [] # additional extractors
# material properties
p['Young'] = 17000. # [MPa] elastic modulus - bone: 17-20 GPa
@@ -135,6 +136,10 @@ def solve(p={}):
for fix in p['fixations']:
create_group(fix['nodes'], nods_no, nods_pos, groups, fix['name'])
+ # create groups of additional extractors
+ for extr in p['extractors']:
+ create_group(extr['nodes'], nods_no, nods_pos, groups, extr['name'])
+
# --------------------------------------------------------------------------
# print('== ADDING MEDIUM...')
print('')
@@ -189,6 +194,7 @@ def solve(p={}):
print('extracting results...')
print(f'\tinternal energy = {solver.int_energy:.3f} N.mm')
+ # fixations (reaction forces)
for fix in p['fixations']:
grpname = fix['name']
v = np.array(post.probe('force_vector', grpname))
@@ -198,6 +204,7 @@ def solve(p={}):
print(f'\t{grpname}_Fy = {v[1]:.3f} N')
print(f'\t{grpname}_Fz = {v[2]:.3f} N')
+ # muscle groups
for name, group in mgroups.items():
v=np.array( [0.0, 0.0, 0.0] )
for tag in group.ntags:
@@ -206,6 +213,46 @@ def solve(p={}):
print(f'\t{name}_Fy = {v[1]:.3f} N')
print(f'\t{name}_Fz = {v[2]:.3f} N')
+ # additional extractors
+
+ var2field = {
+ 'dx': ('X', 'mm'),
+ 'dy': ('Y', 'mm'),
+ 'dz': ('Z', 'mm'),
+ 'sigxx': ('smooth_stress_xx', 'MPa'),
+ 'sigyy': ('smooth_stress_yy', 'MPa'),
+ 'sigzz': ('smooth_stress_zz', 'MPa'),
+ 'sigxy': ('smooth_stress_xy', 'MPa'),
+ 'sigxz': ('smooth_stress_xz', 'MPa'),
+ 'sigyz': ('smooth_stress_yz', 'MPa'),
+ 'epsxx': ('strain_xx', ''),
+ 'epsyy': ('strain_yy', ''),
+ 'epszz': ('strain_zz', ''),
+ 'epsxy': ('strain_xy', ''),
+ 'epsxz': ('strain_xz', ''),
+ 'epsyz': ('strain_yz', '')
+ }
+
+ for extr in p['extractors']:
+ # if only one variable is given, convert it to a list
+ if isinstance(extr['variables'], str):
+ extr['variables'] = [extr['variables']]
+
+ for var in extr['variables']:
+ fname = f'{extr["name"]}_{var}.tsv'
+ grpname = extr['name']
+ if var.lower() in var2field:
+ field = var2field[var.lower()][0]
+ units = var2field[var.lower()][1]
+ v = np.array(post.probe(field, grpname))
+ np.savetxt(fname, v, delimiter='\t', header='value')
+ if len(v) == 1:
+ print(f'\t{grpname}_{var} = {v[0]:.3f} {units}')
+ else:
+ print(f'\t{grpname}_{var} = ({len(v)} values in "{fname}")')
+ else:
+ print(f'unknown variable: "{var}"')
+
end_time = time.perf_counter()
elapsed_time = end_time - start_time
print(f'Total Execution time: {elapsed_time:.1f}s')
--
GitLab