Merge branch 'extractors' into 'master'

Extractors

See merge request !2

cxxfem/src/femPost.cpp

@@ -24,11 +24,45 @@ Post::~Post()
         std::cout << "~Post()\n";
 }
 
-/// rebuild solution vectors for gmsh and add views
+/// rebuild solution vectors for gmsh and add views:
+///
+///  X
+///  Y
+///  Z
+///  displacement_vector
+///  force_vector
+///  stress_tensor
+///      stress_xx
+///      stress_yy
+///      stress_zz
+///      stress_xy
+///      stress_xz
+///      stress_yz
+///  strain_tensor
+///      strain_xx
+///      strain_yy
+///      strain_zz
+///      strain_xy
+///      strain_xz
+///      strain_yz
+///  smooth_stress_tensor
+///      smooth_stress_xx
+///      smooth_stress_yy
+///      smooth_stress_zz
+///      smooth_stress_xy
+///      smooth_stress_xz
+///      smooth_stress_yz
+///  smooth_strain_tensor
+///      smooth_strain_xx
+///      smooth_strain_yy
+///      smooth_strain_zz
+///      smooth_strain_xy
+///      smooth_strain_xz
+///      smooth_strain_yz
 
 void Post::build_views()
 {
-    std::cout << "creating views...\n";
+    std::cout << "creating views..." << std::endl;
     // nodal dofs
     add_dofs_view(solver.pbl, solver.dofs, solver.X, solver.dofPart, solver.dofIdx);
 
@@ -330,6 +364,19 @@ Post::probe(std::string const &name, std::string const &grpname)
     // get physical group from grpname
     Group *grp = solver.pbl.groups_by_name.at(grpname); // TODO: better handling of errors
 
+    std::vector<std::size_t> nodeTags;
+    std::vector<double> coord;
+    gmsh::model::mesh::getNodesForPhysicalGroup(grp->dim, grp->tag, nodeTags, coord);
+    // j'ai essayé de boucler sur les entités du groupe et faire un getNode de l'entité 
+    // mais ça ne me retourne rien (?)
+    // std::cout << "nodeTags.size()=" << nodeTags.size() << '\n';
+
+    // retourne les coordonnées des noeuds du groupe si demandé
+    if(name=="coords")
+    {   
+        return coord;
+    }
+
     // get the view
     auto it = views.find(name);
     if (it == views.end())
@@ -339,14 +386,6 @@ Post::probe(std::string const &name, std::string const &grpname)
     // probe values
     std::vector<double> all_values;
 
-    std::vector<std::size_t> nodeTags;
-    std::vector<double> coord;
-    gmsh::model::mesh::getNodesForPhysicalGroup(grp->dim, grp->tag, nodeTags, coord);
-    // j'ai essayé de boucler sur les entités du groupe et faire un getNode de l'entité 
-    // mais ça ne me retourne rien (?)
-
-    // std::cout << "nodeTags.size()=" << nodeTags.size() << '\n';
-
     for (size_t i = 0; i < nodeTags.size(); ++i)
     {
         std::vector<double> values;
@@ -365,6 +404,21 @@ Post::probe(std::string const &name, std::string const &grpname)
         // std::cout << "\tgot from gmsh: " << values << '\n';
         // std::cout << "\tdistance = " << distance << std::endl;
 
+        if(values.size() == 9) // tensor -> compute j2
+        {
+            // xx, xy, xz, yx, yy, yz, zx, zy, zz
+            //  0,  1,  2,  3,  4,  5,  6,  7,  8
+            double xx = values[0];
+            double yy = values[4];
+            double zz = values[8];
+            double xy = values[1];
+            double zx = values[6];
+            double yz = values[5];
+            double j2 = sqrt( ((xx-yy)*(xx-yy) + (yy-zz)*(yy-zz) + (zz-xx)*(zz-xx) )/2 + 3*(xy*xy+yz*yz+zx*zx) );
+            values.clear();
+            values.push_back(j2);
+        }
+
         all_values.insert(all_values.end(), values.begin(), values.end());
     }
 

fossils-win.spec

@@ -2,7 +2,7 @@
 
 # BUILD DIST:
 #   load environment
-#   pyinstaller --noconfirm fossils.spec
+#   pyinstaller --noconfirm fossils-win.spec
 
 block_cipher = None
 
@@ -63,13 +63,29 @@ exe = EXE(
     entitlements_file=None,
     icon='gui/icons/fossils.ico',
 )
-coll = COLLECT(
-    exe,
-    a.binaries,
-    a.zipfiles,
-    a.datas,
-    strip=False,
-    upx=False, # works well but slower runtime
-    upx_exclude=[],
-    name='fossils',
-)
+
+# when using VBox + mesa drivers installed on c:\Python310\python.exe, 
+# mesa dlls are copied into the installer which makes
+# gmsh display crash at runtime.
+# => we manually remove the libraries related to openGL (they are not present when built outside a vbox)
+# see https://stackoverflow.com/questions/58097858/how-to-exclude-opengl32sw-dll-from-pyqt5-library-when-using-pyinstaller
+to_keep = []
+to_exclude = {'opengl32.dll', 'libglapi.dll', 'libgallium_wgl.dll'}
+for (dest, source, kind) in a.binaries:
+    if os.path.split(dest)[1].lower() in to_exclude:
+        print(f'removing {source}')
+        continue
+    to_keep.append((dest, source, kind))
+a.binaries = to_keep
+
+if 1:
+    coll = COLLECT(
+        exe,
+        a.binaries,
+        a.zipfiles,
+        a.datas,
+        strip=False,
+        upx=False, # works well but slower runtime
+        upx_exclude=[],
+        name='fossils',
+    )

fossils.iss

@@ -2,7 +2,7 @@
 ; SEE THE DOCUMENTATION FOR DETAILS ON CREATING INNO SETUP SCRIPT FILES!
 
 #define MyAppName "Fossils"
-#define MyAppVersion "1.3"
+#define MyAppVersion "1.4"
 #define MyAppPublisher "University of Liege"
 #define MyAppURL "https://gitlab.uliege.be/rboman/fossils"
 #define MyAppExeName "fossils.exe"

models/Clidastes/Clidastes_propython_50k.py

@@ -7,7 +7,7 @@ def parms(d={}):
     p = {}
     import os
     path = os.path.join(os.path.dirname(__file__),'50k')
-    p['bone'] = f'{path}\Clidastes_50k.ply'
+    p['bone'] = f'{path}/Clidastes_50k.ply'
     p['fixations'] = [
         {
             'name': 'contact_pts',

models/bonemodel.py

@@ -318,6 +318,8 @@ def getMetafor(p={}):
                 valuesmanager.add(idx, IFNodalValueExtractor(target, field), fname)
             else:
                 raise Exception(f'Unknown variable: {var}')
+            # if the target is a single node, display the extractor [TSC-EXT]
+            # otherwise, the user should read the ascii file in the workspace.
             if target.getNumberOfMeshPoints() == 1:
                 metafor.getTestSuiteChecker().checkExtractor(idx)
             idx+=1

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,65 @@ 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': ('smooth_strain_xx', ''), 
+        'epsyy': ('smooth_strain_yy', ''), 
+        'epszz': ('smooth_strain_zz', ''), 
+        'epsxy': ('smooth_strain_xy', ''), 
+        'epsxz': ('smooth_strain_xz', ''), 
+        'epsyz': ('smooth_strain_yz', ''),
+        'sigvm': ('smooth_stress_tensor', 'MPa')  # post.probe computes sigvm                  
+    }
+
+    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']
+            v = []
+            units = ""
+            if var.lower() in var2field:
+                # extract values from gmsh views
+                field = var2field[var.lower()][0]
+                units = var2field[var.lower()][1]
+                v = np.array(post.probe(field, grpname))              
+            elif var.lower() in ['x0', 'y0', 'z0']:
+                # extract coordinates of the nodes
+                units = 'mm'
+                v = np.array(post.probe('coords', grpname))
+                v = np.reshape(v, (v.size//3, 3))
+                v = v[:, ['x0', 'y0', 'z0'].index(var.lower())] 
+            elif var.lower() in ['fx', 'fy', 'fz']: 
+                # extract forces
+                units = 'N'
+                v = np.array(post.probe('force_vector', grpname))
+                v = np.reshape(v, (v.size//3, 3))
+                v = v[:, ['fx', 'fy', 'fz'].index(var.lower())]
+            else:
+                print(f'unknown variable: "{var}"')
+
+            if len(v) > 0:
+                np.savetxt(fname, v, delimiter='\t', header=var.lower())
+                if len(v) == 1:
+                    print(f'\t{grpname}_{var} = {v[0]:.3f} {units}')
+                else:
+                    print(f'\t{grpname}_{var} = ({len(v)} values in "{fname}")')                
+
+
     end_time = time.perf_counter()
     elapsed_time = end_time - start_time
     print(f'Total Execution time: {elapsed_time:.1f}s')

models/dolicorhynchops/dolicorhynchops_10k.py

@@ -42,6 +42,23 @@ def parms(d={}):
             'direction': ['x', 'z']
         }
     ]
+    p['extractors'] = [
+        {
+            'name': 'tooth1',
+            'nodes': [-10.20362, -17.46838, -229.9061],
+            'variables': [ 'sigvm' ]
+        },
+        {
+            'name': 'tooth2',
+            'nodes': [-11.92466, 26.3042, -229.5354],
+            'variables': [ 'sigxx', 'sigyy', 'sigzz', 'sigvm' ]
+        },
+        {
+            'name': 'Lmuscle',
+            'nodes': f'{path}/Lmuscle.stl',
+            'variables': [ 'x0', 'y0', 'z0', 'sigvm' ]
+        }
+    ]
 
     # material properties
     p['density'] = 1.850e-9  # [T/mm³]