From 4e62114403088b3d22944d0d52613b4c76ff7f09 Mon Sep 17 00:00:00 2001
From: Romain Boman <r.boman@uliege.be>
Date: Wed, 22 May 2024 12:38:32 +0200
Subject: [PATCH] add additional extractors (Metafor only)

---
 models/Clidastes/50k/Point_5_L.off          |  3 +
 models/Clidastes/Clidastes_propython_50k.py | 23 ++++++
 models/boneload.py                          |  2 +-
 models/bonemodel.py                         | 84 ++++++++++++++++++---
 4 files changed, 102 insertions(+), 10 deletions(-)
 create mode 100644 models/Clidastes/50k/Point_5_L.off

diff --git a/models/Clidastes/50k/Point_5_L.off b/models/Clidastes/50k/Point_5_L.off
new file mode 100644
index 0000000..2e67a62
--- /dev/null
+++ b/models/Clidastes/50k/Point_5_L.off
@@ -0,0 +1,3 @@
+NOFF
+1 0 0
+47.137432 171.505753 9.123382
diff --git a/models/Clidastes/Clidastes_propython_50k.py b/models/Clidastes/Clidastes_propython_50k.py
index 172ea1b..bfe6ccb 100644
--- a/models/Clidastes/Clidastes_propython_50k.py
+++ b/models/Clidastes/Clidastes_propython_50k.py
@@ -100,6 +100,29 @@ def parms(d={}):
             'method': 'T+N'
         }
     ]
+    p['extractors'] = [
+        {
+            'name': 'Point_1_L',
+            'nodes': [9.737800, -107.443642, -41.476234],
+            'variables': [ 'x0', 'y0', 'z0', 'sigvm' ]
+        },
+        {
+            'name': 'Point_2_L',
+            'nodes': [18.720367, -42.742260, -36.978378],
+            'variables': 'sigvm'
+        },
+        {
+            'name': 'Point_3_and_4_L',
+            'nodes': [[33.984665, 40.152878, -24.569260], 
+                      [43.314785, 99.891335, -16.491674]],
+            'variables': 'sigzz'
+        },
+        {
+            'name': 'Point_5_L',
+            'nodes': f'{path}/Point_5_L.off',
+            'variables': [ 'sigvm', 'dx', 'dy', 'dz' ]
+        }
+    ]
 
     # material properties
     p['density'] = 1.662e-9  # [T\mm³]
diff --git a/models/boneload.py b/models/boneload.py
index 50f96a1..e4fb3ae 100644
--- a/models/boneload.py
+++ b/models/boneload.py
@@ -22,7 +22,7 @@ colors = vtk.vtkNamedColors()
 
 
 def load_msh(meshfile):
-    """loads a mesh in .off or .ply format
+    """loads a mesh in .off, .stl or .ply format
     returns: 
         a list of node coordinates 
         a list of triangles defined by 3 node tags
diff --git a/models/bonemodel.py b/models/bonemodel.py
index 00d1c7a..a973806 100644
--- a/models/bonemodel.py
+++ b/models/bonemodel.py
@@ -8,11 +8,6 @@ from wrap import *
 import boneload
 import os
 
-# enable full parallelism
-StrVectorBase.useTBB()
-StrMatrixBase.useTBB()
-ContactInteraction.useTBB()
-
 def parms(d={}):
     """default parameters of the model
     """
@@ -64,6 +59,7 @@ def parms(d={}):
     ]
 
     p['loads'] = [] # additional loads
+    p['extractors'] = [] # additional extractors
 
     p['food'] = 'fixteeth'  # type of food: 'fixteeth' or 'cylinder'
     
@@ -146,6 +142,10 @@ def getMetafor(p={}):
     for fix in p['fixations']:
         create_group(fix['nodes'], nods_no, nods_pos, domain, groups, fix['name'])
 
+    # create groups of additional extractors
+    for extr in p['extractors']:
+        create_group(extr['nodes'], nods_no, nods_pos, domain, groups, extr['name'])
+
     # --------------------------------------------------------------------------
 
     # material 
@@ -242,12 +242,22 @@ def getMetafor(p={}):
     solver = DSSolver()
     solman.setSolver(solver)
 
+    # enable full parallelism
+    StrVectorBase.useTBB()
+    StrMatrixBase.useTBB()
+    ContactInteraction.useTBB()
+
     # history curves
+    # by default, extract all forces on all groups
     valuesmanager = metafor.getValuesManager()
     idx=1
     valuesmanager.add(idx, MiscValueExtractor(metafor,EXT_T),'time'); idx+=1
     for name, group in groups.items():
+        # skip the main mesh and its surface
         if name in ['all','surf']: continue
+        # skip all additional extractors
+        if name in [e['name'] for e in p['extractors']]: continue
+        # muscles and loads should remain
         for key, c in zip([TX,TY,TZ],['x','y','z']):
             fname = f'{name}_F{c}'
             valuesmanager.add(idx, DbNodalValueExtractor(group, Field1D(key,GF1)), SumOperator(), fname)
@@ -259,11 +269,63 @@ def getMetafor(p={}):
             valuesmanager.add(idx, InteractionValueExtractor(ci, key), fname)
             metafor.getTestSuiteChecker().checkExtractor(idx)
             idx+=1
-    
+
+    # add additional extractors
+
+    # conversion dicts
+    var2field_db = {
+        'dx': Field1D(TX, RE),
+        'dy': Field1D(TY, RE),
+        'dz': Field1D(TZ, RE),
+        'x0': Field1D(TX, AB),
+        'y0': Field1D(TY, AB),
+        'z0': Field1D(TZ, AB),
+        'fx': Field1D(TX, GF1),
+        'fy': Field1D(TY, GF1),
+        'fz': Field1D(TZ, GF1)
+    }
+
+    var2field_if = { 
+        'sigvm': IF_EVMS, 
+        'sigxx': IF_SIG_XX, 
+        'sigyy': IF_SIG_YY, 
+        'sigzz': IF_SIG_ZZ,
+        'sigxy': IF_SIG_XY, 
+        'sigyz': IF_SIG_YZ, 
+        'sigxz': IF_SIG_XZ,
+        'epsxx': IF_GL_STRAIN_XX,
+        'epsyy': IF_GL_STRAIN_YY,
+        'epszz': IF_GL_STRAIN_ZZ,
+        'epsxy': IF_GL_STRAIN_XY,
+        'epsyz': IF_GL_STRAIN_YZ,
+        'epsxz': IF_GL_STRAIN_XZ 
+    }
+
+    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}'
+            target = groups[extr['name']]
+            # nodal values
+            if var.lower() in var2field_db:
+                field = var2field_db[var.lower()]
+                valuesmanager.add(idx, DbNodalValueExtractor(target, field), fname)
+            elif var.lower() in var2field_if:
+                field = var2field_if[var.lower()]
+                valuesmanager.add(idx, IFNodalValueExtractor(target, field), fname)
+            else:
+                raise Exception(f'Unknown variable: {var}')
+            if target.getNumberOfMeshPoints() == 1:
+                metafor.getTestSuiteChecker().checkExtractor(idx)
+            idx+=1
+
     return metafor
 
 # ------------------------------------------------------------------------------
-# utility functions involing Metafor objects
+# utility functions involving Metafor objects
 
 
 def import_mesh(domain, filename, use_gmshOld=True):
@@ -342,7 +404,7 @@ def create_group(mesh_file_or_coord, all_no, all_coords, domain, groups, gname=N
     Then, the routine creates a group named "gname" (built from mesh_file_or_coord if None)
 
     mesh_file_or_coord can be either:
-      a filename
+      a filename (.off, .ply, .stl)
       an array of points: [ [x1,y1,z1], [x2,y2,z2], ... ]
       one single point: [x1,y1,z1]
     """  
@@ -373,6 +435,10 @@ def create_group(mesh_file_or_coord, all_no, all_coords, domain, groups, gname=N
     if gname is None:
         gname, _ = os.path.splitext(os.path.basename(mesh_file_or_coord))
     group = groupset.add(Group(groupset.getLargestNo()+1))
+
+    if gname in groups:
+        raise Exception(f'create_group: group "{gname}" already exists!')
+
     groups[gname] = group
 
     # add all the identified vertices to the group
@@ -397,7 +463,7 @@ class MuscleGroup:
 
 
 def ply2mtf(domain, meshfile):
-    """load a surface into Metafor's memory from a ply file
+    """load a surface into Metafor's memory from a mesh file (ply, off or stl)
     """
     fullpath = os.path.join(os.path.dirname(__file__), meshfile)
     nodes, tris = boneload.load_msh(fullpath)
-- 
GitLab