# 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);
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;


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;
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

Doxygen generated

Updated