mlp6/fem

Finite Element Modeling (FEM) Code: Python Tools, Field II Intensity Field Solution, LS-DYNA Pre/Post Processing

Acoustic Radiation Force FEM Tools

Finite Element Method (FEM) tools to simulate acoustic radiation force excitations and resultant shear wave propagation.

All software in this repository is licensed under the Apache v2.0 license, as detailed in the LICENSE file.

If you are using the FEM simulation tools in work that you publish, then please consider citing the following manuscript:

Palmeri ML, Sharma AC, Bouchard RR, Nightingale RW, Nightingale KR. "A Finite-Element Method Model of Soft Tissue Response to Impulsive Acoustic Radiation Force," IEEE UFFC, 52(10): 1699-1712, 2005

Also please cite the following manuscript if you use Field II:

J.A. Jensen and N. B. Svendsen: Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers, IEEE Trans. Ultrason., Ferroelec., Freq. Contr., 39, pp. 262-267, 1992.

Installation

These tools have been developed in a Linux environment, which is the recommended platform for usage. Most of the Python code, and some of the C code, will work on Windows and Mac platforms, but build tools for these systems is not outlined here and robust testing does not exist for those environments (i.e., use at your own risk). The Python tools require a version >= 3.8.

There are three different methods to use/install this package:

1. [RECOMMENDED] Install this with pip in a local virtualenv: pip install --prefer-binary -e git+https://github.com/mlp6/fem.git#egg=fem. If swig is not available on your machine, install with pip install --prefer-binary -e git+https://github.com/mlp6/fem.git@no-swig#egg=fem.
2. Clone this source directory and manually install it in your local virtualenv: pip install -e .
3. Clone this source directory and work directly with it by defining PYTHONPATH to include the root directory this repository was cloned into. This is the easiest method if you are planning to actively edit/develop the codebase. You can do this on-the-fly for a given interactive ipython session with syntax like: PYTHONPATH=$HOME/fem ipython.

Documentation

The latest documentation is automatically generated from package docstrings using Sphinx and can be built in docs/. That directory also contains static HTML files that can be included in the documentation, as defined in docs/index.rst. The documentation synced to the root docs/ directory level are also rendered at https://mlp6.github.io/fem/.

To build the documentation:

1. Make sure that PYTHONPATH includes the fem package.
2. Make sure that the virtualenv in activated, which includes the sphinx-build package.
3. Within the docs/ directory: make html
4. Run docs/rsync_build.sh to bring the newly-built HTML files and associated source files into docs/.

Issues

Please file any bug reports, features requests, etc. using the GitHub Issues.

Contributors

• Mark Palmeri ([email protected])
• Joey Richardson (@jrichardson97)
• Anna Phillips (@annap2277)
• Ningrui Li ([email protected])
• Mallory Selzo (UNC-CH)
• Chris Moore ([email protected])
• David Bradway ([email protected])
• Nick Bottenus ([email protected])
• Brian Bigler ([email protected])
• Carl Herickhoff ([email protected])
• Sam Lipman ([email protected])
• Anna Knight ([email protected])
• Matthew Huber (@matthew-huber)
• Felix Jin (@fqjin)

fem's Issues

update examples to reflect v7.x releases

field2dyna and dynaField requires Field_II_Pro

been a while since i wrote a git issue so apologies if this isn't the correct formatting

1. even though there is a nargin catch in field2dyna it still errors if it does not have the correct number of inputs (there are catches for nargin 7,8, and 11 but 9 and 10 are missing). I think i can fix this by changing the nargin catch for threads to 8 not 7, and for lownslow to 9 not 8. (assuming those without pro will not be running multithreaded or messing with lownslow.
2. in dynaField two lines cause matlab to crash if you do not have field_ii_pro (which students in my class do not):
   set_field('threads', threads);
   disp(sprintf('PARALLEL THREADS: %d', threads));
   If you comment these out it runs fine. but non-pro field doesn't have set_field i think

let me know if you want me to submit my version with these changes? n

add matlab function to read JSON probe files

update thermal documentation

https://github.com/mlp6/fem/blob/master/docs/thermal.rst

field/formatExpImpResp: hard-coded time indices

Time indices for Start/Stop Index are hard-coded to +/- 400 ns, causing errors for very low frequencies probes. Need to more elegantly code these limits.

post-c: reduce double -> single

we don't need the full dynamic range of double-precision numbers, which are eating up a lot of extra space; will convert to single-precision floats with truncates range

rm legacy eval() approach to defining probes

specify subset of timesteps to extract

create_res_sim currently extracts all timesteps... would be nice to limit this if a user just needs some or a sparse set for memory concerns, etc. might be a feature to enable in combinations with issue #83

dig up old thermal sims from 2003 TUFFC manuscript

probes.m files need focus_m

fem/field/field2dyna.m

Line 51 in 390797f

FIELD_PARAMS.focus = focus;

Some files (vf105) in the ultrasound/probes repository use FIELD_PARAMS.focus_m instead of FIELD_PARAMS.focus. Change one to be in concordance with the other?

create_res_sim_mat: LONG precision index error

  File "/dscrhome/mlp6/pkgs/fem/post/create_res_sim_mat.py", line 120, in extract_arfi_data
    arfidata[j, i, t - 1] = disp_scale * zdisp[nodeid]
IndexError: index 25897041 is out of bounds for axis 0 with size 25897041

Seems to happen all the way back to versions 6.5.x. This error is happening only on a HUGE mesh, so not sure what is going on.

add ReadTheDocs badge to README

Parse tprint Output File

Including the *DATABASE_TPRINT card in the .dyn file that controls the simulation results in temperatures at each node being output for the timestep specified. A program should be written to parse this output and export it into a .mat file. The program should have the option of selecting the entire node set temperature data for exporting, or just exporting the data associated with a desired face of the simulated volume.

post-c: check for complete last time step

last time step can sometimes be incompletely written; either need to check node count and use only if full, or just ignore the last timestep

fieldC: high memory usage

Memory usage in my tests is high, even when running single-threaded, which usually is not a memory hog. This memory usage also seems to grow, so this isn't something static that is preallocated.

Need to profile this is more detail.

integrate Field II Pro C

utilize the C version of Field II Pro to run on systems w/o Matlab licenses (but don't drop Matlab support)

examples/vf105/run.sh (bc.py)

The CLI syntax for this example needs to be update to reflect the explicit choice of PML (--pml) or NONREFLECT (--nonreflect) since the latest version does not have a default.

pointloads.py: handle direction of loads in show_image_plane()

The method currently assumes there is a single magnitude at a single, unique node definition in the point loads file; that is not true when there are multiple load components at a node. I think this will currently overwrite an image plane location with whatever the last component magnitude is in the point loads file.

SWIG build (create_disp_dat.c)

Need to find a way to automatically build the SWIG module... maybe setup.py?

fieldC compilation error

Getting the following compilation error:

$ make
gcc -I . -I /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library checkOnAxis.o checkUniform.o correctAxialLens.o dynaField.o dynaWrite.o field2dyna.o gaussPulse.o readMpn.o main.c cJSON.c /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a -lpthread -lm -O
main.c:8:1: warning: return type defaults to ‘int’ [-Wimplicit-int]
 main()
 ^~~~
/usr/bin/x86_64-linux-gnu-ld: field2dyna.o: relocation R_X86_64_32 against `.rodata' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: gaussPulse.o: relocation R_X86_64_32 against `.rodata' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: readMpn.o: relocation R_X86_64_32 against `.rodata' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(field_II_main.o): relocation R_X86_64_32 against `.rodata.str1.1' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(transducers.o): relocation R_X86_64_32 against `.rodata.str1.8' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(calculations_par.o): relocation R_X86_64_32 against symbol `calc_scat_call' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(general.o): relocation R_X86_64_32 against `.rodata.str1.1' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(thread_lib.o): relocation R_X86_64_32 against `.rodata.str1.8' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(calculations.o): relocation R_X86_64_32 against `.rodata.str1.8' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(triangles.o): relocation R_X86_64_32 against `.rodata.str1.8' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: /home/mlp6/Documents/MATLAB/Field_II_Pro/c_library/lib_Field_II.a(lines.o): relocation R_X86_64_32 against `.rodata.str1.1' can not be used when making a PIE object; recompile with -fPIC
/usr/bin/x86_64-linux-gnu-ld: final link failed: Nonrepresentable section on output
collect2: error: ld returned 1 exit status
Makefile:5: recipe for target 'main' failed
make: *** [main] Error 1

@ndanieley did you hit any errors like this?

associated with PR #51

allow test_savepvd to test actual saved files

original attempt at a test failed b/c the py.test.tmdir fixture created a dynamic path written to the PVD file, so the filecmp.cmp() operation would fail with the validation data. Need to change the code to not write the entire path to the PSD file, but instead simply the basename of the VTR file.

create_pointloads_vtk.py: not running (numElem)

Extracting data . . .
Writing node positions
Traceback (most recent call last):
  File "/home/mlp6/projects/fem/mesh/create_pointloads_vtk.py", line 430, in <module>
    main()
  File "/home/mlp6/projects/fem/mesh/create_pointloads_vtk.py", line 56, in main
    create_vtu(args)
  File "/home/mlp6/projects/fem/mesh/create_pointloads_vtk.py", line 136, in create_vtu
    numNodes, numElems = writeNodePositions(loadout, args, 'vtu')
  File "/home/mlp6/projects/fem/mesh/create_pointloads_vtk.py", line 237, in writeNodePositions
    if args.numElem[0] is None:
AttributeError: 'Namespace' object has no attribute 'numElem'

write Field II mesh interpolation function

instead of having to always re-run Field II when using a new mesh, we should have a function that will take an existing 3D Field II intensity output and 3D map it onto a new grid

get sphinx docs to render in GitHub pages (docs/)

what does `xdc_get()` return (C version)

@ndanieley isn't sure what this should be returning and why it appears to be returning nothing (see the fieldC branch code review in progress)

[calc_disp_dat_c] progressive status

Generating the disp.dat file from nodout can take a while, and it would be nice to see progressive output status to know how many timesteps have been completed.

convert rst -> Google-style docstring

have field2dyna.c read input parameters from a JSON file

create_pointloads_vtk.py: clips filenames with '.'

selectively import SWIG create_disp_dat

Utilize a try/except structure to attempt to import the SWIG create_disp_dat, and default to python implementation if it fails.

break out `makeLoadsTemps.m` into thermal-specific python driver

legacy mess of code; use common module that can then write loads or temps

add compilation of C code to setup process

res_sim.py: axes deprecation

MatplotlibDeprecationWarning: Adding an axes using the same argument
s as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this w
arning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.

create_res_sim_mat: extract >=1 displacement component

currently this method is called multiple times when different components of displacement have to be extracted; would be nice to just pass a vector of the components to extract

update documentation to use JSON probe defs

read timesteps in parallel

Attempt to accelerate the post-processing by reading each timestep from nodout in parallel instead of serially.

swig make error

make: *** No rule to make target create_disp_dat_wrap.c, needed by _create_disp_dat.so. Stop.

create VTK from binary dat file

create_disp_dat.py can create the binary data file and the VTK-formatted displacement file from a raw nodout file, but disp.dat files cannot be directly converted to VTK format.

fieldC: field2dyna compilation

formatExpImpResp.c:5:28: fatal error: gsl/gsl_interp.h: No such file or directory
 #include <gsl/gsl_interp.h>

This error was on durmstrang.

dt extraction limited to `DATABASE_NODOUT`

The current extract_dt() function in create_res_sim only looks for the DATABASE_NODOUT keywork field to find the time step increment, but some models only have DATABASE_BINARY_D3PLOT defined. Would be good to have this be more robust. Either look for both and choose one by default, have the user specify the string to look for, or maybe allow the user to just specify it. ??

field2dyna: lens correction

Sam discovered that the axial lens correction is also incorrectly being applied in the lateral and elevation dimensions, creating an incorrect offset for the solved fields.

Integrate Transducer Face Heating with Thermal Sims

Branch mlp6-thermal has two files that have hard-coded adaptations for face heating simulations: field2dynaFaceHeating.m and makeFaceHeating.m. These could either be incorporated into the existing field2dyna.m and makeLoadsTemps.m or simplified to operate as standalone components in a face heating simulation workflow for later combination with the other thermal simulations
The hard-coded transducer dimensions should be revised to read from the probe configuration file.

add unit tests

add `xdc_center_focus()` to field2dyna

Need to add the ability to offset an excitation from the "center" (x = 0 and y = 0) of a probe.

impulseResponse parameters

from @ndanieley

"impulseResponse":
{
        "f0": 4000000.0,
        "phase": 0.0,
        "bw": 55.0,
        "wavetype": "gaussian"
},

I don't use any of those variables in my code, and as far as I can tell, you
don't use them in the matlab code. am I missing something?

Tukey example syntax error

The two variables being passed into format are incorrectly put into a tuple, which creates an IndexError:

fem/examples/tukey_gauss_qsym_pml/run.py

Line 36 in ebe9e1e

system('ls-dyna-d ncpu={} i={}'.format((NTASKS, DYNADECK)))

These inputs should not be in a tuple!

create_res_sim_mat: cannot find ele = 0 image plane

In running a mesh that is half-symmetry, and not centered about lat = 0, the following error occurred when trying to create the res_sim.mat file from disp.dat:

In [3]: create_res_sim_mat(DYNADECK)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-3-a01d3b524844> in <module>()
----> 1 create_res_sim_mat(DYNADECK)

~/fem/post/create_res_sim_mat.py in run(dynadeck, disp_comp, disp_scale, ressim, nodedyn, dispout, legacynodes)
     40     [snic, axes] = fem_mesh.SortNodeIDs(node_id_coords)
     41
---> 42     image_plane = extract_image_plane(snic, axes, ele_pos=0.0)
     43
     44     header = read_header(dispout)

~/fem/post/create_res_sim_mat.py in extract_image_plane(snic, axes, ele_pos)
    196
    197     ele0 = np.min(np.where(axes[0] >= ele_pos))
--> 198     image_plane = np.squeeze(snic[ele0, :, :]).astype(int)
    199
    200     return image_plane

ValueError: Can't cast from structure to non-structure, except if the structure only has a single field.

have setup.py pull down examples

examples/ are not being pulled down since they are not importable parts of the package, but then the example files are not available locally for someone to easily copy and use. Would be worth pulling these down for a normal end user to have accessible.

linear.m

getting a strange class error when trying to use linear.m in place of proprietary probe definitions

scaleFactor not used in disp.dat scaling

scale_disp_dat.cpp doesn't actually use the scaleFactor to scale the displacement data (uses a hard-coded value of 0.1)