3. A simple 3D example


3.1. Drawing

We’ll start by drawing a volume defined by cross sections that are aligned with the cartesian coordinates. Open Inkscape, draw one shape, add a label and a phase to your shape (see previous tutorial) and save as an “Inkscape SVG” file (default file format). The name of the file should follow the format myfilename.orient.svg, where orient indicates the orientation of the drawing plane and is either HZ, NS or EW, respectively, horizontal, North-South or East-west vertical cross section.


Create a new layer. Layer names indicate the coordinate of the plane in the direction normal to the drawing plane (e.g. if in HZ orientation, the normal direction is the vertical). The layer name follows the format <letter numbers> without space, dots or commas. The letter indicate the sign, + (p) or - (m) and the number indicate the value. E.g. p0030, m250 etc... You can also have other layers starting with a # sign. Those layers won’t be read later by geomIO. You may want to add images in your file, be sure to link them rather than to embed them because SVG files with embedded images are cripplingly slow to read with matlab.

Copy the first shape (Ctrl+d) and transfer it to the other layer (Cmd+page up). Without adding new points, move the points defining the rabbit to fit another shape (e.g. a cat).


Later on, LaMEM will connect the dots between the shape on the first plane to the shape on the second (assuming they share the same label).

Next we can create a box to define our model boundary. Just create a rectangle on each layer. Don’t forget to give it a label and a phase.

Finally we are going to add a special layer called ‘Reference’. Only one path can go on this layer: a line that defines coordinate transformation from the drawing coordinate system to another coordinates (the ones of the numerical model for example). This line defines only translation and scaling transformation. Advanced coordinate transformation are gonna be discussed later (e.g spherical coordinates to map projections). This line doesn’t need a label but requires that you add a ‘CoordRef’ attribute with values: x1,y1 x2,y2, where x and y are coordinates in the new coordinate system.


3.2. Creating volumes

Drawing was a lot of work! Fortunately now you can relax and let geomIO do the work. One last thing though. As for the 2D example, initialize and run geomIO:

opt                 = geomIO_Options();
opt.inputFileName   = ['Test.HZ.svg'];
[~, Volumes]        = run_geomIO(opt);

You can visualize the volumes created directly in matlab:


geomIO can also output a set a vtk files (in the directory ./Output/ ), that you can use to visualize the volumes in Paraview. For optimal viewing in Paraview you can slightly shift the volumes so that they don’t overlap. The modified code looks like:

opt                 = geomIO_Options();
opt.inputFileName   = ['Test.HZ.svg'];
opt.writeParaview   = true;
opt.shiftPVobj      = [0 0 1];
[~, Volumes]        = run_geomIO(opt);

You can find details of the opt.shiftPVobj option here.

3.3. Creating Polygon files

Now that we have volumes a nice thing would be to use them to assign material properties to our numerical model. In order to efficiently assign material properties to markers in 3D we will now output a Polygon file. A Polygon file is a binary file containing 2D slices of the 3D volumes. These slices then allow one to perform 2D points in polygon check rather than more expensive 3D points in volume checks. Details of the Polygon file format and pseudo code for reading them in a parallel fashion (using petsc) are available here. To create a Polygons file you need to provide three arrays containing the structure of markers in each direction. For example:

opt.setup.marker_x = 0:1000; opt.setup.marker_y = 0:1000; opt.setup.marker_z = 0:1000;

And to specify the option opt.writePolygons = true. If you use LaMEM you can read the marker information directly from your LaMEM input file.

opt.readLaMEM           = true;
opt.LaMEMinputFileName  = './RabbitAndCat.dat';
opt.writePolygons       = true;

If you drew a a shape with many segments you might experience that creating polygons is slow. You can set up the number of points used to draw each segment with the option opt.DrawCoordRes = aNumber. The default is 21, but if all your segments are relatively short you can get a good resolution with as low as ~5 points per segments.