add basic extractors to cxxfem (to be continued)

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