Index: cpp/frams/Makefile-GDK
===================================================================
--- cpp/frams/Makefile-GDK	(revision 187)
+++ cpp/frams/Makefile-GDK	(revision 191)
@@ -5,5 +5,5 @@
 include frams/Makefile-GDK-files
 
-GDK_TESTS=genoconv_test gdk_test geno_test neuro_test loader_test serial_test multiline_f0_test f0_variants_test full_props part_shapes neuro_layout_test genooper_test
+GDK_TESTS=genoconv_test gdk_test geno_test neuro_test loader_test serial_test multiline_f0_test f0_variants_test full_props part_shapes neuro_layout_test genooper_test geometry_apices_test geometry_info_test geometry_surface_test geometry_volume_test
 
 gdk_tests: $(GDK_TESTS)
@@ -50,7 +50,19 @@
 	$(CXX) $(GENOOPER_TEST_OBJS) -o $@
 
+geometry_apices_test: $(GEOMETRY_APICES_TEST_OBJS)
+	$(CXX) $(GEOMETRY_APICES_TEST_OBJS) -o $@
+
+geometry_info_test: $(GEOMETRY_INFO_TEST_OBJS)
+	$(CXX) $(GEOMETRY_INFO_TEST_OBJS) -o $@
+
+geometry_surface_test: $(GEOMETRY_SURFACE_TEST_OBJS)
+	$(CXX) $(GEOMETRY_SURFACE_TEST_OBJS) -o $@
+
+geometry_volume_test: $(GEOMETRY_VOLUME_TEST_OBJS)
+	$(CXX) $(GEOMETRY_VOLUME_TEST_OBJS) -o $@
+
 #############################################
 
-ALL_OBJS=$(GDK_TEST_OBJS) $(GENOCONV_TEST_OBJS) $(GENO_TEST_OBJS) $(NEURO_TEST_OBJS) $(LOADER_TEST_OBJS) $(SERIAL_TEST_OBJS) $(MULTILINE_F0_OBJS) $(F0_VARIANTS_OBJS) $(FULL_PROPS_OBJS) $(PART_SHAPES_OBJS) $(NEURO_LAYOUT_TEST_OBJS) $(GENOOPER_TEST_OBJS)
+ALL_OBJS=$(GDK_TEST_OBJS) $(GENOCONV_TEST_OBJS) $(GENO_TEST_OBJS) $(NEURO_TEST_OBJS) $(LOADER_TEST_OBJS) $(SERIAL_TEST_OBJS) $(MULTILINE_F0_OBJS) $(F0_VARIANTS_OBJS) $(FULL_PROPS_OBJS) $(PART_SHAPES_OBJS) $(NEURO_LAYOUT_TEST_OBJS) $(GENOOPER_TEST_OBJS) $(GEOMETRY_APICES_TEST_OBJS) $(GEOMETRY_INFO_TEST_OBJS) $(GEOMETRY_SURFACE_TEST_OBJS) $(GEOMETRY_VOLUME_TEST_OBJS)
 
 include frams/Makefile-maintain
Index: cpp/frams/Makefile-GDK-files
===================================================================
--- cpp/frams/Makefile-GDK-files	(revision 187)
+++ cpp/frams/Makefile-GDK-files	(revision 191)
@@ -2,5 +2,5 @@
 
 # ALL_DIRS is later expanded by the shell, no spaces/newlines allowed, or it breaks
-ALL_DIRS={common,frams,frams/canvas,frams/config,frams/errmgr,frams/genetics,frams/genetics/f0,frams/genetics/f1,frams/genetics/f2,frams/genetics/f3,frams/genetics/f4,frams/genetics/f5,frams/genetics/f6,frams/genetics/f7,frams/genetics/f8,frams/genetics/f9,frams/genetics/fF,frams/model,frams/neuro,frams/neuro/impl,frams/param,frams/test,frams/util,frams/vm/classes,frams/virtfile,frams/_demos}
+ALL_DIRS={common,frams,frams/canvas,frams/config,frams/errmgr,frams/genetics,frams/genetics/f0,frams/genetics/f1,frams/genetics/f2,frams/genetics/f3,frams/genetics/f4,frams/genetics/f5,frams/genetics/f6,frams/genetics/f7,frams/genetics/f8,frams/genetics/f9,frams/genetics/fF,frams/model,frams/neuro,frams/neuro/impl,frams/param,frams/test,frams/util,frams/vm/classes,frams/virtfile,frams/_demos,frams/model/geometry,frams/_demos/geometry}
 
 GENMANF9=frams/genetics/f9/oper_f9.o
@@ -18,4 +18,5 @@
 GDK_OBJS=frams/util/list.o frams/util/advlist.o frams/param/param.o frams/util/sstring.o frams/util/sstringutils.o frams/util/3d.o frams/vm/classes/3dobject.o frams/model/model.o frams/model/modelparts.o frams/neuro/neurolibrary.o frams/genetics/geno.o frams/genetics/genoconv.o frams/util/extvalue.o frams/vm/classes/collectionobj.o frams/util/hashtable.o common/framsg.o common/stl-util.o frams/util/callbacks.o frams/param/syntparam.o frams/util/multirange.o frams/util/multimap.o frams/param/paramtabobj.o frams/errmgr/errmanager.o frams/param/paramobj.o frams/genetics/oper_fx.o common/nonstd_math.o frams/errmgr/stderrors.o common/Convert.o frams/util/rndutil.o
 
+GEOMETRY_OBJS=frams/model/geometry/meshbuilder.o frams/model/geometry/modelgeometryinfo.o frams/model/geometry/geometryutils.o
 
 #### gdk_tests
@@ -46,2 +47,10 @@
 
 GENOOPER_TEST_OBJS=frams/_demos/genooper_test.o frams/virtfile/virtfile.o frams/errmgr/stdouterr.o $(GDK_OBJS) $(GENOCONV_GDK_OBJS) $(GENMAN_GDK_OBJS)
+
+GEOMETRY_APICES_TEST_OBJS=frams/_demos/geometry/apices_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/virtfile/virtfile.o frams/param/multiparamload.o frams/virtfile/stdiofile-autoinit.o frams/virtfile/stdiofile.o common/nonstd_stdio.o $(GDK_OBJS) $(GENOCONV_GDK_OBJS) $(GEOMETRY_OBJS)
+
+GEOMETRY_INFO_TEST_OBJS=frams/_demos/geometry/info_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/virtfile/virtfile.o frams/param/multiparamload.o frams/virtfile/stdiofile-autoinit.o frams/virtfile/stdiofile.o common/nonstd_stdio.o $(GDK_OBJS) $(GENOCONV_GDK_OBJS) $(GEOMETRY_OBJS)
+
+GEOMETRY_SURFACE_TEST_OBJS=frams/_demos/geometry/surface_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/virtfile/virtfile.o frams/param/multiparamload.o frams/virtfile/stdiofile-autoinit.o frams/virtfile/stdiofile.o common/nonstd_stdio.o $(GDK_OBJS) $(GENOCONV_GDK_OBJS) $(GEOMETRY_OBJS)
+
+GEOMETRY_VOLUME_TEST_OBJS=frams/_demos/geometry/volume_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/virtfile/virtfile.o frams/param/multiparamload.o frams/virtfile/stdiofile-autoinit.o frams/virtfile/stdiofile.o common/nonstd_stdio.o $(GDK_OBJS) $(GENOCONV_GDK_OBJS) $(GEOMETRY_OBJS)
Index: cpp/frams/_demos/geometry/apices_test.cpp
===================================================================
--- cpp/frams/_demos/geometry/apices_test.cpp	(revision 191)
+++ cpp/frams/_demos/geometry/apices_test.cpp	(revision 191)
@@ -0,0 +1,57 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "geometrytestutils.h"
+#include <frams/model/geometry/meshbuilder.h>
+#include <frams/model/model.h>
+#include <frams/util/3d.h>
+#include <stdio.h>
+
+void test(Model &model, const double density)
+{
+	// Creating and preparing result Model object.
+	Model resultModel;
+	resultModel.open();
+	GeometryTestUtils::addAnchorToModel(resultModel);
+	
+	// Creating instance of Iterator class (MeshBuilder::ModelApices in this case). Object is
+	// placed in stack memory, thus there is no heap memory allocation and creation is fast.
+	MeshBuilder::ModelApices iterator(density);
+	
+	// Once created, iterator can be used many times, but must be initialized before every time.
+	iterator.initialize(&model);
+	
+	// Creating iterator output object in stack memory.
+	Pt3D point;
+	
+	// Method tryGetNext checks if next point exists and optionally updates fields of Pt3D object.
+	while (iterator.tryGetNext(point))
+	{
+		// Processing points created by iterator. In this case, they are added to result model as
+		// small spheres.
+		GeometryTestUtils::addPointToModel(point, resultModel);
+	}
+	
+	// After Iterator initialization there are two alternative ways of using it.
+	//
+	// First adds all points to the specified list:
+	// SListTempl<Pt3D> list;
+	// iterator.addAllPointsToList(list);
+	//
+	// Second executes specified operation for each point:
+	// Operation func;
+	// iterator.forEach(func);
+	// In this case, Operation is derived class of MeshBuilder::Callback class.
+	
+	// Finishing result Model and printing its genotype.
+	resultModel.close();
+	puts((const char*)resultModel.getF0Geno().toString());
+}
+
+int main(int argc, char *argv[])
+{
+	return GeometryTestUtils::execute(
+		"Creates new model built of small spheres at input model apices and prints its genotype.",
+		argc, argv, test);
+}
Index: cpp/frams/_demos/geometry/geometrytestutils.cpp
===================================================================
--- cpp/frams/_demos/geometry/geometrytestutils.cpp	(revision 191)
+++ cpp/frams/_demos/geometry/geometrytestutils.cpp	(revision 191)
@@ -0,0 +1,307 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "geometrytestutils.h"
+
+#include "../genotypeloader.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+int printGenotypesList(const char *file)
+{
+	long count = 0;
+	long totalSize = 0;
+	MiniGenotypeLoader loader(file);
+	MiniGenotype *genotype;
+	
+	while (genotype = loader.loadNextGenotype())
+	{
+		count++;
+		totalSize += genotype->genotype.len();
+		
+		fprintf(stderr, "%d. (%6d chars) %s\n", count, genotype->genotype.len(),
+			(const char*)genotype->name);
+	}
+	
+	if (loader.getStatus() == MiniGenotypeLoader::OnError)
+	{
+		fprintf(stderr, "Error: %s\n", (const char*)loader.getError());
+		return 2;
+	}
+	else
+	{
+		fprintf(stderr, "\ntotal: %d items, %d chars\n", count, totalSize);
+		return 0;
+	}
+}
+
+class TestInvoker
+{
+	public:
+		virtual void operator()(Model &model) = 0;
+};
+
+int executeTestUsingLoadedModel(const char *file, const char *genoId, TestInvoker &test)
+{
+	const char* genoName = genoId;
+	const int genoIndex = isdigit(genoId[0]) ? atol(genoId) : 0;
+	long count = 0;
+	MiniGenotypeLoader loader(file);
+	MiniGenotype *genotype;
+	
+	while (genotype = loader.loadNextGenotype())
+	{
+		count++;
+		
+		if ((genoIndex == count) || (strcmp((const char*)genotype->name, genoName) == 0))
+		{
+			Model model(genotype->genotype);
+			
+			if (!model.isValid())
+			{
+				fprintf(stderr, "Cannot build Model from this genotype!\n");
+				return 4;
+			}
+			
+			test(model);
+			return 0;
+		}
+	}
+	
+	if (loader.getStatus() == MiniGenotypeLoader::OnError)
+	{
+		fprintf(stderr, "Error: %s\n", (const char*)loader.getError());
+		return 2;
+	}
+	else
+	{
+		fprintf(stderr, "Genotype %s not found in %s\n", genoId, file);
+		return 3;
+	}
+}
+
+int executeTestUsingRandomModel(int shape, TestInvoker &test)
+{
+	Model model;
+	model.open();
+
+	if ((shape < 1) || (shape > 3))
+	{
+		shape = (rand()%3) + 1;
+	}
+
+	Part *part = model.addNewPart(Part::Shape(shape));
+	GeometryTestUtils::randomizePositionScaleAndOrient(part);
+
+	model.close();
+	test(model);
+	GeometryTestUtils::describePart(part, stdout);
+}
+
+class ModelBasedTestInvoker: public TestInvoker
+{
+	private:
+		void (*test)(Model &);
+	public:
+		ModelBasedTestInvoker(void (*_test)(Model &)):
+			test(_test)
+		{}
+		void operator()(Model &model)
+		{
+			test(model);
+		}
+};
+
+int GeometryTestUtils::execute(const SString header, int argc, char *argv[], void (*test)(Model &))
+{
+	srand(time(NULL));
+	
+	if ((argc == 3) && (strcmp("-l", argv[1]) == 0))
+	{
+		return printGenotypesList(argv[2]);
+	}
+	
+	if ((argc == 4) && (strcmp("-l", argv[1]) == 0))
+	{
+		ModelBasedTestInvoker invoker(test);
+		return executeTestUsingLoadedModel(argv[2], argv[3], invoker);
+	}
+	
+	if ((argc == 2) && (strcmp("-c", argv[1]) == 0))
+	{
+		ModelBasedTestInvoker invoker(test);
+		return executeTestUsingRandomModel(-1, invoker);
+	}
+	
+	if ((argc == 3) && (strcmp("-c", argv[1]) == 0) && isdigit(argv[2][0]))
+	{
+		int shape = atol(argv[2]);
+		ModelBasedTestInvoker invoker(test);
+		return executeTestUsingRandomModel(shape, invoker);
+	}
+	
+	fprintf(stderr,
+		"%s\n\n"
+		"argument lists:\n"
+		"-l FILENAME            - to print list of models in file\n"
+		"-l FILENAME GENO_ID    - to load model from file and run test\n"
+		"-c [SHAPE]             - to create simple random model and run test\n\n"
+		"FILENAME - name of file containing named f0 genotypes\n"
+		"GENO_ID - either genotype name or index (1-based)\n"
+		"SHAPE - 1=ellipsoid, 2=cuboid, 3=cylinder, others or none=random\n",
+		(const char*)header);
+	return 1;
+}
+
+class ModelAndDensityBasedTestInvoker: public TestInvoker
+{
+	private:
+		void (*test)(Model &, const double);
+		double density;
+	public:
+		ModelAndDensityBasedTestInvoker(void (*_test)(Model &, const double), double _density):
+			test(_test),
+			density(_density)
+		{}
+		
+		void operator()(Model &model)
+		{
+			test(model, density);
+		}
+};
+
+int GeometryTestUtils::execute(const SString header, int argc, char *argv[],
+	void (*test)(Model &, const double))
+{
+	srand(time(NULL));
+	
+	if ((argc == 3) && (strcmp("-l", argv[1]) == 0))
+	{
+		return printGenotypesList(argv[2]);
+	}
+
+	if ((argc == 5) && (strcmp("-l", argv[1]) == 0) && isdigit(argv[4][0]))
+	{
+		double density = atol(argv[4]);
+		ModelAndDensityBasedTestInvoker invoker(test, density);
+		return executeTestUsingLoadedModel(argv[2], argv[3], invoker);
+	}
+	
+	if ((argc == 3) && (strcmp("-c", argv[1]) == 0) && isdigit(argv[2][0]))
+	{
+		double density = atol(argv[2]);
+		ModelAndDensityBasedTestInvoker invoker(test, density);
+		return executeTestUsingRandomModel(-1, invoker);
+	}
+	
+	if ((argc == 4) && (strcmp("-c", argv[1]) == 0) && isdigit(argv[2][0]) && isdigit(argv[3][0]))
+	{
+		double density = atol(argv[2]);
+		int shape = atol(argv[3]);
+		ModelAndDensityBasedTestInvoker invoker(test, density);
+		return executeTestUsingRandomModel(shape, invoker);
+	}
+	
+	fprintf(stderr,
+		"%s\n\n"
+		"argument lists:\n"
+		"-l FILENAME                    - to print list of models in file\n"
+		"-l FILENAME GENO_ID DENSITY    - to load model from file and run test\n"
+		"-c DENSITY [SHAPE]             - to create simple random model and run test\n\n"
+		"FILENAME - name of file containing named f0 genotypes\n"
+		"GENO_ID - either genotype name or index (1-based)\n"
+		"DENSITY - minimal number of samples per unit\n"
+		"SHAPE - 1=ellipsoid, 2=cuboid, 3=cylinder, others or none=random",
+		(const char*)header);
+	return 1;
+}
+
+void GeometryTestUtils::addAnchorToModel(Model &model)
+{
+	Part *part = model.addNewPart(Part::SHAPE_ELLIPSOID);
+	
+	part->p = Pt3D(0);
+	part->scale = Pt3D(0.1);
+	part->vcolor = Pt3D(1.0, 0.0, 1.0);
+	
+	addAxesToModel(Pt3D(0.5), Orient(Orient_1), Pt3D(0.0), model);
+}
+
+void GeometryTestUtils::addPointToModel(const Pt3D &markerLocation, Model &model)
+{
+	Part *anchor = model.getPart(0);
+	Part *part = model.addNewPart(Part::SHAPE_ELLIPSOID);
+	
+	part->p = Pt3D(markerLocation);
+	part->scale = Pt3D(0.05);
+	part->vcolor = Pt3D(1.0, 1.0, 0.0);
+	
+	model.addNewJoint(anchor, part, Joint::SHAPE_SOLID);
+}
+
+void GeometryTestUtils::addAxesToModel(const Pt3D &sizes, const Orient &axes, const Pt3D &center,
+	Model &model)
+{
+	Part *anchor = model.getPart(0);
+	Part *part;
+	
+	part = model.addNewPart(Part::SHAPE_CUBOID);
+	part->scale = Pt3D(sizes.x, 0.05, 0.05);
+	part->setOrient(axes);
+	part->p = center;
+	part->vcolor = Pt3D(1.0, 0.0, 0.0);
+	model.addNewJoint(anchor, part, Joint::SHAPE_SOLID);
+	
+	part = model.addNewPart(Part::SHAPE_CUBOID);
+	part->scale = Pt3D(0.05, sizes.y, 0.05);
+	part->setOrient(axes);
+	part->p = center;
+	part->vcolor = Pt3D(0.0, 1.0, 0.0);
+	model.addNewJoint(anchor, part, Joint::SHAPE_SOLID);
+	
+	part = model.addNewPart(Part::SHAPE_CUBOID);
+	part->scale = Pt3D(0.05, 0.05, sizes.z);
+	part->setOrient(axes);
+	part->p = center;
+	part->vcolor = Pt3D(0.0, 0.0, 1.0);
+	model.addNewJoint(anchor, part, Joint::SHAPE_SOLID);
+}
+
+void GeometryTestUtils::mergeModels(Model &target, Model &source)
+{
+	Part *targetAnchor = target.getPart(0);
+	Part *sourceAnchor = source.getPart(0);
+	
+	target.moveElementsFrom(source);
+	
+	target.addNewJoint(targetAnchor, sourceAnchor, Joint::SHAPE_SOLID);
+}
+
+double frand(double from, double width)
+{
+	return from + width * ((rand()%10000) / 10000.0);
+}
+
+void GeometryTestUtils::randomizePositionScaleAndOrient(Part *part)
+{
+	part->p = Pt3D(frand(1.5, 1.0), frand(1.5, 1.0), frand(1.5, 1.0));
+	part->scale = Pt3D(frand(0.1, 0.9), frand(0.1, 0.9), frand(0.1, 0.9));
+	part->setRot(Pt3D(frand(0.0, M_PI), frand(0.0, M_PI), frand(0.0, M_PI)));
+}
+
+void GeometryTestUtils::describePart(const Part *part, FILE *output)
+{
+	fprintf(output, "# shape=%d\n", part->shape);
+	fprintf(output, "# x=%f\n", part->p.x);
+	fprintf(output, "# y=%f\n", part->p.y);
+	fprintf(output, "# z=%f\n", part->p.z);
+	fprintf(output, "# sx=%f\n", part->scale.x);
+	fprintf(output, "# sy=%f\n", part->scale.y);
+	fprintf(output, "# sz=%f\n", part->scale.z);
+	fprintf(output, "# rx=%f\n", part->rot.x);
+	fprintf(output, "# ry=%f\n", part->rot.y);
+	fprintf(output, "# rz=%f\n", part->rot.z);
+}
Index: cpp/frams/_demos/geometry/geometrytestutils.h
===================================================================
--- cpp/frams/_demos/geometry/geometrytestutils.h	(revision 191)
+++ cpp/frams/_demos/geometry/geometrytestutils.h	(revision 191)
@@ -0,0 +1,113 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#ifndef _GEOMETRYTESTUTILS_H_
+#define _GEOMETRYTESTUTILS_H_
+
+#include <frams/model/model.h>
+#include <frams/model/modelparts.h>
+#include <frams/util/3d.h>
+#include <frams/util/sstring.h>
+
+/**
+ * @brief Geometry module testing utilities.
+ * @details Contains functions frequently used by tests, such as loading or generating random Model,
+ *     or adding markers of points to a Model.
+ * @author Radosław Gołębiewski
+ */
+namespace GeometryTestUtils
+{
+	/**
+	 * @brief Execution entry point for test based on specified input model.
+	 * @details Parses call arguments and performs one of following tasks: printing usage message;
+	 *     printing list of genotypes in specified file; creating Model from loaded genotype and
+	 *     executing test; creating random Model and executing test.
+	 * @param[in] header Header of usage message. This text should contain test description.
+	 * @param[in] argc Number of call arguments.
+	 * @param[in] argv Call arguments.
+	 * @param[in] test Pointer to function implementing test.
+	 * @returns Result code: 0 - no error, 1 - incorrect arguments (usage message printed),
+	 *     2 - cannot open file, 3 - cannot find genotype, 4 - incorrect genotype.
+	 */
+	int execute(const SString header, int argc, char *argv[], void (*test)(Model &));
+	
+	/**
+	 * @brief Execution entry point for test based on specified input model and density.
+	 * @details Parses call arguments and performs one of following tasks: printing usage message;
+	 *     printing list of genotypes in specified file; creating Model from loaded genotype and
+	 *     executing test using specified density; creating random Model and executing test using
+	 *     specified density.
+	 * @param[in] header Header of usage message. This text should contain test description.
+	 * @param[in] argc Number of call arguments.
+	 * @param[in] argv Call arguments.
+	 * @param[in] test Pointer to function implementing test.
+	 * @returns Result code: 0 - no error, 1 - incorrect arguments (usage message printed),
+	 *     2 - cannot open file, 3 - cannot find genotype, 4 - incorrect genotype.
+	 */
+	int execute(const SString header, int argc, char *argv[], void (*test)(Model &, const double));
+	
+	/**
+	 * @brief Adds anchor to the specified Model.
+	 * @details An anchor has two functions. First is to provide Model consistency. Some functions in
+	 *     GeometryTestUtils namespace requires Model passed to them as an argument to contain at
+	 *     least one Part. All new Parts are bonded to the rest of Model using Joint connecting them
+	 *     with first Part of Model. Second is to provide reference which helps to understand Model
+	 *     position, scale and orientation. Anchor is built from four Parts: small sphere placed in
+	 *     global coordinate system origin and three cuboids visualising global coordinate system
+	 *     axes.
+	 * @see addAxesToModel.
+	 * @param[in] model Owner of Parts to be created.
+	 */
+	void addAnchorToModel(Model &model);
+	
+	/**
+	 * @brief Adds point marker to Model.
+	 * @details Marker of point is a small sphere (radius = 0.05).
+	 * @param[in] point Location of marker.
+	 * @param[in] model Owner of Part to be created, must contain at least one part.
+	 */
+	void addPointToModel(const Pt3D &point, Model &model);
+	
+	/**
+	 * @brief Adds axes markers to Model.
+	 * @details Axes markers are three streched (one of scales = 0.5, others = 0.05) and colored
+	 *     cuboids. Cuboid visualising OX axis is red, OY - green, and OZ - blue.
+	 * @param[in] sizes Axes visual lengths.
+	 * @param[in] axes Axes orientation, relatively to global coordinate system axes.
+	 * @param[in] center Axes intersection point, relatively to global coordinate system origin.
+	 * @param[in] model Owner of Parts to be created, must contain at least one part.
+	 */
+	void addAxesToModel(const Pt3D &sizes, const Orient &axes, const Pt3D &center, Model &model);
+	
+	/**
+	 * @brief Merges two Models.
+	 * @details Moves all parts from source Model to target Model and - to provide Model
+	 *     consistency - creates Joint between firsts Parts of each of them. Each model must contain
+	 *     at least one Part.
+	 * @param[in] target Target Model, must contain at least one part.
+	 * @param[in] source Source Model, must contain at least one part.
+	 */
+	void mergeModels(Model &target, Model &source);
+	
+	/**
+	 * @brief Randomizes position, scale and rotations of Part.
+	 * @details Sets coords of Part position to random values from range (1.5, 2.5), scales to
+	 *     random values from range (0.1, 1.0), and rotations around each axis to random values from
+	 *     range (0, M_PI).
+	 * @param[in] part Part which position, scale and orient should be randomized.
+	 */
+	void randomizePositionScaleAndOrient(Part *part);
+	
+	/**
+	 * @brief Prints description of given Part to specified file.
+	 * @details Prints shape ("shape"), position ("x", "y", "z"), scale ("sx", "sy", "sz") and
+	 *     rotations ("rx", "ry", "rz") of given Part. Each value is printed in separate line using
+	 *     following syntax: "# NAME=VALUE".
+	 * @param[in] part Part to be described.
+	 * @param[in] output Output file.
+	 */
+	void describePart(const Part *part, FILE *output);
+}
+
+#endif
Index: cpp/frams/_demos/geometry/info_test.cpp
===================================================================
--- cpp/frams/_demos/geometry/info_test.cpp	(revision 191)
+++ cpp/frams/_demos/geometry/info_test.cpp	(revision 191)
@@ -0,0 +1,93 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "geometrytestutils.h"
+#include <frams/model/geometry/geometryutils.h>
+#include <frams/model/geometry/modelgeometryinfo.h>
+#include <frams/model/model.h>
+#include <frams/util/3d.h>
+#include <stdio.h>
+
+void test(Model &model, const double density)
+{
+	// EXAMPLES ------------------------------------------------------------------------------------
+	
+	// Calculating approximate volume of model. Parameter density describes minimal number of points
+	// created for each unit of length. It means, that for each cubic unit of bounding box of model,
+	// at least density^3 points will be created.
+	double volume = ModelGeometryInfo::volume(model, density);
+	
+	
+	// Calculating approximate surface area of model. Parameter density describes minimal number of
+	// points created for each unit of length. It means, that for each square unit of surface area
+	// of parts of model, at least density^2 points will be created.
+	double area = ModelGeometryInfo::area(model, density);
+	
+	
+	// Creating output variables.
+	Pt3D sizes;
+	Orient axes;
+	// Calculating approximate sizes and axes of model. These values are calculated together because
+	// they depends on each other. Sizes are distances between furthest points on model surface
+	// measured along corresponding axes. Fields sizes.x and axes.x describes the 'length' of 3D
+	// model. Fields sizes.y and axes.y desctibes 'width' of 2D model created by projecting original
+	// 3D model onto plane defined by axes.x vector. Fields sizes.z and axes.z destribes 'height' of
+	// model. Axes are perpendicular to each other. Parameter density describes minimal number of
+	// points created for each unit of length. It means, that for each square unit of surface area
+	// of parts of model, at least density^2 points will be created.
+	ModelGeometryInfo::findSizesAndAxesOfModel(model, density, sizes, axes);
+	
+	
+	// Creating output variables.
+	Pt3D lowerBoundary, upperBoundary;
+	// Calculationg bounding box of model. Result is stored in variables lowerBoundary and
+	// upperBoundary.
+	ModelGeometryInfo::boundingBox(model, lowerBoundary, upperBoundary);
+	
+	// PRINTING RESULT -----------------------------------------------------------------------------
+	
+	// Creating and preparing result Model object.
+	Model resultModel;
+	resultModel.open();
+	GeometryTestUtils::addAnchorToModel(resultModel);
+	
+	// Adding bounding markers of bounding box apices.
+	for (Octants::Octant o = Octants::FIRST; o < Octants::NUMBER; o = Octants::Octant(o+1))
+	{
+		Pt3D apex;
+		apex.x = Octants::isPositiveX(o) ? upperBoundary.x : lowerBoundary.x;
+		apex.y = Octants::isPositiveY(o) ? upperBoundary.y : lowerBoundary.y;
+		apex.z = Octants::isPositiveZ(o) ? upperBoundary.z : lowerBoundary.z;
+		GeometryTestUtils::addPointToModel(apex, resultModel);
+	}
+	
+	// Adding markers of axes (intersection of axes is in the center of bounding box).
+	Pt3D intersection = (lowerBoundary + upperBoundary) * 0.5;
+	GeometryTestUtils::addAxesToModel(sizes, axes, intersection, resultModel);
+	
+	// Merging with original model.
+	GeometryTestUtils::mergeModels(resultModel, model);
+	
+	// Finishing result Model and printing its genotype.
+	resultModel.close();
+	puts((const char*)resultModel.getF0Geno().toString());
+	
+	// Printing calculated values.
+	printf("# volume=%f\n", volume);
+	printf("# area=%f\n", area);
+	printf("# sizes.x=%f sizes.y=%f sizes.z=%f\n", sizes.x, sizes.y, sizes.z);
+	printf("# axes.x=(%f, %f, %f)\n", axes.x.x, axes.x.y, axes.x.z);
+	printf("# axes.y=(%f, %f, %f)\n", axes.y.x, axes.y.y, axes.y.z);
+	printf("# axes.z=(%f, %f, %f)\n", axes.z.x, axes.z.y, axes.z.z);
+	printf("# box.x=[%f, %f]\n", lowerBoundary.x, upperBoundary.x);
+	printf("# box.y=[%f, %f]\n", lowerBoundary.y, upperBoundary.y);
+	printf("# box.z=[%f, %f]\n", lowerBoundary.z, upperBoundary.z);
+}
+
+int main(int argc, char *argv[])
+{
+	return GeometryTestUtils::execute(
+		"Calculates approximate values of surface area, volume, sizes and axes of model. Adds to "
+		"model markers of bounding box apices and axes.", argc, argv, test);
+}
Index: cpp/frams/_demos/geometry/surface_test.cpp
===================================================================
--- cpp/frams/_demos/geometry/surface_test.cpp	(revision 191)
+++ cpp/frams/_demos/geometry/surface_test.cpp	(revision 191)
@@ -0,0 +1,57 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "geometrytestutils.h"
+#include <frams/model/geometry/meshbuilder.h>
+#include <frams/model/model.h>
+#include <frams/util/3d.h>
+#include <stdio.h>
+
+void test(Model &model, const double density)
+{
+	// Creating and preparing result Model object.
+	Model resultModel;
+	resultModel.open();
+	GeometryTestUtils::addAnchorToModel(resultModel);
+	
+	// Creating instance of Iterator class (MeshBuilder::ModelSurface in this case). Object is
+	// placed in stack memory, thus there is no heap memory allocation and creation is fast.
+	MeshBuilder::ModelSurface iterator(density);
+	
+	// Once created, iterator can be used many times, but must be initialized before every time.
+	iterator.initialize(&model);
+	
+	// Creating iterator output object in stack memory.
+	Pt3D point;
+	
+	// Method tryGetNext checks if next point exists and optionally updates fields of Pt3D object.
+	while (iterator.tryGetNext(point))
+	{
+		// Processing points created by iterator. In this case, they are added to result model as
+		// small spheres.
+		GeometryTestUtils::addPointToModel(point, resultModel);
+	}
+	
+	// After Iterator initialization there are two alternative ways of using it.
+	//
+	// First adds all points to the specified list:
+	// SListTempl<Pt3D> list;
+	// iterator.addAllPointsToList(list);
+	//
+	// Second executes specified operation for each point:
+	// Operation func;
+	// iterator.forEach(func);
+	// In this case, Operation is derived class of MeshBuilder::Callback class.
+	
+	// Finishing result Model and printing its genotype.
+	resultModel.close();
+	puts((const char*)resultModel.getF0Geno().toString());
+}
+
+int main(int argc, char *argv[])
+{
+	return GeometryTestUtils::execute(
+		"Creates new model built of small spheres covering surface of input model and prints its "
+		"genotype.", argc, argv, test);
+}
Index: cpp/frams/_demos/geometry/volume_test.cpp
===================================================================
--- cpp/frams/_demos/geometry/volume_test.cpp	(revision 191)
+++ cpp/frams/_demos/geometry/volume_test.cpp	(revision 191)
@@ -0,0 +1,49 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "geometrytestutils.h"
+#include <frams/model/geometry/geometryutils.h>
+#include <frams/model/geometry/meshbuilder.h>
+#include <frams/model/model.h>
+#include <stdio.h>
+
+void test(Model &model, const double density)
+{
+	// Creating and preparing result Model object.
+	Model resultModel;
+	resultModel.open();
+	GeometryTestUtils::addAnchorToModel(resultModel);
+	
+	// Creating instance of Iterator class (MeshBuilder::ModelApices in this case). Object is
+	// placed in stack memory, thus there is no heap memory allocation and creation is fast.
+	MeshBuilder::BoundingBoxVolume iterator(density);
+	
+	// Once created, iterator can be used many times, but must be initialized before every time.
+	iterator.initialize(model);
+	
+	// Creating output variable.
+	Pt3D point;
+	
+	// Method tryGetNext checks if next point exists and optionally updates fields of Pt3D object.
+	while (iterator.tryGetNext(point))
+	{
+		// Processing points created by iterator. In this case, they are added to result model as
+		// small spheres only if they are placed inside or on surface of input model.
+		if (GeometryUtils::isPointInsideModel(point, model))
+		{
+			GeometryTestUtils::addPointToModel(point, resultModel);
+		}
+	}
+	
+	// Finishing result Model and printing its genotype.
+	resultModel.close();
+	puts((const char*)resultModel.getF0Geno().toString());
+}
+
+int main(int argc, char *argv[])
+{
+	return GeometryTestUtils::execute(
+		"Creates new model built of small spheres filling volume of input model and prints its "
+		"genotype.", argc, argv, test);
+}
Index: cpp/frams/model/geometry/geometryutils.cpp
===================================================================
--- cpp/frams/model/geometry/geometryutils.cpp	(revision 191)
+++ cpp/frams/model/geometry/geometryutils.cpp	(revision 191)
@@ -0,0 +1,351 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "geometryutils.h"
+
+#include <math.h>
+
+double GeometryUtils::pointPosition(const int pointIndex, const int numberOfPoints)
+{
+	return pointIndex / (numberOfPoints-1.0);
+}
+
+double GeometryUtils::pointOnAxis(const double scale, const double position)
+{
+	return (position-0.5) * scale;
+}
+
+double GeometryUtils::pointOnAxis(const double scale, const int pointIndex, const int numberOfPoints)
+{
+	return pointOnAxis(scale, pointPosition(pointIndex, numberOfPoints));
+}
+
+double GeometryUtils::combination(const double value1, const double value2, const double position)
+{
+	return value1 + position * (value2-value1);
+}
+
+double GeometryUtils::combination(const double value1, const double value2, const int pointIndex, const int numberOfPoints)
+{
+	return combination(value1, value2, pointPosition(pointIndex, numberOfPoints));
+}
+
+bool GeometryUtils::isPointInsideModelExcludingPart(const Pt3D &point, const Model *model, const int excludedPartIndex)
+{
+	for (int i = 0; i < excludedPartIndex; i++)
+	{
+		if (isPointInsidePart(point, model->getPart(i)))
+		{
+			return true;
+		}
+	}
+	
+	for (int i = excludedPartIndex+1; i < model->getPartCount(); i++)
+	{
+		if (isPointStrictlyInsidePart(point, model->getPart(i)))
+		{
+			return true;
+		}
+	}
+	
+	return false;
+}
+
+bool GeometryUtils::isPointInsideModel(const Pt3D &point, const Model &model)
+{
+	for (int i = 0; i < model.getPartCount(); i++)
+	{
+		if (isPointInsidePart(point, model.getPart(i)))
+		{
+			return true;
+		}
+	}
+	
+	return false;
+}
+
+bool GeometryUtils::isPointInsidePart(const Pt3D &point, const Part *part)
+{
+	switch (part->shape)
+	{
+		case Part::SHAPE_ELLIPSOID:
+			return isPointInsideEllipsoid(point, part);
+			break;
+			
+		case Part::SHAPE_CUBOID:
+			return isPointInsideCuboid(point, part);
+			break;
+			
+		case Part::SHAPE_CYLINDER:
+			return isPointInsideCylinder(point, part);
+			break;
+	}
+}
+
+bool GeometryUtils::isPointStrictlyInsidePart(const Pt3D &point, const Part *part)
+{
+	switch (part->shape)
+	{
+		case Part::SHAPE_ELLIPSOID:
+			return isPointStrictlyInsideEllipsoid(point, part);
+			break;
+			
+		case Part::SHAPE_CUBOID:
+			return isPointStrictlyInsideCuboid(point, part);
+			break;
+			
+		case Part::SHAPE_CYLINDER:
+			return isPointStrictlyInsideCylinder(point, part);
+			break;
+	}
+}
+
+bool GeometryUtils::isPointInsideEllipsoid(const Pt3D &point, const Part *part)
+{
+	Pt3D moved = point - part->p;
+	Pt3D rotated;
+	part->o.revTransform(rotated, moved);
+	
+	double r
+		= (pow(rotated.x, 2.0) / pow(part->scale.x, 2.0))
+		+ (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
+		+ (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
+	
+	return r <= 1.0;
+}
+
+bool GeometryUtils::isPointStrictlyInsideEllipsoid(const Pt3D &point, const Part *part)
+{
+	Pt3D moved = point - part->p;
+	Pt3D rotated;
+	part->o.revTransform(rotated, moved);
+	
+	double r
+		= (pow(rotated.x, 2.0) / pow(part->scale.x, 2.0))
+		+ (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
+		+ (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
+	
+	return r < 1.0;
+}
+
+bool GeometryUtils::isPointInsideCuboid(const Pt3D &point, const Part *part)
+{
+	Pt3D moved = point - part->p;
+	Pt3D rotated;
+	part->o.revTransform(rotated, moved);
+	
+	return (fabs(rotated.x) <= part->scale.x)
+		&& (fabs(rotated.y) <= part->scale.y)
+		&& (fabs(rotated.z) <= part->scale.z);
+}
+
+bool GeometryUtils::isPointStrictlyInsideCuboid(const Pt3D &point, const Part *part)
+{
+	Pt3D moved = point - part->p;
+	Pt3D rotated;
+	part->o.revTransform(rotated, moved);
+	
+	return (fabs(rotated.x) < part->scale.x)
+		&& (fabs(rotated.y) < part->scale.y)
+		&& (fabs(rotated.z) < part->scale.z);
+}
+
+bool GeometryUtils::isPointInsideCylinder(const Pt3D &point, const Part *part)
+{
+	Pt3D moved = point - part->p;
+	Pt3D rotated;
+	part->o.revTransform(rotated, moved);
+	
+	double r
+		= (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
+		+ (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
+	
+	return (fabs(rotated.x) <= part->scale.x) && (r <= 1.0);
+}
+
+bool GeometryUtils::isPointStrictlyInsideCylinder(const Pt3D &point, const Part *part)
+{
+	Pt3D moved = point - part->p;
+	Pt3D rotated;
+	part->o.revTransform(rotated, moved);
+	
+	double r
+		= (pow(rotated.y, 2.0) / pow(part->scale.y, 2.0))
+		+ (pow(rotated.z, 2.0) / pow(part->scale.z, 2.0));
+	
+	return (fabs(rotated.x) < part->scale.x) && (r < 1.0);
+}
+
+void GeometryUtils::findSizesAndAxesOfPointsGroup(SListTempl<Pt3D> &points, Pt3D &sizes,
+	Orient &axes)
+{
+	findSizeAndAxisOfPointsGroup(points, sizes.x, axes.x);
+	orthographicProjectionToPlane(points, axes.x);
+	findSizeAndAxisOfPointsGroup(points, sizes.y, axes.y);
+	orthographicProjectionToPlane(points, axes.y);
+	
+	Pt3D min(points.get(0)), max(points.get(0));
+	
+	for (int i = 1; i < points.size(); i++)
+	{
+		min.getMin(points.get(i));
+		max.getMax(points.get(i));
+	}
+	
+	sizes.z = min.distanceTo(max);
+	axes.z.vectorProduct(axes.x, axes.y);
+}
+
+void GeometryUtils::findSizeAndAxisOfPointsGroup(const SListTempl<Pt3D> &points, double &size,
+	Pt3D &axis)
+{
+	int index1, index2;
+	size = findTwoFurthestPoints(points, index1, index2);
+	createAxisFromTwoPoints(axis, points.get(index1), points.get(index2));
+}
+
+double GeometryUtils::findTwoFurthestPoints(const SListTempl<Pt3D> &points, int &index1,
+	int &index2)
+{
+	double distance = 0;
+	index1 = index2 = 0;
+	
+	for (int i = 0; i < points.size()-1; i++)
+	{
+		Pt3D p1 = points.get(i);
+		
+		for (int j = i+1; j < points.size(); j++)
+		{
+			Pt3D p2 = points.get(j);
+			double d = p1.distanceTo(p2);
+			
+			if (d > distance)
+			{
+				distance = d;
+				index1 = i;
+				index2 = j;
+			}
+		}
+	}
+	
+	return distance;
+}
+
+void GeometryUtils::createAxisFromTwoPoints(Pt3D &axis, const Pt3D &point1, const Pt3D &point2)
+{
+	Pt3D vector = point2 - point1;
+	vector.normalize();
+	
+	axis.x = vector.x;
+	axis.y = vector.y;
+	axis.z = vector.z;
+}
+
+void GeometryUtils::orthographicProjectionToPlane(SListTempl<Pt3D> &points,
+	const Pt3D &planeNormalVector)
+{
+	for (int i = 0; i < points.size(); i++)
+	{
+		Pt3D &point = points.get(i);
+		
+		double distance = pointDistanceToPlane(point, planeNormalVector);
+		
+		point.x -= planeNormalVector.x * distance;
+		point.y -= planeNormalVector.y * distance;
+		point.z -= planeNormalVector.z * distance;
+	}
+}
+
+double GeometryUtils::pointDistanceToPlane(const Pt3D &point, const Pt3D &planeNormalVector)
+{
+	return planeNormalVector.x*point.x + planeNormalVector.y*point.y + planeNormalVector.z*point.z;
+}
+
+void GeometryUtils::getRectangleApicesFromCuboid(const Part *part, const CuboidFaces::Face face, Pt3D &apex1, Pt3D &apex2, Pt3D &apex3, Pt3D &apex4)
+{
+	Pt3D temp1(part->scale), temp2(part->scale), temp3(part->scale), temp4(part->scale);
+	
+	if (CuboidFaces::isX(face))
+	{
+		temp2.z *= -1;
+		temp3.y *= -1;
+		temp4.z *= -1;
+		temp4.y *= -1;
+	}
+	else if (CuboidFaces::isY(face))
+	{
+		temp2.x *= -1;
+		temp3.z *= -1;
+		temp4.x *= -1;
+		temp4.z *= -1;
+	}
+	else if (CuboidFaces::isZ(face))
+	{
+		temp2.y *= -1;
+		temp3.x *= -1;
+		temp4.y *= -1;
+		temp4.x *= -1;
+	}
+	
+	if (CuboidFaces::isNegative(face))
+	{
+		temp1 *= -1;
+		temp2 *= -1;
+		temp3 *= -1;
+		temp4 *= -1;
+	}
+
+	part->o.transform(apex1, temp1);
+	part->o.transform(apex2, temp2);
+	part->o.transform(apex3, temp3);
+	part->o.transform(apex4, temp4);
+
+	apex1 += part->p;
+	apex2 += part->p;
+	apex3 += part->p;
+	apex4 += part->p;
+}
+
+void GeometryUtils::getRectangleApices(const double width, const double height, const Pt3D &position, const Orient &orient, Pt3D &apex1, Pt3D &apex2, Pt3D &apex3, Pt3D &apex4)
+{
+	Pt3D temp1(0.0, +width, +height);
+	Pt3D temp2(0.0, +width, -height);
+	Pt3D temp3(0.0, -width, +height);
+	Pt3D temp4(0.0, -width, -height);
+
+	orient.transform(apex1, temp1);
+	orient.transform(apex2, temp2);
+	orient.transform(apex3, temp3);
+	orient.transform(apex4, temp4);
+
+	apex1 += position;
+	apex2 += position;
+	apex3 += position;
+	apex4 += position;
+}
+
+void GeometryUtils::getNextEllipseSegmentationPoint(const double d, const double a, const double b, double &x, double &y)
+{
+	x += d / sqrt(1.0 + (b*b * x*x) / (a*a * (a*a - x*x)));
+	y = b * sqrt(1.0 - (x*x) / (a*a));
+}
+
+double GeometryUtils::ellipsoidArea(const Pt3D &sizes)
+{
+	return ellipsoidArea(sizes.x, sizes.y, sizes.z);
+}
+
+double GeometryUtils::ellipsoidArea(const double a, const double b, const double c)
+{
+	double p = 1.6075;
+	double ap = pow(a, p);
+	double bp = pow(b, p);
+	double cp = pow(c, p);
+	return 4*M_PI * pow((ap*bp + bp*cp + cp*ap) / 3.0, 1.0 / p);
+}
+
+double GeometryUtils::ellipsePerimeter(const double a, const double b)
+{
+	return M_PI * ((3 * (a+b)) - sqrt((3*a + b) * (a + 3*b)));
+}
Index: cpp/frams/model/geometry/geometryutils.h
===================================================================
--- cpp/frams/model/geometry/geometryutils.h	(revision 191)
+++ cpp/frams/model/geometry/geometryutils.h	(revision 191)
@@ -0,0 +1,157 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#ifndef _GEOMETRYUTILS_H_
+#define _GEOMETRYUTILS_H_
+
+#include <frams/model/model.h>
+#include <frams/model/modelparts.h>
+#include <frams/util/3d.h>
+#include <frams/util/list.h>
+
+namespace CuboidFaces
+{
+	enum Face
+	{
+		NEGATIVE_X = 0,
+		POSITIVE_X = 1,
+		NEGATIVE_Y = 2,
+		POSITIVE_Y = 3,
+		NEGATIVE_Z = 4,
+		POSITIVE_Z = 5,
+		FIRST = 0,
+		NUMBER = 6
+	};
+	
+	inline bool isPositive(Face f) { return f & 0b001; }
+	inline bool isNegative(Face f) { return !isPositive(f); }
+	inline bool isX(Face f) { return (f & 0b110) == 0b000; }
+	inline bool isY(Face f) { return (f & 0b110) == 0b010; }
+	inline bool isZ(Face f) { return (f & 0b110) == 0b100; }
+}
+
+namespace CylinderBases
+{
+	enum Base
+	{
+		NEGATIVE_X = 0,
+		POSITIVE_X = 1,
+		FIRST = 0,
+		NUMBER = 2
+	};
+	
+	inline bool isPositive(Base b) { return b & 0b001; }
+	inline bool isNegative(Base b) { return !isPositive(b); }
+}
+
+namespace QuadrantsXY
+{
+	enum QuadrantXY
+	{
+		NEGATIVE_X_NEGATIVE_Y = 0,
+		NEGATIVE_X_POSITIVE_Y = 1,
+		POSITIVE_X_NEGATIVE_Y = 2,
+		POSITIVE_X_POSITIVE_Y = 3,
+		FIRST = 0,
+		NUMBER = 4
+	};
+	
+	inline bool isPositiveX(QuadrantXY q) { return q & 0b10; }
+	inline bool isNegativeX(QuadrantXY q) { return !isPositiveX(q); }
+	inline bool isPositiveY(QuadrantXY q) { return q & 0b01; }
+	inline bool isNegativeY(QuadrantXY q) { return !isPositiveY(q); }
+}
+
+namespace QuadrantsYZ
+{
+	enum QuadrantYZ
+	{
+		NEGATIVE_Y_NEGATIVE_Z = 0,
+		NEGATIVE_Y_POSITIVE_Z = 1,
+		POSITIVE_Y_NEGATIVE_Z = 2,
+		POSITIVE_Y_POSITIVE_Z = 3,
+		FIRST = 0,
+		NUMBER = 4
+	};
+	
+	inline bool isPositiveY(QuadrantYZ q) { return q & 0b10; }
+	inline bool isNegativeY(QuadrantYZ q) { return !isPositiveY(q); }
+	inline bool isPositiveZ(QuadrantYZ q) { return q & 0b01; }
+	inline bool isNegativeZ(QuadrantYZ q) { return !isPositiveZ(q); }
+}
+
+namespace QuadrantsZX
+{
+	enum QuadrantZX
+	{
+		NEGATIVE_Z_NEGATIVE_X = 0,
+		NEGATIVE_Z_POSITIVE_X = 1,
+		POSITIVE_Z_NEGATIVE_X = 2,
+		POSITIVE_Z_POSITIVE_X = 3,
+		FIRST = 0,
+		NUMBER = 4
+	};
+	
+	inline bool isPositiveZ(QuadrantZX q) { return q & 0b10; }
+	inline bool isNegativeZ(QuadrantZX q) { return !isPositiveZ(q); }
+	inline bool isPositiveX(QuadrantZX q) { return q & 0b01; }
+	inline bool isNegativeX(QuadrantZX q) { return !isPositiveX(q); }
+}
+
+namespace Octants
+{
+	enum Octant
+	{
+		NEGATIVE_X_NEGATIVE_Y_NEGATIVE_Z = 0,
+		NEGATIVE_X_NEGATIVE_Y_POSITIVE_Z = 1,
+		NEGATIVE_X_POSITIVE_Y_NEGATIVE_Z = 2,
+		NEGATIVE_X_POSITIVE_Y_POSITIVE_Z = 3,
+		POSITIVE_X_NEGATIVE_Y_NEGATIVE_Z = 4,
+		POSITIVE_X_NEGATIVE_Y_POSITIVE_Z = 5,
+		POSITIVE_X_POSITIVE_Y_NEGATIVE_Z = 6,
+		POSITIVE_X_POSITIVE_Y_POSITIVE_Z = 7,
+		FIRST = 0,
+		NUMBER = 8
+	};
+	
+	inline bool isPositiveX(Octant o) { return o & 0b100; }
+	inline bool isNegativeX(Octant o) { return !isPositiveX(o); }
+	inline bool isPositiveY(Octant o) { return o & 0b010; }
+	inline bool isNegativeY(Octant o) { return !isPositiveY(o); }
+	inline bool isPositiveZ(Octant o) { return o & 0b001; }
+	inline bool isNegativeZ(Octant o) { return !isPositiveZ(o); }
+}
+
+namespace GeometryUtils
+{
+	double pointPosition(const int pointIndex, const int numberOfPoints);
+	double pointOnAxis(const double scale, const double position);
+	double pointOnAxis(const double scale, const int pointIndex, const int numberOfPoints);
+	double combination(const double value1, const double value2, const double position);
+	double combination(const double value1, const double value2, const int pointIndex, const int numberOfPoints);
+	bool isPointInsideModelExcludingPart(const Pt3D &point, const Model *model, const int excludedPartIndex);
+	bool isPointInsideModel(const Pt3D &point, const Model &model);
+	bool isPointInsidePart(const Pt3D &point, const Part *part);
+	bool isPointStrictlyInsidePart(const Pt3D &point, const Part *part);
+	bool isPointInsideEllipsoid(const Pt3D &point, const Part *part);
+	bool isPointStrictlyInsideEllipsoid(const Pt3D &point, const Part *part);
+	bool isPointInsideCuboid(const Pt3D &point, const Part *part);
+	bool isPointStrictlyInsideCuboid(const Pt3D &point, const Part *part);
+	bool isPointInsideCylinder(const Pt3D &point, const Part *part);
+	bool isPointStrictlyInsideCylinder(const Pt3D &point, const Part *part);
+	void findSizesAndAxesOfPointsGroup(SListTempl<Pt3D> &points, Pt3D &sizes, Orient &axes);
+	void findSizeAndAxisOfPointsGroup(const SListTempl<Pt3D> &points, double &size, Pt3D &axis);
+	double findTwoFurthestPoints(const SListTempl<Pt3D> &points, int &index1, int &index2);
+	void createAxisFromTwoPoints(Pt3D &axis, const Pt3D &point1, const Pt3D &point2);
+	void orthographicProjectionToPlane(SListTempl<Pt3D> &points, const Pt3D &planeNormalVector);
+	double pointDistanceToPlane(const Pt3D &point, const Pt3D &planeNormalVector);
+	void getRectangleApicesFromCuboid(const Part *part, const CuboidFaces::Face face, Pt3D &apex1, Pt3D &apex2, Pt3D &apex3, Pt3D &apex4);
+	void getRectangleApices(const double width, const double height, const Pt3D &position, const Orient &orient, Pt3D &apex1, Pt3D &apex2, Pt3D &apex3, Pt3D &apex4);
+	void getNextEllipseSegmentationPoint(const double d, const double a, const double b, double &x, double &y);
+	double ellipsoidArea(const Pt3D &sizes);
+	double ellipsoidArea(const double a, const double b, const double c);
+	double ellipsePerimeter(const double a, const double b);
+}
+
+#endif
Index: cpp/frams/model/geometry/meshbuilder.cpp
===================================================================
--- cpp/frams/model/geometry/meshbuilder.cpp	(revision 191)
+++ cpp/frams/model/geometry/meshbuilder.cpp	(revision 191)
@@ -0,0 +1,826 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "meshbuilder.h"
+
+#include <frams/model/geometry/modelgeometryinfo.h>
+#include <stdio.h>
+#include <math.h>
+
+MeshBuilder::Iterator::Iterator(const double _density):
+	density(_density)
+{}
+
+double MeshBuilder::Iterator::getDensity() const
+{
+	return density;
+}
+
+double MeshBuilder::Iterator::setDensity(const double _density)
+{
+	density = _density;
+}
+
+void MeshBuilder::Iterator::forEach(Callback &callback)
+{
+	Pt3D point;
+	
+	while (tryGetNext(point))
+	{
+		callback(point);
+	}
+}
+
+void MeshBuilder::Iterator::addAllPointsToList(SListTempl<Pt3D> &list)
+{
+	Pt3D point;
+	
+	while (tryGetNext(point))
+	{
+		list += point;
+	}
+}
+
+
+
+MeshBuilder::Segment::Segment(const double _density):
+	Iterator(_density), point1(0), point2(0), i(0), I(0), numberOfPoints(0)
+{}
+
+void MeshBuilder::Segment::initialize(const Pt3D &point1, const Pt3D &point2, const bool skipFirst, const bool skipLast, const bool forceOddNumberOfPoints)
+{
+	this->point1 = point1;
+	this->point2 = point2;
+	numberOfPoints = (int)ceil(density * point1.distanceTo(point2)) + 1;
+	numberOfPoints += forceOddNumberOfPoints && (numberOfPoints%2 == 0) ? 1 : 0;
+	i = skipFirst ? 1 : 0;
+	I = numberOfPoints - (skipLast ? 1 : 0);
+}
+
+bool MeshBuilder::Segment::tryGetNext(Pt3D &point)
+{
+	if (i >= I)
+	{
+		return false;
+	}
+	
+	double position = GeometryUtils::pointPosition(i, numberOfPoints);
+	point.x = GeometryUtils::combination(point1.x, point2.x, position);
+	point.y = GeometryUtils::combination(point1.y, point2.y, position);
+	point.z = GeometryUtils::combination(point1.z, point2.z, position);
+	
+	i += 1;
+	return true;
+}
+
+
+
+MeshBuilder::RectangleSurface::RectangleSurface(const double _density):
+	Iterator(_density), edge1(_density), edge2(_density), area(_density), skipVerticalEdges(false), forceOddNumberOfPoints(false)
+{}
+
+void MeshBuilder::RectangleSurface::initialize(const Part *part, const CuboidFaces::Face face, const bool skipVerticalEdges, const bool skipHorizontalEdges, const bool forceOddNumberOfPoints)
+{
+	Pt3D apex1, apex2, apex3, apex4;
+	GeometryUtils::getRectangleApicesFromCuboid(part, face, apex1, apex2, apex3, apex4);
+	initialize(apex1, apex2, apex3, apex4, skipVerticalEdges, skipHorizontalEdges, forceOddNumberOfPoints);
+}
+
+void MeshBuilder::RectangleSurface::initialize(const double width, const double height, const Pt3D &position, const Orient &orient, const bool skipVerticalEdges, const bool skipHorizontalEdges, const bool forceOddNumberOfPoints)
+{
+	Pt3D apex1, apex2, apex3, apex4;
+	GeometryUtils::getRectangleApices(width, height, position, orient, apex1, apex2, apex3, apex4);
+	initialize(apex1, apex2, apex3, apex4, skipVerticalEdges, skipHorizontalEdges, forceOddNumberOfPoints);
+}
+
+void MeshBuilder::RectangleSurface::initialize(const Pt3D &apex1, const Pt3D &apex2, const Pt3D &apex3, const Pt3D &apex4, const bool skipVerticalEdges, const bool skipHorizontalEdges, const bool forceOddNumberOfPoints)
+{
+	this->skipVerticalEdges = skipVerticalEdges;
+	this->forceOddNumberOfPoints = forceOddNumberOfPoints;
+	
+	edge1.setDensity(density);
+	edge1.initialize(apex1, apex2, skipHorizontalEdges, skipHorizontalEdges, forceOddNumberOfPoints);
+	
+	edge2.setDensity(density);
+	edge2.initialize(apex3, apex4, skipHorizontalEdges, skipHorizontalEdges, forceOddNumberOfPoints);
+	
+	area = Segment(density);
+}
+
+bool MeshBuilder::RectangleSurface::tryGetNext(Pt3D &point)
+{
+	while (!area.tryGetNext(point))
+	{
+		Pt3D point1, point2;
+		if (edge1.tryGetNext(point1) && edge2.tryGetNext(point2))
+		{
+			area.initialize(point1, point2, skipVerticalEdges, skipVerticalEdges, forceOddNumberOfPoints);
+		}
+		else
+		{
+			return false;
+		}
+	}
+	
+	return true;
+}
+
+
+
+MeshBuilder::CuboidApices::CuboidApices():
+	Iterator(0), scale(0), position(0), orient(Orient_1), octant(Octants::NUMBER)
+{}
+
+void MeshBuilder::CuboidApices::initialize(const Part *part)
+{
+	initialize(part->scale, part->p, part->o);
+}
+
+void MeshBuilder::CuboidApices::initialize(const Pt3D &scale, const Pt3D &position, const Orient &orient)
+{
+	this->scale = scale;
+	this->position = position;
+	this->orient = orient;
+	octant = Octants::FIRST;
+}
+
+bool MeshBuilder::CuboidApices::tryGetNext(Pt3D &point)
+{
+	if (octant == Octants::NUMBER)
+	{
+		return false;
+	}
+	
+	Pt3D temp(scale);
+	
+	temp.x *= Octants::isPositiveX(octant) ? 1 : -1;
+	temp.y *= Octants::isPositiveY(octant) ? 1 : -1;
+	temp.z *= Octants::isPositiveZ(octant) ? 1 : -1;
+	
+	orient.transform(point, temp);
+	point += position;
+	
+	octant = Octants::Octant(octant+1);
+	
+	return true;
+}
+
+
+
+MeshBuilder::CuboidSurface::CuboidSurface(const double _density):
+	Iterator(_density), apices(), rectangle(_density), iterator(NULL), part(NULL), face(CuboidFaces::NUMBER)
+{}
+
+void MeshBuilder::CuboidSurface::initialize(const Part *part)
+{
+	this->part = part;
+	face = CuboidFaces::FIRST;
+	apices.setDensity(density);
+	apices.initialize(part);
+	iterator = &apices;
+	rectangle.setDensity(density);
+}
+
+bool MeshBuilder::CuboidSurface::tryGetNext(Pt3D &point)
+{
+	if (iterator == NULL)
+	{
+		return false;
+	}
+	
+	while(!iterator->tryGetNext(point))
+	{
+		if (face < CuboidFaces::NUMBER)
+		{
+			iterator = &rectangle;
+			rectangle.initialize(part, face, true, false);
+			face = CuboidFaces::Face(face+1);
+		}
+		else
+		{
+			return false;
+		}
+	}
+	
+	return true;
+}
+
+
+
+MeshBuilder::EllipseSurface::EllipseSurface(const double _density):
+	Iterator(_density), rectangle(_density), position(0), orient(Orient_1), inversedSquaredWidth(0.0), inversedSquaredHeight(0.0)
+{}
+
+void MeshBuilder::EllipseSurface::initialize(const Part *part, const CylinderBases::Base base)
+{
+	double relativePositionX = part->scale.x;
+	relativePositionX *= CylinderBases::isPositive(base) ? +1 : -1;
+	Pt3D relativePosition;
+	part->o.transform(relativePosition, Pt3D(relativePositionX, 0.0, 0.0));
+	initialize(part->scale.y, part->scale.z, part->p + relativePosition, part->o);
+}
+
+void MeshBuilder::EllipseSurface::initialize(const double width, const double height, const Pt3D &position, const Orient &orient)
+{
+	this->position = position;
+	this->orient = orient;
+	rectangle.setDensity(density);
+	rectangle.initialize(width, height, Pt3D(0), Orient_1, false, false, true);
+	inversedSquaredWidth = 1.0 / pow(width, 2.0);
+	inversedSquaredHeight = 1.0 / pow(height, 2.0);
+}
+
+bool MeshBuilder::EllipseSurface::tryGetNext(Pt3D &point)
+{
+	Pt3D temp;
+	
+	while (rectangle.tryGetNext(temp))
+	{
+		double r
+			= (pow(temp.y, 2.0) * inversedSquaredWidth)
+			+ (pow(temp.z, 2.0) * inversedSquaredHeight);
+		if (r <= 1.0)
+		{
+			orient.transform(point, temp);
+			point += position;
+			return true;
+		}
+	}
+	
+	return false;
+}
+
+
+
+MeshBuilder::Ellipse::Ellipse(const double _density):
+	Iterator(_density), phase(Done), position(0), orient(Orient_1), width(0.0), height(0.0), d(0.0), p(0.0), q(0.0), P(0.0), Q(0.0), D(0.0), a(0.0), b(0.0), quadrant(QuadrantsYZ::FIRST)
+{}
+
+void MeshBuilder::Ellipse::initialize(const Part *part, const CylinderBases::Base base)
+{
+	double relativePositionX = part->scale.x;
+	relativePositionX *= CylinderBases::isPositive(base) ? +1 : -1;
+	Pt3D relativePosition;
+	part->o.transform(relativePosition, Pt3D(relativePositionX, 0.0, 0.0));
+	initialize(part->scale.y, part->scale.z, part->p + relativePosition, part->o);
+}
+
+void MeshBuilder::Ellipse::initialize(const double width, const double height, const Pt3D &position, const Orient &orient)
+{
+	this->width = width;
+	this->height = height;
+	this->position = position;
+	this->orient = orient;
+	calculateAndSetSegmentationDistance();
+	initializePhase(yToZ);
+	D = 0;
+}
+
+bool MeshBuilder::Ellipse::tryGetNext(Pt3D &point)
+{
+	if (phase == Done)
+	{
+		return false;
+	}
+	
+	setPoint(point);
+	
+	quadrant = QuadrantsYZ::QuadrantYZ(quadrant+1);
+	
+	if (quadrant == QuadrantsYZ::NUMBER)
+	{
+		quadrant = QuadrantsYZ::FIRST;
+		findNextPQAndPhase();
+	}
+	
+	return true;
+}
+
+void MeshBuilder::Ellipse::calculateAndSetSegmentationDistance()
+{
+	d = 1.0 / density;
+	int i = 0;
+	
+	Pt3D p1(0);
+	i += findLastPQOfPhase(yToZ, p1.z, p1.y);
+	
+	Pt3D p2(0);
+	i += findLastPQOfPhase(zToY, p2.y, p2.z);
+		
+	double ld = p1.distanceTo(p2);
+	double td = i * d + ld;
+	int n = (int)ceil(td * density);
+	d = td / n;
+}
+
+void MeshBuilder::Ellipse::initializePhase(Phase ph)
+{
+	phase = ph;
+	
+	if (phase != Done)
+	{
+		a = phase == yToZ ? height : width;
+		b = phase == yToZ ? width : height;
+		
+		p = 0.0;
+		q = b;
+		P = a*a / sqrt(a*a + b*b);
+		Q = b*b / sqrt(a*a + b*b);
+		
+		quadrant = QuadrantsYZ::POSITIVE_Y_NEGATIVE_Z;
+	}
+}
+
+void MeshBuilder::Ellipse::setPoint(Pt3D &point)
+{
+	double np = p * (QuadrantsYZ::isPositiveY(quadrant) ? +1 : -1);
+	double nq = q * (QuadrantsYZ::isPositiveZ(quadrant) ? +1 : -1);
+	
+	double y = phase == yToZ ? nq : np;
+	double z = phase == yToZ ? np : nq;
+	
+	Pt3D temp(0.0, y, z);
+	orient.transform(point, temp);
+	point += position;
+}
+
+void MeshBuilder::Ellipse::findNextPQAndPhase()
+{
+	double _p = p;
+	double _q = q;
+	
+	GeometryUtils::getNextEllipseSegmentationPoint(d, a, b, p, q);
+	
+	if (p > P)
+	{
+		D += sqrt(pow(P-_p, 2.0) + pow(Q-_q, 2.0));
+		
+		if ((phase == zToY) && (D > 1.5 * d))
+		{
+			p = P;
+			q = Q;
+			D = 0;
+		}
+		else
+		{
+			initializePhase(Phase(phase+1));
+		}
+	}
+}
+
+int MeshBuilder::Ellipse::findLastPQOfPhase(Phase ph, double &lp, double &lq)
+{
+	initializePhase(ph);
+	int i = 0;
+	
+	do
+	{
+		lp = p;
+		lq = q;
+		
+		GeometryUtils::getNextEllipseSegmentationPoint(d, a, b, p, q);
+		
+		i += p <= P ? 1 : 0;
+	} while (p <= P);
+	
+	return i;
+}
+
+MeshBuilder::CylinderEdges::CylinderEdges(const double _density):
+	Iterator(_density), ellipse(_density), edge(0.0), length(0.0), base(CylinderBases::NUMBER)
+{}
+
+void MeshBuilder::CylinderEdges::initialize(const Part *part)
+{
+	initialize(part->scale, part->p, part->o);
+}
+
+void MeshBuilder::CylinderEdges::initialize(const Pt3D &scale, const Pt3D &position, const Orient &orient)
+{
+	initialize(scale.x, scale.y, scale.z, position, orient);
+}
+
+void MeshBuilder::CylinderEdges::initialize(const double length, const double width, const double height, const Pt3D &position, const Orient &orient)
+{
+	ellipse.setDensity(density);
+	ellipse.initialize(width, height, position, orient);
+	
+	orient.transform(this->length, Pt3D(length, 0.0, 0.0));
+	
+	base = CylinderBases::NUMBER;
+}
+
+bool MeshBuilder::CylinderEdges::tryGetNext(Pt3D &point)
+{
+	if (base == CylinderBases::NUMBER)
+	{
+		if (!ellipse.tryGetNext(edge))
+		{
+			return false;
+		}
+		
+		base = CylinderBases::FIRST;
+	}
+	
+	point.x = edge.x + length.x * (CylinderBases::isPositive(base) ? +1 : -1);
+	point.y = edge.y + length.y * (CylinderBases::isPositive(base) ? +1 : -1);
+	point.z = edge.z + length.z * (CylinderBases::isPositive(base) ? +1 : -1);
+	
+	base = CylinderBases::Base(base+1);
+	
+	return true;
+}
+
+
+
+MeshBuilder::CylinderWallSurface::CylinderWallSurface(const double _density):
+	Iterator(_density), edge(_density), length(0.0), area(_density)
+{}
+
+void MeshBuilder::CylinderWallSurface::initialize(const Part *part)
+{
+	initialize(part->scale, part->p, part->o);
+}
+
+void MeshBuilder::CylinderWallSurface::initialize(const Pt3D &scale, const Pt3D &position, const Orient &orient)
+{
+	initialize(scale.x, scale.y, scale.z, position, orient);
+}
+
+void MeshBuilder::CylinderWallSurface::initialize(const double length, const double width, const double height, const Pt3D &position, const Orient &orient)
+{
+	edge.setDensity(density);
+	edge.initialize(width, height, position, orient);
+	
+	orient.transform(this->length, Pt3D(length, 0.0, 0.0));
+	
+	area.setDensity(density);
+}
+
+bool MeshBuilder::CylinderWallSurface::tryGetNext(Pt3D &point)
+{
+	while (!area.tryGetNext(point))
+	{
+		if (edge.tryGetNext(point))
+		{
+			area.initialize(point + length, point - length);
+		}
+		else
+		{
+			return false;
+		}
+	}
+	
+	return true;
+}
+
+
+
+MeshBuilder::CylinderSurface::CylinderSurface(const double _density):
+	Iterator(_density), wall(_density), ellipse(_density), iterator(NULL), part(NULL), base(CylinderBases::NUMBER)
+{}
+
+void MeshBuilder::CylinderSurface::initialize(const Part *part)
+{
+	this->part = part;
+	base = CylinderBases::FIRST;
+	wall.setDensity(density);
+	wall.initialize(part);
+	iterator = &wall;
+	ellipse.setDensity(density);
+}
+
+bool MeshBuilder::CylinderSurface::tryGetNext(Pt3D &point)
+{
+	if (iterator == NULL)
+	{
+		return false;
+	}
+	
+	while(!iterator->tryGetNext(point))
+	{
+		if (base < CylinderBases::NUMBER)
+		{
+			iterator = &ellipse;
+			ellipse.initialize(part, base);
+			base = CylinderBases::Base(base+1);
+		}
+		else
+		{
+			return false;
+		}
+	}
+	
+	return true;
+}
+
+
+
+MeshBuilder::EllipsoidSurface::EllipsoidSurface(const double _density):
+	Iterator(_density), phase(Done), edge(0.0), area(0.0), scale(0.0), limit(0.0), d(0.0), part(NULL), octant(Octants::NUMBER)
+{}
+
+void MeshBuilder::EllipsoidSurface::initialize(const Part *part)
+{
+	this->part = part;
+	d = 1.0 / density;
+	initializePhase(X);
+}
+
+bool MeshBuilder::EllipsoidSurface::tryGetNext(Pt3D &point)
+{
+	if (phase == Done)
+	{
+		return false;
+	}
+	
+	setPoint(point);
+	proceedToNextOctant();
+	
+	if (octant == Octants::NUMBER)
+	{
+		octant = Octants::FIRST;
+		findNextAreaEdgeAndPhase();
+	}
+	
+	return true;
+}
+
+void MeshBuilder::EllipsoidSurface::initializePhase(Phase ph)
+{
+	phase = ph;
+	
+	if (phase != Done)
+	{
+		scale.x = phase == X ? part->scale.x : phase == Y ? part->scale.y : part->scale.z;
+		scale.y = phase == X ? part->scale.y : phase == Y ? part->scale.z : part->scale.x;
+		scale.z = phase == X ? part->scale.z : phase == Y ? part->scale.x : part->scale.y;
+		
+		edge = Pt3D(scale.x, 0.0, scale.z);
+		limit.y = scale.y*scale.y / sqrt(scale.x*scale.x + scale.y*scale.y);
+		limit.z = edge.z*edge.z / sqrt(edge.x*edge.x + edge.z*edge.z);
+		area = Pt3D(edge.x, edge.y, 0.0);
+		
+		octant = Octants::FIRST;
+	}
+}
+
+void MeshBuilder::EllipsoidSurface::setPoint(Pt3D &point)
+{
+	Pt3D temp1(area);
+	temp1.x *= Octants::isPositiveX(octant) ? +1 : -1;
+	temp1.y *= Octants::isPositiveY(octant) ? +1 : -1;
+	temp1.z *= Octants::isPositiveZ(octant) ? +1 : -1;
+	
+	Pt3D temp2;
+	temp2.x = phase == X ? temp1.x : phase == Y ? temp1.z : temp1.y;
+	temp2.y = phase == X ? temp1.y : phase == Y ? temp1.x : temp1.z;
+	temp2.z = phase == X ? temp1.z : phase == Y ? temp1.y : temp1.x;
+	
+	part->o.transform(point, temp2);
+	point += part->p;
+}
+
+void MeshBuilder::EllipsoidSurface::proceedToNextOctant()
+{
+	int step = 1;
+	step += (Octants::isNegativeZ(octant) && (area.z == 0.0)) ? 1 : 0;
+	step += (Octants::isNegativeY(octant) && (area.y == 0.0)) ? 2 : 0;
+	octant = Octants::Octant(octant + step);
+}
+
+void MeshBuilder::EllipsoidSurface::findNextAreaEdgeAndPhase()
+{
+	GeometryUtils::getNextEllipseSegmentationPoint(d, edge.z, edge.x, area.z, area.x);
+		
+	if ((area.z > limit.z) || (area.y > limit.y * sqrt(1.0 - (area.z*area.z) / (scale.z*scale.z))))
+	{
+		GeometryUtils::getNextEllipseSegmentationPoint(d, scale.y, scale.x, edge.y, edge.x);
+		edge.z = scale.z * sqrt(1.0 - (edge.y*edge.y) / (scale.y*scale.y));
+		limit.z = edge.z*edge.z / sqrt(edge.x*edge.x + edge.z*edge.z);
+		area = Pt3D(edge.x, edge.y, 0.0);
+		
+		if (edge.y > limit.y)
+		{
+			initializePhase(Phase(phase+1));
+		}
+	}
+}
+
+
+
+MeshBuilder::PartSurface::PartSurface(const double _density):
+	Iterator(_density), cuboid(_density), cylinder(_density), ellipsoid(_density), iterator(NULL)
+{}
+
+void MeshBuilder::PartSurface::initialize(const Part *part)
+{
+	switch (part->shape)
+	{
+		case Part::SHAPE_ELLIPSOID:
+			ellipsoid.setDensity(density);
+			ellipsoid.initialize(part);
+			iterator = &ellipsoid;
+			break;
+			
+		case Part::SHAPE_CUBOID:
+			cuboid.setDensity(density);
+			cuboid.initialize(part);
+			iterator = &cuboid;
+			break;
+			
+		case Part::SHAPE_CYLINDER:
+			cylinder.setDensity(density);
+			cylinder.initialize(part);
+			iterator = &cylinder;
+			break;
+	}
+}
+
+bool MeshBuilder::PartSurface::tryGetNext(Pt3D &point)
+{
+	if (iterator == NULL)
+	{
+		return false;
+	}
+
+	return iterator->tryGetNext(point);
+}
+
+
+
+MeshBuilder::ModelSurface::ModelSurface(const double _density):
+	Iterator(_density), surface(_density), model(NULL), index(0)
+{}
+
+void MeshBuilder::ModelSurface::initialize(const Model *model)
+{
+	this->model = model;
+	surface.setDensity(density);
+	
+	index = 0;
+	if (model->getPartCount() > index)
+	{
+		surface.initialize(model->getPart(index));
+	}
+}
+
+bool MeshBuilder::ModelSurface::tryGetNext(Pt3D &point)
+{
+	if (model == NULL)
+	{
+		return false;
+	}
+	
+	do
+	{
+		while (!surface.tryGetNext(point))
+		{
+			index += 1;
+			if (model->getPartCount() > index)
+			{
+				surface.initialize(model->getPart(index));
+			}
+			else
+			{
+				return false;
+			}
+		}
+	} while (GeometryUtils::isPointInsideModelExcludingPart(point, model, index));
+	
+	return true;
+}
+
+
+
+MeshBuilder::PartApices::PartApices(const double _density):
+	Iterator(_density), cuboid(), cylinder(_density), ellipsoid(_density), iterator(NULL)
+{}
+
+void MeshBuilder::PartApices::initialize(const Part *part)
+{
+	switch (part->shape)
+	{
+		case Part::SHAPE_ELLIPSOID:
+			ellipsoid.setDensity(density);
+			ellipsoid.initialize(part);
+			iterator = &ellipsoid;
+			break;
+			
+		case Part::SHAPE_CUBOID:
+			cuboid.initialize(part);
+			iterator = &cuboid;
+			break;
+			
+		case Part::SHAPE_CYLINDER:
+			cylinder.setDensity(density);
+			cylinder.initialize(part);
+			iterator = &cylinder;
+			break;
+	}
+}
+
+bool MeshBuilder::PartApices::tryGetNext(Pt3D &point)
+{
+	if (iterator == NULL)
+	{
+		return false;
+	}
+
+	return iterator->tryGetNext(point);
+}
+
+
+
+MeshBuilder::ModelApices::ModelApices(const double _density):
+	Iterator(_density), surface(_density), model(NULL), index(0)
+{}
+
+void MeshBuilder::ModelApices::initialize(const Model *model)
+{
+	this->model = model;
+	surface.setDensity(density);
+	
+	index = 0;
+	if (model->getPartCount() > index)
+	{
+		surface.initialize(model->getPart(index));
+	}
+}
+
+bool MeshBuilder::ModelApices::tryGetNext(Pt3D &point)
+{
+	if (model == NULL)
+	{
+		return false;
+	}
+	
+	do
+	{
+		while (!surface.tryGetNext(point))
+		{
+			index += 1;
+			if (model->getPartCount() > index)
+			{
+				surface.initialize(model->getPart(index));
+			}
+			else
+			{
+				return false;
+			}
+		}
+	} while (GeometryUtils::isPointInsideModelExcludingPart(point, model, index));
+	
+	return true;
+}
+
+
+
+MeshBuilder::BoundingBoxVolume::BoundingBoxVolume(const double _density):
+	Iterator(_density), edge(_density), area(_density), volume(_density), length(0.0), width(0.0), height(0.0)
+{}
+
+bool MeshBuilder::BoundingBoxVolume::initialize(const Model &model)
+{
+	Pt3D lowerBoundary, upperBoundary;
+	ModelGeometryInfo::boundingBox(model, lowerBoundary, upperBoundary);
+	initialize(lowerBoundary, upperBoundary);
+}
+
+bool MeshBuilder::BoundingBoxVolume::initialize(const Pt3D &lowerBoundary, const Pt3D &upperBoundary)
+{
+	length = Pt3D(upperBoundary.x - lowerBoundary.x, 0.0, 0.0);
+	width = Pt3D(0.0, upperBoundary.y - lowerBoundary.y, 0.0);
+	height = Pt3D(0.0, 0.0, upperBoundary.z - lowerBoundary.z);
+
+	edge.setDensity(density);
+	edge.initialize(lowerBoundary, lowerBoundary + length);
+	
+	area.setDensity(density);
+	area.initialize(lowerBoundary, lowerBoundary + width);
+	
+	volume = Segment(density);
+}
+
+bool MeshBuilder::BoundingBoxVolume::tryGetNext(Pt3D &point)
+{
+	while (!volume.tryGetNext(point))
+	{
+		while (!area.tryGetNext(point))
+		{
+			if (edge.tryGetNext(point))
+			{
+				area.initialize(point, point + width);
+			}
+			else
+			{
+				return false;
+			}
+		}
+		
+		volume.initialize(point, point + height);
+	}
+	
+	return true;
+}
Index: cpp/frams/model/geometry/meshbuilder.h
===================================================================
--- cpp/frams/model/geometry/meshbuilder.h	(revision 191)
+++ cpp/frams/model/geometry/meshbuilder.h	(revision 191)
@@ -0,0 +1,248 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#ifndef _MESHBUILDER_H_
+#define _MESHBUILDER_H_
+
+#include <frams/model/model.h>
+#include <frams/model/modelparts.h>
+#include <frams/util/3d.h>
+#include <frams/util/list.h>
+#include "geometryutils.h"
+
+namespace MeshBuilder
+{
+	class Callback
+	{
+		public:
+			virtual void operator()(const Pt3D &point) = 0;
+	};
+	
+	class Iterator
+	{
+		protected:
+			double density;
+		public:
+			Iterator(const double _density);
+			double getDensity() const;
+			double setDensity(const double _density);
+			virtual bool tryGetNext(Pt3D &point) = 0;
+			void forEach(Callback &callback);
+			void addAllPointsToList(SListTempl<Pt3D> &list);
+	};
+
+	class Segment: public Iterator
+	{
+		private:
+			Pt3D point1, point2;
+			int i, I, numberOfPoints;
+		public:
+			Segment(const double _density);
+			void initialize(const Pt3D &point1, const Pt3D &point2, const bool skipFirst = false, const bool skipLast = false, const bool forceOddNumberOfPoints = false);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class RectangleSurface: public Iterator
+	{
+		private:
+			Segment edge1, edge2, area;
+			bool skipVerticalEdges, forceOddNumberOfPoints;
+		public:
+			RectangleSurface(const double _density);
+			void initialize(const Part *part, const CuboidFaces::Face face, const bool skipVerticalEdges = false, const bool skipHorizontalEdges = false, const bool forceOddNumberOfPoints = false);
+			void initialize(const double width, const double height, const Pt3D &position, const Orient &orient, const bool skipVerticalEdges = false, const bool skipHorizontalEdges = false, const bool forceOddNumberOfPoints = false);
+			void initialize(const Pt3D &apex1, const Pt3D &apex2, const Pt3D &apex3, const Pt3D &apex4, const bool skipVerticalEdges = false, const bool skipHorizontalEdges = false, const bool forceOddNumberOfPoints = false);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class CuboidApices: public Iterator
+	{
+		private:
+			Pt3D scale;
+			Pt3D position;
+			Orient orient;
+			Octants::Octant octant;
+		public:
+			CuboidApices();
+			void initialize(const Part *part);
+			void initialize(const Pt3D &scale, const Pt3D &position, const Orient &orient);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class CuboidSurface: public Iterator
+	{
+		private:
+			CuboidApices apices;
+			RectangleSurface rectangle;
+			Iterator *iterator;
+			const Part *part;
+			CuboidFaces::Face face;
+		public:
+			CuboidSurface(const double _density);
+			void initialize(const Part *part);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class EllipseSurface: public Iterator
+	{
+		private:
+			RectangleSurface rectangle;
+			Pt3D position;
+			Orient orient;
+			double inversedSquaredWidth, inversedSquaredHeight;
+		public:
+			EllipseSurface(const double _density);
+			void initialize(const Part *part, const CylinderBases::Base base);
+			void initialize(const double width, const double height, const Pt3D &position, const Orient &orient);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class Ellipse: public Iterator
+	{
+		private:
+			enum Phase { yToZ, zToY, Done };
+			Phase phase;
+			Pt3D position;
+			Orient orient;
+			double width, height, d, p, q, P, Q, D, a, b;
+			QuadrantsYZ::QuadrantYZ quadrant;
+		public:
+			Ellipse(const double _density);
+			void initialize(const Part *part, const CylinderBases::Base base);
+			void initialize(const double width, const double height, const Pt3D &position, const Orient &orient);
+			bool tryGetNext(Pt3D &point);
+		private:
+			void calculateAndSetSegmentationDistance();
+			void initializePhase(Phase ph);
+			void setPoint(Pt3D &point);
+			void findNextPQAndPhase();
+			int findLastPQOfPhase(Phase ph, double &lp, double &lq);
+	};
+
+	class CylinderEdges: public Iterator
+	{
+		private:
+			Ellipse ellipse;
+			Pt3D edge, length;
+			CylinderBases::Base base;
+		public:
+			CylinderEdges(const double _density);
+			void initialize(const Part *part);
+			void initialize(const Pt3D &scale, const Pt3D &position, const Orient &orient);
+			void initialize(const double length, const double width, const double height, const Pt3D &position, const Orient &orient);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class CylinderWallSurface: public Iterator
+	{
+		private:
+			Ellipse edge;
+			Pt3D length;
+			Segment area;
+		public:
+			CylinderWallSurface(const double _density);
+			void initialize(const Part *part);
+			void initialize(const Pt3D &scale, const Pt3D &position, const Orient &orient);
+			void initialize(const double length, const double width, const double height, const Pt3D &position, const Orient &orient);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class CylinderSurface: public Iterator
+	{
+		private:
+			CylinderWallSurface wall;
+			EllipseSurface ellipse;
+			Iterator *iterator;
+			const Part *part;
+			CylinderBases::Base base;
+		public:
+			CylinderSurface(const double _density);
+			void initialize(const Part *part);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class EllipsoidSurface: public Iterator
+	{
+		private:
+			enum Phase { X, Y, Z, Done };
+			Phase phase;
+			Pt3D edge, area, scale, limit;
+			double d;
+			const Part *part;
+			Octants::Octant octant;
+		public:
+			EllipsoidSurface(const double _density);
+			void initialize(const Part *part);
+			bool tryGetNext(Pt3D &point);
+		private:
+			void initializePhase(Phase ph);
+			void setPoint(Pt3D &point);
+			void proceedToNextOctant();
+			void findNextAreaEdgeAndPhase();
+	};
+
+	class PartSurface: public Iterator
+	{
+		private:
+			CuboidSurface cuboid;
+			CylinderSurface cylinder;
+			EllipsoidSurface ellipsoid;
+			Iterator *iterator;
+		public:
+			PartSurface(const double _density);
+			void initialize(const Part *part);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class ModelSurface: public Iterator
+	{
+		private:
+			PartSurface surface;
+			const Model *model;
+			int index;
+		public:
+			ModelSurface(const double _density);
+			void initialize(const Model *model);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class PartApices: public Iterator
+	{
+		private:
+			CuboidApices cuboid;
+			CylinderEdges cylinder;
+			EllipsoidSurface ellipsoid;
+			Iterator *iterator;
+		public:
+			PartApices(const double _density);
+			void initialize(const Part *part);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class ModelApices: public Iterator
+	{
+		private:
+			PartApices surface;
+			const Model *model;
+			int index;
+		public:
+			ModelApices(const double _density);
+			void initialize(const Model *model);
+			bool tryGetNext(Pt3D &point);
+	};
+
+	class BoundingBoxVolume: Iterator
+	{
+		private:
+			Segment edge, area, volume;
+			Pt3D length, width, height;
+		public:
+			BoundingBoxVolume(const double _density);
+			bool initialize(const Model &model);
+			bool initialize(const Pt3D &lowerBoundary, const Pt3D &upperBoundary);
+			bool tryGetNext(Pt3D &point);
+	};
+};
+
+#endif
Index: cpp/frams/model/geometry/modelgeometryinfo.cpp
===================================================================
--- cpp/frams/model/geometry/modelgeometryinfo.cpp	(revision 191)
+++ cpp/frams/model/geometry/modelgeometryinfo.cpp	(revision 191)
@@ -0,0 +1,186 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#include "modelgeometryinfo.h"
+#include <frams/model/geometry/geometryutils.h>
+#include <frams/model/geometry/meshbuilder.h>
+
+void ModelGeometryInfo::findSizesAndAxesOfModel(const Model &model, const double density,
+	Pt3D &sizes, Orient &axes)
+{
+	SListTempl<Pt3D> points;
+	MeshBuilder::ModelApices apices(density);
+	apices.initialize(&model);
+	apices.addAllPointsToList(points);
+	GeometryUtils::findSizesAndAxesOfPointsGroup(points, sizes, axes);
+}
+
+bool ModelGeometryInfo::boundingBox(const Model &model, Pt3D &lowerBoundary, Pt3D &upperBoundary)
+{
+	if (model.getPartCount() == 0)
+	{
+		return false;
+	}
+	
+	boundingBox(model.getPart(0), lowerBoundary, upperBoundary);
+	
+	for (int i = 1; i < model.getPartCount(); i++)
+	{
+		Pt3D partLowerBoundary, partUpperBoundary;
+		boundingBox(model.getPart(i), partLowerBoundary, partUpperBoundary);
+		
+		lowerBoundary.getMin(partLowerBoundary);
+		upperBoundary.getMax(partUpperBoundary);
+	}
+	
+	return true;
+}
+
+void ModelGeometryInfo::boundingBox(const Part *part, Pt3D &lowerBoundary, Pt3D &upperBoundary)
+{
+	lowerBoundary.x = upperBoundary.x = part->p.x;
+	lowerBoundary.y = upperBoundary.y = part->p.y;
+	lowerBoundary.z = upperBoundary.z = part->p.z;
+	
+	for (Octants::Octant o = Octants::FIRST; o < Octants::NUMBER; o = Octants::Octant(o+1))
+	{
+		Pt3D vertex = part->scale;
+		vertex.x *= Octants::isPositiveX(o) ? +1 : -1;
+		vertex.y *= Octants::isPositiveY(o) ? +1 : -1;
+		vertex.z *= Octants::isPositiveZ(o) ? +1 : -1;
+		
+		vertex = part->o.transform(vertex) + part->p;
+		
+		lowerBoundary.getMin(vertex);
+		upperBoundary.getMax(vertex);
+	}
+}
+
+double ModelGeometryInfo::volume(const Model &model, const double density)
+{
+	Pt3D lowerBoundary, upperBoundary;
+	boundingBox(model, lowerBoundary, upperBoundary);
+	
+	MeshBuilder::BoundingBoxVolume iterator(density);
+	iterator.initialize(lowerBoundary, upperBoundary);
+	
+	Pt3D point;
+	int allPoints = 0, pointsInsideModel = 0;
+	
+	while (iterator.tryGetNext(point))
+	{
+		allPoints += 1;
+		pointsInsideModel += GeometryUtils::isPointInsideModel(point, model) ? 1 : 0;
+	}
+	
+	double boundingBoxVolume
+		= (upperBoundary.x - lowerBoundary.x)
+		* (upperBoundary.y - lowerBoundary.y)
+		* (upperBoundary.z - lowerBoundary.z);
+	
+	return boundingBoxVolume * (double)pointsInsideModel / (double)allPoints;
+}
+
+double ModelGeometryInfo::area(const Model &model, const double density)
+{
+	double area = 0.0;
+	
+	for (int partIndex = 0; partIndex < model.getPartCount(); partIndex+=1)
+	{
+		area += externalAreaOfPart(model, partIndex, density);
+	}
+	
+	return area;
+}
+
+double ModelGeometryInfo::externalAreaOfPart(const Model &model, const int partIndex, const double density)
+{
+	switch (model.getPart(partIndex)->shape)
+	{
+		case Part::SHAPE_ELLIPSOID:
+			return externalAreaOfEllipsoid(model, partIndex, density);
+			
+		case Part::SHAPE_CUBOID:
+			return externalAreaOfCuboid(model, partIndex, density);
+			
+		case Part::SHAPE_CYLINDER:
+			return externalAreaOfCylinder(model, partIndex, density);
+	}
+}
+
+double externalAreaOfSurface(const Model &model, const int partIndex, MeshBuilder::Iterator &surface, const double area)
+{
+	Pt3D point;
+	int all = 0, sur = 0;
+	
+	while (surface.tryGetNext(point))
+	{
+		all += 1;
+		sur += GeometryUtils::isPointInsideModelExcludingPart(point, &model, partIndex) ? 0 : 1;
+	}
+	
+	return sur * area / all;
+}
+
+double ModelGeometryInfo::externalAreaOfEllipsoid(const Model &model, const int partIndex, const double density)
+{
+	Part *part = model.getPart(partIndex);
+	
+	MeshBuilder::EllipsoidSurface ellipsoid(density);
+	ellipsoid.initialize(part);
+	
+	double area = GeometryUtils::ellipsoidArea(part->scale);
+	
+	return externalAreaOfSurface(model, partIndex, ellipsoid, area);
+}
+
+double ModelGeometryInfo::externalAreaOfCuboid(const Model &model, const int partIndex, const double density)
+{
+	Part *part = model.getPart(partIndex);
+	
+	double externalArea = 0.0;
+	
+	MeshBuilder::RectangleSurface rectangle(density);
+	
+	for (CuboidFaces::Face f = CuboidFaces::FIRST; f < CuboidFaces::NUMBER; f = CuboidFaces::Face(f+1))
+	{
+		rectangle.initialize(part, f);
+		
+		double area = 4.0;
+		area *= !CuboidFaces::isX(f) ? part->scale.x : 1.0;
+		area *= !CuboidFaces::isY(f) ? part->scale.y : 1.0;
+		area *= !CuboidFaces::isZ(f) ? part->scale.z : 1.0;
+		
+		externalArea += externalAreaOfSurface(model, partIndex, rectangle, area);
+	}
+	
+	return externalArea;
+}
+
+double ModelGeometryInfo::externalAreaOfCylinder(const Model &model, const int partIndex, const double density)
+{
+	Part *part = model.getPart(partIndex);
+	
+	double externalArea = 0.0;
+	
+	MeshBuilder::EllipseSurface ellipse(density);
+	
+	double area = M_PI * part->scale.y * part->scale.z;
+	
+	for (CylinderBases::Base b = CylinderBases::FIRST; b < CylinderBases::NUMBER; b = CylinderBases::Base(b+1))
+	{
+		ellipse.initialize(part, b);
+		
+		externalArea += externalAreaOfSurface(model, partIndex, ellipse, area);
+	}
+	
+	MeshBuilder::CylinderWallSurface cylinderWall(density);
+	cylinderWall.initialize(part);
+	
+	area = 2.0*part->scale.x * GeometryUtils::ellipsePerimeter(part->scale.y, part->scale.z);
+	
+	externalArea += externalAreaOfSurface(model, partIndex, cylinderWall, area);
+	
+	return externalArea;
+}
Index: cpp/frams/model/geometry/modelgeometryinfo.h
===================================================================
--- cpp/frams/model/geometry/modelgeometryinfo.h	(revision 191)
+++ cpp/frams/model/geometry/modelgeometryinfo.h	(revision 191)
@@ -0,0 +1,26 @@
+// This file is a part of the Framsticks GDK.
+// Copyright (C) 2002-2014  Maciej Komosinski and Szymon Ulatowski.  See LICENSE.txt for details.
+// Refer to http://www.framsticks.com/ for further information.
+
+#ifndef _MODELGEOMETRYINFO_H_
+#define _MODELGEOMETRYINFO_H_
+
+#include <frams/model/model.h>
+#include <frams/util/3d.h>
+#include <frams/util/list.h>
+
+namespace ModelGeometryInfo
+{
+	void findSizesAndAxesOfModel(const Model &model, const double density, Pt3D &sizes,
+		Orient &axes);
+	bool boundingBox(const Model &model, Pt3D &lowerBoundary, Pt3D &upperBoundary);
+	void boundingBox(const Part *part, Pt3D &lowerBoundary, Pt3D &upperBoundary);
+	double volume(const Model &model, const double density);
+	double area(const Model &model, const double density);
+	double externalAreaOfPart(const Model &model, const int partIndex, const double density);
+	double externalAreaOfEllipsoid(const Model &model, const int partIndex, const double density);
+	double externalAreaOfCuboid(const Model &model, const int partIndex, const double density);
+	double externalAreaOfCylinder(const Model &model, const int partIndex, const double density);
+}
+
+#endif
