bugpyFormex - Bugs: bug #45061, Faster computation of element...


bug #45061: Faster computation of element volumes for hex8 in Mesh

Submitter:  francesco <francio>
Submitted:  Sat 09 May 2015 08:58:02 PM UTC
Category:  Core Severity:  3 - Normal
Item Group:  Feature request Status:  None
Privacy:  Public Assigned to:  bverheg
Open/Closed:  Open
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Wed 20 May 2015 06:32:22 AM UTC, comment #15: 

Very nice. this was the sign I was expecting.
Anyway for next time you can just give directions so that we can implement things the way you like, speeding up the improvements and needed changes.

At this point I suppose the C or Fortran implementations are not needed anymore.
So now the function can be included in geomtools.levelVolumes,
be called in Mesh.levelVolumes or just be called directly?

francesco <francio>
Group Member
Tue 19 May 2015 09:38:04 PM UTC, comment #14: 

In 7e4122a a new function geomtools.hexVolume was added computing the volume of hexahedral elements using the formulas mentioned in https://savannah.nongnu.org/bugs/?45061#comment12
This results in a tremendous speedup of the volume computation in large models (probably the conversion takes a lot of time?):

# hex:        1; old: 0.015084s; new: 0.000328s; maxdiff: 6.5193e-09
# hex:       10; old: 0.004398s; new: 0.000333s; maxdiff: 7.4506e-09
# hex:      100; old: 0.007259s; new: 0.000394s; maxdiff: 1.8626e-08
# hex:     1000; old: 0.031803s; new: 0.000478s; maxdiff: 2.9802e-08
# hex:    10000; old: 0.224211s; new: 0.001738s; maxdiff: 3.7253e-08
# hex:   100000; old: 2.401208s; new: 0.016304s; maxdiff: 4.4703e-08
# hex:  1000000; old: 25.634602s; new: 0.175133s; maxdiff: 5.9605e-08

The new function should be used to compute the volume for all hex type elements. Maybe more accurate formulas can be derived for hex20/hex27.

Benedict Verhegghe <bverheg>
Group administrator
Tue 19 May 2015 08:16:20 PM UTC, comment #13: 

In b1a81c7 arraytools.vectorTripleProduct was changed to use a determinant formula, resulting in an important speedup on large datasets. Here are results for 1 to 10 million tetrahedrons.

# tet:        1; old: 0.000174s; new: 0.000109s; maxdiff: 0.0000e+00
# tet:       10; old: 0.000145s; new: 0.000109s; maxdiff: 0.0000e+00
# tet:      100; old: 0.000173s; new: 0.000097s; maxdiff: 0.0000e+00
# tet:     1000; old: 0.000529s; new: 0.000142s; maxdiff: 0.0000e+00
# tet:    10000; old: 0.003892s; new: 0.000369s; maxdiff: 0.0000e+00
# tet:   100000; old: 0.031370s; new: 0.002073s; maxdiff: 0.0000e+00
# tet:  1000000; old: 0.252968s; new: 0.030414s; maxdiff: 0.0000e+00
# tet: 10000000; old: 2.574564s; new: 0.379788s; maxdiff: 0.0000e+00

Benedict Verhegghe <bverheg>
Group administrator
Tue 19 May 2015 04:18:24 PM UTC, comment #12: 

What sign do you expect?

I first need (a lot) more information on what you want/need.
- individual result for every element?
- how accurate?
- how much faster?

Then we need to profile the current code and see where the most time is spent.

Possible ways for speedup:
- Using formulas like in http://www.osti.gov/scitech/servlets/purl/632793
- Implement these (or the current) algorithms in C
- Use numerical integration (with low or high accuracy). Implement these in Fortran (as I have them already, just need to adjust and construct interface).

Benedict Verhegghe <bverheg>
Group administrator
Tue 19 May 2015 06:40:14 AM UTC, comment #11: 

In commit 36adca0 the function meshQuality has been added to compute all vtk quality measures.

Now if we like we can call this function for some specific measures (like vtkVolumes ) as Gianluca suggested , but I think we should not add many.

For the volumes I still think we should find a way not to use vtk.

Benedict can you give a sign on this?

francesco <francio>
Group Member
Wed 13 May 2015 01:34:52 PM UTC, comment #10: 

Gianluca regarding to teh vtk*Quality

I would use the MeshQuality instead of the CellQuality.
the Mesh quality return derived measures as well which can be useful.

we can just use quality parameter as a string and use the eval function with no checks ( up to the user to select the correct one) to keep the function very short, something like this pseudo code

def meshQuality(mesh,quality):
   convert to vtkobj
   elname = get the elemnt vtk name (Hex,Triangle ect )
   q = vtkMeshQuality(vtkobj)

   eval('q.Set' + elname + 'QualityMeasureTo' + quality +'()')
   get the quality array
and we avoid to add as many functions as the quality measures.

francesco <francio>
Group Member
Tue 12 May 2015 11:23:34 AM UTC, comment #9: 

I would not do like you said

I have already implemented the mesh quality functions based on the jacobian and the norms. I will put them as soon as possible after I finish the tests and I write the paper on the mesh quality improvemnet based on that.

this vtk things would be just a temporary solution for some stuff not included in pyformex, a general evaluateMesh will be more compact and we will not add so many different functions. this way teh user ask directly what he needs otherwise eveytime we need something from the  _vtkCellQuality we will keep on adding specific functions that call it. beside I think it is a
easy way just to use a evaluateMesh(mesh, quality). and we keep it as you said with no converisons as teh user can do it outside the function.

For teh volumes in vtk I remember I read something once but I might be wrong or maybe it is not a problem anymore

francesco <francio>
Group Member
Tue 12 May 2015 10:40:23 AM UTC, comment #8: 

Hi Francesco,

I would propose the following:

- create a function in vtk_itf.py to use the vtkCellQuality with tri/quad/tet/hex and with all types of quality (including volume).
I suggest the name "_vtkCellQuality", with underscore
In this function I would not introduce conversions, but leave as it is in vtk, because conversion are planar approximations.

- create more specific functions to check a specific metric
I suggest the name "vtkVolumes", "vtkScaledJacobians"

Btw, why do you think that vtk volume is wrong? If you sum all the vtk volumes you get the correct total volume. I think that vtk uses a different way to compute the volumes which is not necessarily wrong, but it is different in case of warped quads. Splitting quads into triangles is only an approximation. The read quad surface is probably the ruled surface, which is curved.

gianluca de santis <gianlucadesi>
Group Member
Tue 12 May 2015 10:16:39 AM UTC, comment #7: 

I think with our own c++ functions with can have faster computation of volumes,

The correct results is the one in pyformex but
with the conversion to tet4 and then using the vtkTet* instead of the vtkHex* methods the results are the same. I did not cjheck for the quad

I think we can introduce this function in the vtk plugin for all the additional information still lacking in pyFormex.
Maybe we can introduce a function that wraps the vtk one like

evaluateQuality(mesh, targetQuality):

and make some due conversions of quad and hex first in pyformex to and tri3-x and  tet4

Still an improvement of the pyformex functions is needed.
Gianluca will you do the additin to vtk?

francesco <francio>
Group Member
Tue 12 May 2015 09:42:25 AM UTC, comment #6: 


Very interesting topic. We really need some quality metrics for (large) meshes. It seems that a lot of metrics are ready for us in VTK. Please run the attached script (volume of hex elements is 3x faster with vtk). I think we should create one general pyFormex interface function to vtk.vtkMeshQuality()

Regarding hex volumes:

Volume is indeed strange for hex elements because of the non-planarity: vtk and pyformex return different volumes in case of non-planar faces (including negative volume elements).

  • Minor issue in the current pyFormex implementation:

hex.volume() and sum(hex.volumes()) return different values if there are non-planar faces on the border, because

- volumes() goes through tet4-24 to a tet3-x border
- volume() uses a tet3-u border

(file #33984)

gianluca de santis <gianlucadesi>
Group Member
Tue 12 May 2015 08:00:14 AM UTC, comment #5: 

following the comment of Gianluca
indeed the computation with tet-5 will be accuarate only with plane faces so indeed we need another solution.

 Though I need this to be improved, if you can give me the input to improve it with additional libraries please do, I will work on it and send it back to you. I have a code that needs many calls to volumes and this call slows it down a lot. So please give mehints on how to proceed

francesco <francio>
Group Member
Tue 12 May 2015 07:46:04 AM UTC, comment #4: 

I suggest we leave the volumes computation like it is (for now).
The educated user can use the approximate hex8->tet4-5 speedup himself when he wants.

Benedict Verhegghe <bverheg>
Group administrator
Mon 11 May 2015 07:30:45 AM UTC, comment #3: 

Hi Francesco,

- Regarding the element volumes:

if you use hex8->tet4-5 and you sum up the element volumes you do not get the correct total volume. The reason is that the face shared by two face-adjacent hex elements could be split in 2 tets using two different diagonals, producing a non-conformal tri3 mesh at the interface. If the shared quad4 face would be warping the tri3 meshes of the two hex may have gaps or intersect each other. The only subdivision which guarantees conformity is the tet-24. Of course, if you need a fast computation of all volumes to assess the element quality you do not care about warping, so an approximated volume would suffice, using the hex8->tet4-5, but you should be aware that the sum of the volumes is not the total volume

- Regarding the mesh volume:

what Benedict suggested is already the standard in mesh.py and it is indeed very fast:
3dmesh -> getbordermesh -> convert to tri3 -> surface_volume

gianluca de santis <gianlucadesi>
Group Member
Sun 10 May 2015 01:03:04 PM UTC, comment #2: 

Hi Benedict I dont even know what calpy is, so if you can start to include the libraries it would be easier.
Later we can start implementing our own and you could send me the functions to be converted in c++ to help you in this

francesco <francio>
Group Member
Sun 10 May 2015 08:43:02 AM UTC, comment #1: 

For the total volume, I would use

3dmesh -> getbordermesh -> convert to tri3 -> surface_volume

For individual volumes of elements, we can best implement direct formulas based on the shape functions. We have them available in calpy and calix(fortran), could be ported to pyformex, possibly even in accelerated form. Or we re-implement them in our C-library.

Benedict Verhegghe <bverheg>
Group administrator
Sat 09 May 2015 08:58:02 PM UTC, original submission:  

The default conversion for to tet4 is tet4-24
if we use tet4-5 teh computational cost will be reduced of 5 times. the problem is that the function uses as default the conversion to the base name 'tet4'.

Possible solutions
1 Either we add a new element conversion to all elements that calls the conversion with the minimum number of tets (eg. 'tet4-min' ) so if the element type is not 'tet4' we always convert to 'tet4-min'

2 add exception in the Mesh.volumes
if self.elName() == 'hex8':
    base = 'tet4-5'

3 use the inertia.surface_volume for all 3D elements. though I tested it and it works, the results are less accurate on an element base and it is a bit slower than the 'tet4-5' conversion

Additionally we could use the inertia.surface_volume to compute the total volume of the 3D mesh, which would much faster

francesco <francio>
Group Member


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by gianlucadesi (Posted a comment)
  • -email is unavailable- added by bverheg (Posted a comment)
  • -email is unavailable- added by francio (Submitted the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only logged-in users can vote.


    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2015-05-12 gianlucadesi Attached File- Added vtk_meshquality120515.py, #33984

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code