# Wiki

# meshlib / Home

## Overview

*MeshLib* is a C++ library for simple calculations with
polygon surface and volume meshes.
The library includes classes for basic geometry primitives
like *Line*, *Segment*, *Plane*, *Triangle*, *Quadrangle*,
*BoundingBox*, *Polyhedron*, *Transformation*, as well as
surface mesh classes *SurfaceMesh* and *VolumeMesh*, and
a class for a set of surface meshes *SurfaceMeshSet*.
The library implements algorithms for intersection tests between
a line segment and a plane/triangle/quadrangle and for
surface integration along polygons.
The *SurfaceMeshSet* class provides a grid
accelerator for intersection tests with many polygons.

## Examples

The library has three major parts:
geometrical primitives, surface and volume meshes
and a set of surface meshes with grid accelerator.
Besides, there are help classes used internally, e.g.
*BaseId*, *ContainerOfIdItems*, *CartesianGrid*,
*GridCache*, *Integration*, *RandomGenerator*.

All classes of the library are in the namespace *meshlib*.
In the examples hereafter it is assumed that the following
namespace are used: *linalg*, *meshlib*, *std*.

### Geometrical primitives

The primitive classes include:
*Point2D*, *Point3D*, *Line*, *Segment*, *Plane*, *Triangle*,
*Quadrangle*, *BoundingBox*, *Polyhedron*, *Transformation*.

**Point.** Cartesian points *Point2D* and *Point3D* are straightforward to use.
They wrap *Vector* class from LinAlg
and provide convenience methods for directly setting the coordinates and
for getting distance and tolerance checks.

```
Point3D p1, p2;
p1[1] = 2;
p2.getCoordinates() = 2;
double distance = p1.getDistance(p2);
if (p1.isEqual(p2, 1e-3))
cout<<"oops";
```

**Line.** A line is given by a point on the line and a tangential vector. The
tangential vector is not necessarily normalized. There are methods
to find a distance between a point and the line and to transform line
with a given transformation.

```
Point3D point, a;
Vector3D direction;
Line line(point, direction);
line.getDistance(a);
Transformation t;
line.transform(t);
```

**Segment.** A line segment is given by its start and stop point. The *Segment*
class overloads the *getDistance* method taking the segment
limits into account.

```
Point3D start, stop, a;
Segment segment(start, stop);
segment.getDistance(a);
```

**Plane.** A plane is given by a point in the plane and a normal vector. The normal
vector is normalized to 1. The *Plane* class has methods to
find an intersection with a line segment, to project a point into the plane
and to transform the plane with a given transformation.

```
Point3D point, a, projection, intersection;
Vector3D normal;
normal[2] = 1.0;
Segment segment;
Plane plane(point, normal);
double distance = plane.project(a, projection);
double t;
if (plane.intersect(segment, intersection, t)
cout<<"intersects";
Transformation t;
plane.transform(t);
```

**Bounding box.** A bounding box is considered axis aligned.
An existing bounding box
can be updated by points, polygons or meshes. Two bounding
boxes can be combined to produce a bounding box including
both of them. There are methods to check if another bounding
box or a point is inside. It is also possible to find intersections
of a line with bounding box.

```
BoundingBox box1, box2;
vector<Point3D > vertices;
box1.update(vertices);
Triangle triangle;
box2.update(triangle);
box1.isInside(box2);
box1.combine(box2);
Point3D point;
box1.isInside(point);
Line line;
double t1, t2;
box1.intersect(line, t1, t2);
```

**Transformation.** A space transformation
consists of a matrix, scaling factor and a translation.

```
v = scale * M * u + shift
```

Different parts of the transformation can be set independently without any restrictions. There are also methods to create widely used types of transformations, e.g. rotation, center of symmetry, identity, inverse.

```
Transformation t;
Matrix3D mat;
t.setMarix(mat);
t.setScale(2.0);
t.makeUnit(); // identity
t.makeTranslation(shift); // only translation
Vector3D ax1, ax2;
t.makeRotation(ax1, ax2); // rotation from ax1 to ax2
t.makeRotation(ax1, 0.5); // rotation around ax1 by 0.5 ccw
```

*Transformation* has methods to be applied to a point or a vector.
Other geometrical primitives usually have an own method for
transformation.

```
Transformation t;
Point3D point;
Vector3D vec;
t.apply(point);
t.apply(vec);
t.applyWithoutTranslation(vec);
Line line;
line.transform(t);
```

It is also possible to multiply transformations in a normal mathematical sense:

```
Transformation t1, t2, t3;
t3 = t1*t2;
t1 *= t3;
```

**Polygon** A polygon is assumed to be flat and convex. All polygons have an abstract
base *AbstractPolygon* defining the necessary interface and providing some
additional convenience methods like *intersect* for segment-polygon intersection,
*triangulate* for a simplistic triangulation, finding some help integrals of
vectors *(x,y,z)* and *(x^2, y^2, z^2)*, etc. Integration of *(x,y,z)*
gives geometrical center of a polygon, whereas integration *(x^2, y^2, z^2)*
can be used to computed volume of a volume enclosed by a set of polygons.
Triangulation is done in a fan way. The polygon required interface include
methods for getting vertices, a random point inside, plane of the polygon,
surface area, changing polygon orientation, finding polygon-plane intersection,
measuring distance to a point, determining if a point is inside polygon,
finding surface integrals, applying a geometrical transformation, creating a clone
of a polygon. There are two implementations of the interface available:
*Triangle* and *Quadrangle*. The *Triangle* class
checks if the point is inside by using skew basis consisting of two
adjacent side vectors.

```
Point3D v1, v2, v3, v4;
Triangle tri(v1, v2, v3);
Quadrangle quad(v1, v2, v3, v4);
vector<Point3D> chain;
chain.resize(10);
vector<Triangle> triset;
AbstractPolygin::triangulate(chain, triset); // fan triangulation
AbstractPolygon *poly = &tri;
Point3D point, intersection;
double t;
Segment segment;
// find segment-polygon intersection: find segment-plane intersection
// and check if the intersection point is inside polygon
if (poly->intersect(segment, intresection, t)) // segment-polygon intersection
cout<<"intresects";
const Vector3D &normal = poly->getNormal(); // normal of the polygon plane
Plane plane;
poly->getPlane(plane); // polygon plane
cout<<poly->getArea(); // surface area of the polygon
vector<Point3D> res;
int numPoints = poly->intersect(plane, res); // plane-polygon intersection
double dist = poly->getDistance(point);
if (poly->isInside(point)) // point is assumed to be on the polygon plane
cout<<"point inside";
RandomGenerator rand;
poly->getRandomPoint(rand, point); // a random point in the polygon
Transformation t;
poly->transform(t);
AbstractPolygon *poly2 = poly->clone();
```

In addition to the polygon interface, the *Triangle* class has a method
for breaking the triangle into two subtriangles along the longest
side if its length is longer than a given limit. Note that the original
triangle is not changed.

```
Triangle tri(v1, v2, v3);
Triangle tri1, tri2;
if (tri->split(tri1, tri2, 10))
cout << "was splited";
```

It is possible to find an integral over a polygon surface area.
Integration can be done either for a scalar function of position
or for a vector function. In the latter case the normal
multiplication is not performed, i.e. the returned result is
a vector. To integrate one has first to declare a class
inheriting *ScalarField* or *VectorField* correspondingly and
implementing *operator()*.

```
class ScalarExmpl : public ScalarField
{
public:
void operator()(const Vector3D &r, double &res) const
{ res = 1.0; }
};
class VectorExmpl : public VectorField
{
public:
void operator()(const Vector3D &r, Vector3D &res) const
{ res = r; }
};
ScalarExmpl sc;
VectorExmpl vec;
double scalres;
Vector3D vecres;
Triangle tri;
tri->getSurfaceIntegral(sc, 1e-5, 1e-6, scalres);
tri->getSurfaceIntegral(vec, 1e-5, 1e-6, vecres);
```

The integration methods require absolute and relative tolerance of the result. Internally trapezoidal integration with a decreasing step is performed until the difference of two consecutive approximation is below the required accuracy:

```
|R_i - R_{i-1}| <= absTolerance + relTolerance * |R_i|
```

For a triangle an integration is performed in a skew basis based upon two sides of the triangle.

**Polyhedron**. A polyhedron is a finite volume limited by flat polygon
faces. In spite of being a primitive geometrical body for *VolumeMesh*,
a polyhedron object itself contains a *SurfaceMesh* describing its
boundary. For this reason, examples for using polyhedrons are given
later.

### Meshes

Meshes are collections of corresponding geometrical primitives:
*SurfaceMesh* is a collection of polygons and *VolumeMesh*
is a collection of polyhedrons.

**Surface mesh.** A surface mesh is a collection of polygons that
can represent, for example, a single body. There are methods for
construction of meshes, combining meshes, creating a copy of them,
refining meshes, integration, and writing/reading to files.
A mesh can be set to polygons created externally or to a triangulation
of a polygon line.

```
vector<AbstractPolygon*> list;
SurfaceMesh m1(list); // provide externally created polygons
vector<Point3D> chain;
SurfaceMesh m2(chain); // set to a fan triangulation of the polygon chain
vector<Point3D> chain_in, chain_out;
SurfaceMesh m3(chain_in, chain_out); // triangles between two polygon chains
```

A mesh can be constructed or set by specifying a list of vertices of all polygons and polygons by references into the vertices. The vertices and polygons can be provided in different ways: as vectors or as a map with additional ids, etc.

```
vector<Point3D> nodes; // all vertices
vector<size_t> elements; // flat list of all polygons as references
vector<short int> numVertices; // number of vertices per polygon
SurfaceMesh m1(nodes, elements, numVertices); // number of polygons = numVertices.size()
map<size_t, Point3D> nodes_map; // map vertex id - vertex coordinates
map<size_t, vector<size_t> > elements_map; // map polygon id - polygon vertices
SurfaceMesh m2(nodes_map, elements_map); // number of polygons = elements_map.size()
...
SurfaceMesh *m3 = m2.triangulate().release(); // convert all polygons into triangles
SurfaceMesh *m4 = m3->refineMesh(10).release(); // refine triangles sizes
m4->deleteSmallAreaPolygons(10);
```

It is also possible to construct common types of geometrical bodies with a given resolution.

```
SurfaceMesh m;
Point3D v1, v2;
m.makeRectangularBox(v1, v2); // cube
vector<Point3D> face1, face2;
m.makeConnected(face1, face2); // face 1 and 2 - top bottom + faces on every side between
m.makeSphere(1, 11, 11); // sphere of radius 1 and with 10 polygons in each angle direction
m.makeCylinder(1, 0, 1, 11, 11); // cylinder, R = 1, z = [0, 1], 10 elements in both directions
m.makeCircle(1, 11); // plain circle of radius 1
```

It is possible to do simple operations with a mesh.

```
class ScalarExmpl : public ScalarField
{
public:
void operator()(const Vector3D &r, double &res) const
{ res = 1.0; }
};
class VectorExmpl : public VectorField
{
public:
void operator()(const Vector3D &r, Vector3D &res) const
{ res = r; }
};
ScalarExmpl sc;
VectorExmpl vec;
double scalres, vecres2;
Vector3D vecres;
SurfaceMesh m;
m->getSurfaceIntegral(sc, 1e-5, 1e-6, scalres);
m->getSurfaceIntegral(vec, 1e-5, 1e-6, vecres); // no normal multiplication
m->getSurfaceIntegral(vec, 1e-5, 1e-6, vecres2); // with normal multiplication - \int \vec{F}\cdot\vec{dS}
m->getAres(); // total surface area
AbstractPolygon *poly;
if (m->containsPolygon(poly))
cout<<"polygon in the mesh"; // the same polygon, i.e. the same object
Transformation t;
m->transform(t);
```

It is possible to save a mesh into a file and to read from a file.
Two files are created *"filename.nodes"* and *"filename.elements"*.
The first one contains list of nodes in the form *(id x y z)*,
the second one contains polygons referencing the nodes *(id ref1 ref2 ref3 ...)*.

```
SurfaceMesh m1;
m1->dumpToFile("mymesh");
SurfaceMesh m2;
m2->readFromFile("mymesh");
```

A surface mesh keeps polygons by pointers and frees the memory in the destructor.
It is possible to add/remove external polygons and to exchange polygons
between meshes, see also *ContainerOfIdItems*.

```
SurfaceMesh m1, m2;
AbstractPolygon *poly;
m1.add(auto_ptr<AbstractPolygon>(poly)); // m1 mesh now owns the polygon
poly = new Triangle();
m1.add(&poly); // m1 has a copy of the polygon, it is safe to delete it here
delete poly;
m1.add(m2); // m1 now has a copy of all polygons in m2
m1.transfer(m2); // m1 now has all polygons from m2 and m2 is empty
m1.remove(10); // remove polygon 10
m1.drop(); // drop all polygons without deleting them, m1 is empty
m1.makeCircle(...);
m1.clear(); // delete all polygons
m1 = m2; // create copy of polygons in m2
```

**Polyhedron.** A polyhedron is described by its surface boundary. This can
be supplied as a straight list of polygons or as *SurfaceMesh*.
A *Polyhedron* instance keeps its surface by pointer and deletes it in
the destructor. Polyhedron reorients all boundary faces in such
a way that their normals point outwards from the volume.

```
SurfaceMesh *m1, *m2;
Polyhedron p1(auto_ptr<SurfaceMesh>(m1)); // p1 owns the mesh
m1 = p1.dropSurfaceMesh().release(); // p1 does not own m1 any more
p1.setSurfaceMesh(&m2); // p1 has a copy of m2
```

A polyhedron instance can be used to get surface integrals in the same way as its surface mesh boundary. In addition there are methods to get volume, typical size and center of the polyhedron.

```
Polyhedron p1;
double vol = p1.getVolume();
double s = p1.getTypicalSize(); // an order of magnitude estimation of the size
Point3D center = p1.getCenter(); // center of mass
Transformation t;
p1.transform(t);
```

**VolumeMesh.** Volume mesh is a collection of polyhedrons forming, for example,
volume of a body. There are different methods for creating a volume mesh.
A mesh can be set by a list of polyhedrons.

```
vector<Polyhedron*> vols;
VolumeMesh m(vols); // the new mesh owns all the polyhedrons
```

It is also possible to set up a mesh by specifying boundary polygons
and separate polyhedrons as references into faces. The polygon
faces should be provided as a single *SurfaceMesh* instance
containing. Polyhedrons can be supplied in several forms, which
are similar to the ways of constructing a *SurfaceMesh* instance.

```
SurfaceMesh faces; // polygon faces for all polyhedrons of the volume mesh
vector<short int> numFaces; // number of boundary faces per polygon
// create a volume mesh where all faces are references consecutively
VolumeMesh m(faces, numFaces); // number of polyhderons == numFaces.size()
// a flat array of references into faces
vector<size_t> polyhedrons; // references into faces, can be not ordered
// numFaces is used to split polyhedrons vector into single instances
m.setMesh(faces, polyhedrons, numFaces);
// if number of faces per polyhedron is constant
m.setMesh(faces, polyhedrons, 6);
vector<vector<size_t > > polyhedrons_2d; // 2d array of polyhedron references
m.setMesh(faces, polyhedrons_2d); // number of polyhedrons = polyhedrons_2d.size()
map<size_t, vector<size_t> > polymap; // map polyhedron id -> list of boundary faces
m.setMesh(faces, polymap); // number of polyhedrons == polymap.size()
```

Is is also possible to create volume mesh of common geometrical primitives.

```
// polyhedrons filling an axis aligned box with given resolution
Point3D v1, v2;
m.makeRectangularBox(v1, v2, 11, 11, 11); // 10 intervals in x,y,z, total 1000
// make a cylindrical shell, i.e. filling cylinder between minimal and maximal radii
// 4 cells in radius, 10 cells in angle, 10 cells in z
m.makeCylinderShell(0.5, 1.0, 0, 1.0, 5, 11, 11); // 0.5 < r < 1, 0 < z < 1
// make a spherical shell, i.e. filling sphere between given radii
m.makeSphereShell(0.5, 1.0, 11, 11, 11); // 0.5 < r < 1.0, 10 cells in all directions
// it is possible to create a set cubes filling an arbitrary volume,
// where the volume is specified by a predicate. A predicate
// is given the current polyhedron and should decide if it belongs to the
// volume or not.
bool sphere(const Polyhedron &box)
{
Point3D center; // by default (0,0,0)
return (box.getCenter().getDistance(p) < 1.0);
}
// min/max points of a bounding box where cubes are generated for the predicate check
Point3D p1, p2;
p1.getCoordinates() = -1; // point (-1,-1,-1)
p2.getCoordinates() = 1; // point (1,1,1)
m.makeBoxedVolume(p1, p2, 31, 31, 31, sphere); // a unit sphere with 30 per direction
cout << m.getVolume(); // total volume of the mesh
```

A *VolumeMesh* instance owns its polyhedrons in the same way as a
*SurfaceMesh* instance. All exchange/combine operations are the same as with
*SurfaceMesh*, see examples above.
It is possible to save a volume mesh into a file and to read from a file.
Three files are created *"filename.nodes"*, *"filename.elements"* and
*"filename.solids"*.
The first one contains list of nodes in the form *(id x y z)*,
the second one contains polygons referencing the nodes *(id ref1 ref2 ref3 ...)*,
and the third one contains polyhedrons referencing the boundary
faces *(id face1 face2 face3 face4 ...)*.

```
VolumeMesh m1;
m1.dumpToFile("mymesh");
VolumeMesh m2;
m2.readFromFile("mymesh");
Transformation t;
m2.transform(t);
```

### Set of surface meshes

*SurfaceMeshSet* is a collection of *SurfaceMesh* instances with
a grid accelerator to speed up segment-mesh intersection testing.
A *SurfaceMeshSet* instance subdivides the bounding box
of all meshes into a given number of volume cells, voxels, and
for every voxel keeps a list of polygons that cross it, or are
inside. On a request for the closest segment-polygon test
only polygons in the volumes covered by the segment are considered.
In such a way the intersection testing can be significantly
accelerated for large polygon meshes. To create
a *SurfaceMeshSet* instance one has to specify the number of
volume cells in all three directions, whereas the space that
is subdivided is fitted to the region of added meshes.
If a new mesh is added to the set, which is beyond the previously
partitioned space, the space is extended and the all polygons
are reassigned. The size of the current space, volume cell size
can be accessed by the corresponding methods.

```
SurfaceMeshSet set(10, 10, 10); // 10 cells in every direction, total 1000
SurfaceMesh *m1, *m2;
// add meshes, the polygons of the meshed are assigned to volume cells
set.add(m1); // set owns the mesh
set.add(*m2); // set has a copy
cout<<set.getBoundingBox(); // get the total partitioned space
cout<<set.getGridSteps();
```

A *SurfaceMeshSet* instance owns added meshes and will delete
them in destructor. The exchange/transfer operations are the
same as for *SurfaceMesh*, see exampled above.
Segment-mesh intersection testing is straightforward,
the intersection closest to the start of the segment is returned.

```
SurfaceMeshSet set(10, 10, 10); // 10 cells in every direction, total 1000
Segment segment;
Point3D intersection;
const AbstractPolygon *poly = set.intersect(segment, intersection);
if (poly)
{
cout<<"intersects";
const SuraceMesh *m1 = set.findComponent(poly); // mesh the polygon belongs to
}
```

It is possible to check many segments at with one request. If Klubok library is present such a request is parallelized. It is also possible to perform the intersection cheks with incrementing an event counter in encountered polygons.

```
SurfaceMeshSet set(10, 10, 10); // 10 cells in every direction, total 1000
vector<Segment> segments;
vector<Point3D> intersections;
vector<const AbstractPolygon* > polys;
set.intersect(segments, polys, intersections);
// with event counter
set.intersect_withEventCounter(segments, polys, intersections);
// check for the number of hits in a polygon
AbstractPolygon *polygon1;
cout<<polygon1->getEventCounter();
```

### Help classes

Help classes are used internally in the library.

*BaseId* - base class having an id and an event counter,
it is used for polygons, meshes, etc.

*CartesianGrid* - is a 3d grid in an axis aligned box,
can be resized, it is used for grid accelerator in
*SurfaceMeshSet*.

*GridCahce* - keeps a list of geometrical elements
that are at least partially inside volume cells of a 3d
grid, it is used for grid accelerator in *SurfaceMeshSet*.

*ContainerOfIdItems* - a template container class that owns
its elements and has operations for their exchange and transfer,
it is used in meshes and mesh set.

*Integration* - several template classes implementing integration.

*RandomGenerator* - a wrapper to produce random numbers, if
GSL is present a Mersenne
twister is used, otherwise the standard C++ generator is used.

## Build

For installation of the library,
LinAlg is required.
For parallelization of multiple segments testing in
*SurfaceMeshSet*,
Klubok is necessary,
if it is not present than the simple looping is done.
For a better random number generator
GSL is recommended,
if it is not present the standard C++ one is used.

Installation is standard for cmake projects.

```
>>> mkdir build && cd build
>>> cmake ..
>>> make
>>> make install
```

To build tests Boost.Tests is required.

```
>>> make tests
>>> tests/tests --log_level=test_suite
```

## API

Updated