diff --git a/cxxfem/src/femPost.cpp b/cxxfem/src/femPost.cpp
index 1c4119ef30d1fed3832df50c139a937a69ba2940..be8d5e18d5c4a14200dc000fb50baf83a72e9dbd 100644
--- a/cxxfem/src/femPost.cpp
+++ b/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());
}
diff --git a/fossils-win.spec b/fossils-win.spec
index 3862edbbb060c5584c048ce5339e8d21fc1b6303..4884ff97b07782ff844f7d467b5df2549bbae11b 100644
--- a/fossils-win.spec
+++ b/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',
+ )
diff --git a/fossils.iss b/fossils.iss
index 265ca0db683c1371d9dfceba3b4fb43725c2b4ab..d2ed60614b990a1f8348ed3c88c992b8356805da 100644
--- a/fossils.iss
+++ b/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"
diff --git a/models/Clidastes/Clidastes_propython_50k.py b/models/Clidastes/Clidastes_propython_50k.py
index bfe6ccb0679d1d7b1b502bc35dfd68469a236a83..15b001ff350644da4c176beac54f3f2f362d2c64 100644
--- a/models/Clidastes/Clidastes_propython_50k.py
+++ b/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',
diff --git a/models/bonemodel.py b/models/bonemodel.py
index a973806c9ad25e7474eb5602b7af8e5a7c314394..6193e0bfdd2e1c53e3fafbff30f5b8817cbc6174 100644
--- a/models/bonemodel.py
+++ b/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
diff --git a/models/bonemodel2.py b/models/bonemodel2.py
index 417ae5240d77eb4446919856e178f52823e7cc08..25ade80a5748e54b97da8a4ac9757c166d4da3b8 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,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')
diff --git a/models/dolicorhynchops/dolicorhynchops_10k.py b/models/dolicorhynchops/dolicorhynchops_10k.py
index 99e195cd8468bc6c925c8192e455e8ed43f569d5..4ef2f960f0853d48d63c62e0f2fd6f95def3fbbd 100644
--- a/models/dolicorhynchops/dolicorhynchops_10k.py
+++ b/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³]