Making a mesh and calculating intersections of mesh boundaries with a straight line

Hi,

I am considering starting to program a so-called MOC (method of characteristics) solver for neutron transport. I have some programming experience, mainly in FORTRAN and Perl. The goal of my solver would be to try various approximations pertaining to MOC solvers, and thus the program does not need to be blisteringly fast. If certain experiments are deemed worthwhile, they will be reprogrammed in FORTRAN for speed.

I'd like to use Python because I think it allows an easier code development, at the expense of being a bit slower than FORTRAN, but much easier to program, and the post-processing of results into graphics will probably be easier. I have never written one single line of Python so it will be a steep learning curve :-))

For the solver, I need to solve a differential equation along a straight line (the path of neutron travel). The geometry is 2D (x-y plane). For instance, there will be mostly cylinders oriented along the z-axis, in a regular square or rectangular lattice. For the integration along a path, I need to know the intersections of the straight line with each of the geometrical shapes in the (x-y) plane (thus, for instance, the intersections of a straight line with circles). It should be noted that because of the angular dependence, the inside of a cylinder would be divided into several "sectors".

I have been thinking about some ways to program this and it seems like a lot of work. There is a famous Monte Carlo code (MCNP) which use the concepts of cells bounded by curved surfaces, but it would seem like a lot of work to re-implement such a scheme. I would also like to extend the solver to non-uniform lattices, which would require a more "generic" mesh generator. Thus I went online and found MeshPy as a mesh generator for 2D and 3D geometries.

My question is the following: is there a package available somewhere to easily calculate the intersections of a straight line with the mesh elements? Or is there an easy mathematical trick? I have zero experience with meshing, ray tracing etc, so I don't really know where to start to learn. The information I need is: for each mesh element, I need to know the length of the straight line in the mesh element. I also need to know the "material" of the mesh element. Furthermore, I need to know the surface area of each mesh element. Also, the algorithm to calculate the intersections has to be fast, because in a typical calculation there may be between 10 and 100 mesh elements, and there would be several 100 straight lines at various angles for which to determine the intersections.

Regards, Wilfred van Rooijen

‹ Problems running the test files Problems installing MeshPy on gentoo (AMD64) ›

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
Submitted by wvanrooijen on Wed, 28/10/2009 - 07:36.

Hi,

OK, I understand that MeshPy is a "wrapper" around two mesh generators, and the data structures returned to the python program are basically the same as coming from the original mesh generators. Once I get MeshPy up and running I'll see how to best use the data.

Cheers, Wilfred

Submitted by inducer on Tue, 27/10/2009 - 21:37.

First of all, I think you might be misunderstanding the purpose of MeshPy. It does not attempt to provide algorithms and data structures for everything you might want to do with unstructured meshes. Instead, it aims to be a interface to a number of mesh generators that allows you to get the mesh out and into your own data structures that are assumed to be better-adapted to the task at hand.

Second, I think the question you are asking depends on a lot of things, most of all the geometric situation of the rays and the intersecting bodies. IMO, only a benchmark of both will be able to tell for sure what's better.

Third, MeshPy so far is Py2.X-only, because Boost.Python (on which MeshPy is built) is Py2.X-only. There was a Google Summer-of-Code project to bring BP to Py3, so I anticipate that MeshPy+Py3 will be possible in the not-so-distant future.

HTH, Andreas

Submitted by wvanrooijen on Tue, 27/10/2009 - 20:57.

Hi,

I am not quite sure I understand what you mean. For my solver, I need to know not just the fact that a given segment of a ray is inside of a particalur mesh element, but the length of the segment on this mesh element.

I don't know what would be the easiest way to calculate the intersections:

  1. For each mesh element, determine which rays intersect this element and calculate the lengths

  2. For each ray, find out which mesh elements are crossed by the ray and determine the segment lengths.

I think the second option is cleaner. If you determine the ingoing and outgoing intersections of the ray with mesh element I, then if you know the neighbour J of this element, you only need to calculate the outgoing intersection with J. I assume that meshpy will give a data structure that expresses "boundary X is shared by mesh elements I and J".

Another question: does meshpy work with Python 3.0? I have bought the new edition of Mark Lutz' book "Learning Python" and it focuses on Python 3.0. On my gentoo computer I have installed Python 3.0 but noticed that nothing works with Python 3.0. I guess it has something to do with the syntax of Python 3.0 being somewhat different, and that many modules are not written using this syntax. I would like to start learning Python 3.0 as it is the newest version, but if meshpy only works with 2.6 that is fine as well.

Cheers, Wilfred

Submitted by Anonymous on Mon, 26/10/2009 - 22:32.

Some comments:

  • Use a cartesian mesh (or even an Octree) as a hash to quickly find elements near a point.
  • For ray-element intersection, consider only the planes containing the faces. Then use barycentric coordinates to check for containment. Once you've checked all 3 (4) surfaces, you have an intersection check.
  • For surface finding, see this code

HTH, Andreas