search

svalinn / r2s-act

Rigorous 2 Step Activation Workflow
4 stars 7 forks source link

Calculate volumes for arbitrary tetrahedra and hexahedra #40

Closed erelson closed 11 years ago

erelson commented 11 years ago

This is needed for both ALARA, and for calculating the source normalization factor prior to photon sampling.

For arbitrary tetrahedra, I am currently using the method from Wikipedia: http://en.wikipedia.org/wiki/Tetrahedron#Volume. One note is that the abs() function may be unnecessary, depending on the vertex ordering that MOAB uses?

For arbitrary hexahedrons, the method described by Grandy in 'Efficient Computation of Volume of Hexahedrons' (also in Zotero) could be used, with adjustments made for the different vertex ordering used in MOAB. (this is for hexahedra with non-planar facets)

MOAB/mbsize already do similar volume calculations, so their methods should be investigated as well.

gonuke commented 11 years ago

See MOAB/iTaps methods for this: I think they are called ''measure"

erelson commented 11 years ago

Looks like the relevant code is measure_volume() in MOAB/src/trunk/tools/dagmc/DagMC.cpp. This function supports only triangular-faceted polyhedra.

For now, I will use the code I have working within Fortran (a local branch), which makes testing and inclusions simpler to deal with.

gonuke commented 11 years ago

I'm not sure why it lives in MOAB/src/trunk/tools/ but measure.[hc]pp does it all (I think).

erelson commented 11 years ago

Ahh, I see that now. I was tracing from what DAG-MCNP calls, and so did not stumble upon measure().

The challenge in my case is that measure() is not part of the iTaps API, and is thus trickier to call from Fortran.

erelson commented 11 years ago

Current implementation for hexahedra duplicates the approach in measure() - calculate volumes of 5 tetrahedra within the hex. Since the hexahedra are limited to being parallelepipeds (#46), there may be a slightly more efficient volume method to be found, since parallelepipeds are a specific subset of quad-faceted hexahedra.

erelson commented 11 years ago

Determined today that the tetrahedral volume formula from Wikipedia, V=abs( (a-d) . ((b-d) x (c-d)) ) / 6 seems to work except when it gives a volume of zero (namely a=1,0,0; b=0,1,0; c=0,0,1; d=1,1,1;). I will look more closely at the method used in measure.cpp. I got Mathematica's cross product operator confused, and so was incorrectly calculating the volume.

erelson commented 11 years ago

In Python, the PyTaps example gives us the necessary code: http://pythonhosted.org/PyTAPS/tutorial/example.html?highlight=volume

The three relevant functions have been added to write_alara_geom.py.