Index: cpp/frams/Makefile-SDK-files
===================================================================
--- cpp/frams/Makefile-SDK-files	(revision 1012)
+++ cpp/frams/Makefile-SDK-files	(revision 1017)
@@ -3,4 +3,6 @@
 # ALL_DIRS is later expanded by the shell, no spaces/newlines allowed, or it breaks
 ALL_DIRS={common,PrintFloat,frams,frams/canvas,frams/config,common/loggers,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/fn,frams/genetics/fF,frams/genetics/fT,frams/genetics/fB,frams/genetics/fH,frams/genetics/fL,frams/genetics/fS,frams/model,frams/neuro,frams/neuro/impl,frams/param,frams/test,frams/util,frams/vm/classes,common/virtfile,frams/_demos,frams/model/geometry,frams/_demos/geometry,frams/model/similarity,frams/model/similarity/hungarian,frams/model/similarity/SVD}
+
+GEOMETRY_OBJS=frams/model/geometry/meshbuilder.o frams/model/geometry/modelgeometryinfo.o frams/model/geometry/geometryutils.o
 
 GENMANF4=frams/genetics/f4/f4_oper.o
@@ -22,5 +24,5 @@
 CONVFH=frams/genetics/fH/fH_conv.o frams/genetics/fH/fH_general.o frams/param/mutableparam.o
 CONVFL=frams/genetics/fL/fL_conv.o frams/genetics/fL/fL_general.o frams/genetics/fL/fL_matheval.o
-CONVFS=frams/genetics/fS/fS_conv.o frams/genetics/fS/fS_general.o 
+CONVFS=frams/genetics/fS/fS_conv.o frams/genetics/fS/fS_general.o $(GEOMETRY_OBJS)
 
 # $(sort - remove duplicates
@@ -40,6 +42,4 @@
 
 SDK_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/log.o common/util-string.o common/util-file.o common/nonstd_stdio.o frams/util/callbacks.o frams/param/syntparam.o frams/util/multirange.o frams/util/multimap.o frams/param/paramtabobj.o common/loggers/loggers.o frams/param/paramobj.o frams/genetics/genooperators.o common/nonstd_math.o frams/util/validitychecks.o common/Convert.o frams/util/rndutil.o common/virtfile/stringfile.o common/virtfile/stdinoutfilesystem.o $(PRINTFLOAT_OBJS)
-
-GEOMETRY_OBJS=frams/model/geometry/meshbuilder.o frams/model/geometry/modelgeometryinfo.o frams/model/geometry/geometryutils.o
 
 STDOUT_LOGGER_OBJS=common/virtfile/virtfile.o common/loggers/loggertostdout.o common/console.o 
@@ -84,11 +84,11 @@
 NEURO_LAYOUT_TEST_OBJS= frams/_demos/neuro_layout_test.o $(STDOUT_LOGGER_OBJS) $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(NN_LAYOUT_OBJS)
 
-GEOMETRY_INFO_TEST_OBJS=frams/_demos/geometry/info_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS)
+GEOMETRY_INFO_TEST_OBJS=$(sort frams/_demos/geometry/info_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS))
 
-GEOMETRY_SURFACE_TEST_OBJS=frams/_demos/geometry/surface_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS)
+GEOMETRY_SURFACE_TEST_OBJS=$(sort frams/_demos/geometry/surface_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS))
 
-GEOMETRY_VOLUME_TEST_OBJS=frams/_demos/geometry/volume_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS)
+GEOMETRY_VOLUME_TEST_OBJS=$(sort frams/_demos/geometry/volume_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS))
 
-GEOMETRY_APICES_TEST_OBJS=frams/_demos/geometry/apices_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS)
+GEOMETRY_APICES_TEST_OBJS=$(sort frams/_demos/geometry/apices_test.o frams/_demos/geometry/geometrytestutils.o frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) $(GEOMETRY_OBJS) $(STDOUT_LOGGER_OBJS))
 
 SIMIL_TEST_OBJS=frams/_demos/simil_test.o $(SIMILARITY_OBJS) frams/_demos/genotypeloader.o frams/_demos/genotypemini.o frams/param/multiparamload.o common/virtfile/stdiofile.o $(STDOUT_LOGGER_OBJS) $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS)
@@ -98,5 +98,5 @@
 PARAMTREE_STDIN_TEST_OBJS=frams/_demos/paramtree_stdin_test.o frams/_demos/paramtree_print.o $(STDOUT_LOGGER_OBJS) common/virtfile/stdiofile.o $(SDK_OBJS) frams/param/paramtree.o
 
-PARAMTREE_PARAMLIST_TEST_OBJS=frams/_demos/paramtree_paramlist_test.o frams/_demos/paramtree_print.o $(STDOUT_LOGGER_OBJS) common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) frams/neuro/neuroimpl.o frams/neuro/neurofactory.o frams/neuro/impl/neuroimpl-simple.o frams/neuro/impl/neuroimpl-channels.o frams/neuro/impl/neuroimpl-fuzzy.o frams/neuro/impl/neuroimpl-fuzzy-f0.o $(GENOTYPE_LOADER_OBJS) $(GEOMETRY_OBJS) frams/model/geometry/modelgeoclass.o frams/model/modelobj.o frams/param/paramtree.o
+PARAMTREE_PARAMLIST_TEST_OBJS=$(sort frams/_demos/paramtree_paramlist_test.o frams/_demos/paramtree_print.o $(STDOUT_LOGGER_OBJS) common/virtfile/stdiofile.o $(SDK_OBJS) $(GENOCONV_AND_GENMAN_SDK_OBJS) frams/neuro/neuroimpl.o frams/neuro/neurofactory.o frams/neuro/impl/neuroimpl-simple.o frams/neuro/impl/neuroimpl-channels.o frams/neuro/impl/neuroimpl-fuzzy.o frams/neuro/impl/neuroimpl-fuzzy-f0.o $(GENOTYPE_LOADER_OBJS) $(GEOMETRY_OBJS) frams/model/geometry/modelgeoclass.o frams/model/modelobj.o frams/param/paramtree.o)
 
 MUTABLEPARAM_TEST_OBJS=frams/_demos/mutableparam_test.o $(STDOUT_LOGGER_OBJS) common/virtfile/stdiofile.o $(SDK_OBJS) frams/param/mutableparam.o frams/param/mutparamlist.o frams/param/paramtrans.o
Index: cpp/frams/genetics/defgenoconv.cpp
===================================================================
--- cpp/frams/genetics/defgenoconv.cpp	(revision 1012)
+++ cpp/frams/genetics/defgenoconv.cpp	(revision 1017)
@@ -104,5 +104,5 @@
 #endif
 #ifdef USE_GENCONV_fS0
-	addConverter(new GenoConv_fS0); //solids
+	addConverter(new GenoConv_fS0s); //solids
 #endif
 
Index: cpp/frams/genetics/fS/fS_conv.cpp
===================================================================
--- cpp/frams/genetics/fS/fS_conv.cpp	(revision 1012)
+++ cpp/frams/genetics/fS/fS_conv.cpp	(revision 1017)
@@ -5,5 +5,5 @@
 #include "fS_conv.h"
 
-SString GenoConv_fS0::convert(SString &i, MultiMap *map, bool using_checkpoints)
+SString GenoConv_fS0s::convert(SString &i, MultiMap *map, bool using_checkpoints)
 {
 	fS_Genotype *genotype;
@@ -14,5 +14,5 @@
 	catch (fS_Exception &e)
 	{
-		logPrintf("GenoConv_fS0", "convert", LOG_ERROR, e.what());
+		logPrintf("GenoConv_fS0s", "convert", LOG_ERROR, e.what());
 		return SString();
 	}
Index: cpp/frams/genetics/fS/fS_conv.h
===================================================================
--- cpp/frams/genetics/fS/fS_conv.h	(revision 1012)
+++ cpp/frams/genetics/fS/fS_conv.h	(revision 1017)
@@ -12,8 +12,8 @@
  * Genotype converter from fS to f0s.
  */
-class GenoConv_fS0 : public GenoConverter
+class GenoConv_fS0s : public GenoConverter
 {
 public:
-	GenoConv_fS0() : GenoConverter()
+	GenoConv_fS0s() : GenoConverter()
 	{
 		name = "Solid encoding";
@@ -27,5 +27,5 @@
 	SString convert(SString &i, MultiMap *map, bool using_checkpoints);
 
-	~GenoConv_fS0()
+	~GenoConv_fS0s()
 	{};
 };
Index: cpp/frams/genetics/fS/fS_general.cpp
===================================================================
--- cpp/frams/genetics/fS/fS_general.cpp	(revision 1012)
+++ cpp/frams/genetics/fS/fS_general.cpp	(revision 1017)
@@ -16,22 +16,62 @@
 int fS_Genotype::precision = 4;
 bool fS_Genotype::TURN_WITH_ROTATION = false;
+std::map<string, double> Node::minValues;
+std::map<string, double> Node::defaultValues;
+std::map<string, double> Node::maxValues;
 
 void Node::prepareParams()
 {
-	defaultValues = {
-			{INGESTION, Model::getDefPart().ingest},
-			{FRICTION,  Model::getDefPart().friction},
-			{STIFFNESS, Model::getDefJoint().stif},
-			{ROT_X,     0.0},
-			{ROT_Y,     0.0},
-			{ROT_Z,     0.0},
-			{RX,        0.0},
-			{RY,        0.0},
-			{RZ,        0.0},
-			{SIZE,      1.0},
-			{SIZE_X,    Model::getDefPart().scale.x},
-			{SIZE_Y,    Model::getDefPart().scale.y},
-			{SIZE_Z,    Model::getDefPart().scale.z}
-	};
+	if(minValues.empty())
+	{
+		minValues = {
+				{INGESTION, Model::getMinPart().ingest},
+				{FRICTION,  Model::getMinPart().friction},
+				{ROT_X,     -M_PI},
+				{ROT_Y,     -M_PI},
+				{ROT_Z,     -M_PI},
+				{RX,        -M_PI},
+				{RY,        -M_PI},
+				{RZ,        -M_PI},
+				{SIZE,      0.01},
+				{SIZE_X,    Model::getMinPart().scale.x},
+				{SIZE_Y,    Model::getMinPart().scale.y},
+				{SIZE_Z,    Model::getMinPart().scale.z}
+		};
+	}
+
+	if(maxValues.empty())
+	{
+		maxValues = {
+				{INGESTION, Model::getMaxPart().ingest},
+				{FRICTION,  Model::getMaxPart().friction},
+				{ROT_X,     M_PI},
+				{ROT_Y,     M_PI},
+				{ROT_Z,     M_PI},
+				{RX,        M_PI},
+				{RY,        M_PI},
+				{RZ,        M_PI},
+				{SIZE,      100.0},
+				{SIZE_X,    Model::getMaxPart().scale.x},
+				{SIZE_Y,    Model::getMaxPart().scale.y},
+				{SIZE_Z,    Model::getMaxPart().scale.z}
+		};
+	}
+	if(defaultValues.empty())
+	{
+		defaultValues = {
+				{INGESTION, Model::getDefPart().ingest},
+				{FRICTION,  Model::getDefPart().friction},
+				{ROT_X,     0.0},
+				{ROT_Y,     0.0},
+				{ROT_Z,     0.0},
+				{RX,        0.0},
+				{RY,        0.0},
+				{RZ,        0.0},
+				{SIZE,      1.0},
+				{SIZE_X,    Model::getDefPart().scale.x},
+				{SIZE_Y,    Model::getDefPart().scale.y},
+				{SIZE_Z,    Model::getDefPart().scale.z}
+		};
+	}
 }
 
@@ -58,5 +98,4 @@
 	fr = _state->fr;
 	s = _state->s;
-	stif = _state->stif;
 }
 
@@ -259,4 +298,5 @@
 	if(paramsEndIndex == -1)
 		throw fS_Exception("Lacking param end sign", restOfGenotype.start);
+
 	for (int i = 0; i < int(separators.size()) - 1; i++)
 	{
@@ -297,15 +337,22 @@
 }
 
-double Node::getParam(string key)
+double Node::getParam(const string &key)
 {
 	auto item = params.find(key);
 	if (item != params.end())
 		return item->second;
-	else
-		return defaultValues.at(key);
-}
-
-
-void Node::getState(State *_state)
+	return defaultValues.at(key);
+}
+
+double Node::getParam(const string &key, double defaultValue)
+{
+	auto item = params.find(key);
+	if (item != params.end())
+		return item->second;
+	return defaultValue;
+}
+
+
+void Node::getState(State *_state, bool calculateLocation)
 {
 	if (state != nullptr)
@@ -328,9 +375,7 @@
 		else if (mod == MODIFIERS[2])
 			state->s *= multiplier;
-		else if (mod == MODIFIERS[3])
-			state->stif *= multiplier;
-	}
-
-	if (parent != nullptr)
+	}
+
+	if (parent != nullptr && calculateLocation)
 	{
 		// Rotate
@@ -341,5 +386,5 @@
 	}
 	for (int i = 0; i < int(children.size()); i++)
-		children[i]->getState(state);
+		children[i]->getState(state, calculateLocation);
 }
 
@@ -388,11 +433,10 @@
 }
 
-Pt3D Node::calculateSize()
+void Node::calculateSize(Pt3D &scale)
 {
 	double sizeMultiplier = getParam(SIZE) * state->s;
-	double sx = getParam(SIZE_X) * sizeMultiplier;
-	double sy = getParam(SIZE_Y) * sizeMultiplier;
-	double sz = getParam(SIZE_Z) * sizeMultiplier;
-	return Pt3D(sx, sy, sz);
+	scale.x = getParam(SIZE_X) * sizeMultiplier;
+	scale.y = getParam(SIZE_Y) * sizeMultiplier;
+	scale.z = getParam(SIZE_Z) * sizeMultiplier;
 }
 
@@ -400,5 +444,6 @@
 {
 	double result;
-	Pt3D size = calculateSize();
+	Pt3D size;
+	calculateSize(size);
 	double radiiProduct = size.x * size.y * size.z;
 	switch (partType)
@@ -421,5 +466,6 @@
 bool Node::isPartSizeValid()
 {
-	Pt3D size = calculateSize();
+	Pt3D size;
+	calculateSize(size);
 	double volume = calculateVolume();
 	Part_MinMaxDef minP = Model::getMinPart();
@@ -449,10 +495,10 @@
 Pt3D Node::getVectorRotation()
 {
-	return Pt3D(getParam(ROT_X), getParam(ROT_Y), getParam(ROT_Z));
+	return Pt3D(getParam(ROT_X, 0.0), getParam(ROT_Y, 0.0), getParam(ROT_Z, 0.0));
 }
 
 Pt3D Node::getRotation()
 {
-	Pt3D rotation = Pt3D(getParam(RX), getParam(RY), getParam(RZ));
+	Pt3D rotation = Pt3D(getParam(RX, 0.0), getParam(RY, 0.0), getParam(RZ, 0.0));
 	if(fS_Genotype::TURN_WITH_ROTATION)
 		rotation += getVectorRotation();
@@ -492,14 +538,9 @@
 {
 	part = new Part(partType);
-	part->p = Pt3D(state->location.x,
-				   state->location.y,
-				   state->location.z);
+	part->p = Pt3D(state->location);
 
 	part->friction = getParam(FRICTION) * state->fr;
 	part->ingest = getParam(INGESTION) * state->ing;
-	Pt3D size = calculateSize();
-	part->scale.x = size.x;
-	part->scale.y = size.y;
-	part->scale.z = size.z;
+	calculateSize(part->scale);
 	part->setRot(getRotation());
 }
@@ -508,7 +549,4 @@
 {
 	Joint *j = new Joint();
-	j->stif = getParam(STIFFNESS) * state->stif;
-	j->rotstif = j->stif;
-
 	j->attachToParts(parent->part, part);
 	switch (joint)
@@ -650,8 +688,8 @@
 }
 
-void fS_Genotype::getState()
+void fS_Genotype::getState(bool calculateLocation)
 {
 	State *initialState = new State(Pt3D(0), Pt3D(1, 0, 0));
-	startNode->getState(initialState);
+	startNode->getState(initialState, calculateLocation);
 }
 
@@ -662,5 +700,5 @@
 	model.open(using_checkpoints);
 
-	getState();
+	getState(true);
 	startNode->buildModel(model, nullptr);
 	buildNeuroConnections(model);
@@ -808,5 +846,5 @@
 int fS_Genotype::checkValidityOfPartSizes()
 {
-	getState();
+	getState(false);
 	vector<Node*> nodes = getAllNodes();
 	for (int i = 0; i < int(nodes.size()); i++)
Index: cpp/frams/genetics/fS/fS_general.h
===================================================================
--- cpp/frams/genetics/fS/fS_general.h	(revision 1012)
+++ cpp/frams/genetics/fS/fS_general.h	(revision 1017)
@@ -54,10 +54,5 @@
 #define RZ "rz"
 //@}
-/** @name Macros and values used in collision detection */
-//@{
-#define DISJOINT 0
-#define COLLISION 1
-#define ADJACENT 2
-//@}
+
 
 #define HINGE_X 'b'
@@ -87,8 +82,7 @@
 const string ALL_JOINTS = "abc";
 const int JOINT_COUNT = JOINTS.length();
-const string MODIFIERS = "IFST";
+const string MODIFIERS = "IFS";
 const char SIZE_MODIFIER = 's';
-const vector<string> PARAMS {INGESTION, FRICTION, ROT_X, ROT_Y, ROT_Z, RX, RY, RZ, SIZE, SIZE_X, SIZE_Y, SIZE_Z,
-							 STIFFNESS};
+const vector<string> PARAMS {INGESTION, FRICTION, ROT_X, ROT_Y, ROT_Z, RX, RY, RZ, SIZE, SIZE_X, SIZE_Y, SIZE_Z};
 const vector<string> SIZE_PARAMS {SIZE, SIZE_X, SIZE_Y, SIZE_Z};
 
@@ -214,5 +208,4 @@
 	double ing = 1.0;      /// Ingestion multiplier
 	double s = 1.0;      /// Size multipliers
-	double stif = 1.0;	/// Stiffness multipliers
 
 	State(State *_state); /// Derive the state from parent
@@ -272,5 +265,4 @@
 	Part *part;     /// A part object built from node. Used in building the Model
 	int partCodeLen; /// The length of substring that directly describes the corresponding part
-	std::map<string, double> defaultValues;
 	GenotypeParams genotypeParams;
 
@@ -335,5 +327,5 @@
 	 * @param _state state of the parent
 	 */
-	void getState(State *_state);
+	void getState(State *_state, bool calculateLocation);
 
 	/**
@@ -370,4 +362,7 @@
 
 public:
+	static std::map<string, double> minValues;
+	static std::map<string, double> defaultValues;
+	static std::map<string, double> maxValues;
 	char joint = DEFAULT_JOINT;           /// Set of all joints
 	Part::Shape partType;  /// The type of the part
@@ -389,5 +384,5 @@
 	 * @return The effective size
 	 */
-	Pt3D calculateSize();
+	void calculateSize(Pt3D &scale);
 
 	/**
@@ -407,5 +402,6 @@
 	 * @return the param value
 	 */
-	double getParam(string key);
+	double getParam(const string &key);
+	double getParam(const string &key, double defaultValue);
 };
 
@@ -454,5 +450,5 @@
 	~fS_Genotype();
 
-	void getState();
+	void getState(bool calculateLocation);
 
 	/**
Index: cpp/frams/genetics/fS/fS_oper.cpp
===================================================================
--- cpp/frams/genetics/fS/fS_oper.cpp	(revision 1012)
+++ cpp/frams/genetics/fS/fS_oper.cpp	(revision 1017)
@@ -24,5 +24,5 @@
 				{"fS_mut_mod_neuro_conn",   0, 0, "Modify neuron connection",    "f 0 100 10", FIELD(prob[FS_MOD_NEURO_CONNECTION]), "mutation: probability of changing a neuron connection",},
 				{"fS_mut_add_neuro_conn",   0, 0, "Add neuron connection",       "f 0 100 10", FIELD(prob[FS_ADD_NEURO_CONNECTION]), "mutation: probability of adding a neuron connection",},
-				{"fS_mut_rem neuro_conn",   0, 0, "Remove neuron connection",    "f 0 100 10", FIELD(prob[FS_REM_NEURO_CONNECTION]), "mutation: probability of removing a neuron connection",},
+				{"fS_mut_rem_neuro_conn",   0, 0, "Remove neuron connection",    "f 0 100 10", FIELD(prob[FS_REM_NEURO_CONNECTION]), "mutation: probability of removing a neuron connection",},
 				{"fS_mut_mod_neuro_params", 0, 0, "Modify neuron params",        "f 0 100 10", FIELD(prob[FS_MOD_NEURO_PARAMS]),     "mutation: probability of changing a neuron param",},
 				{"fS_circle_section",       0, 0, "Ensure circle section",       "d 0 1 1",    FIELD(ensureCircleSection),           "Ensure that ellipsoids and cylinders have circle cross-section"},
@@ -35,43 +35,6 @@
 #undef FIELDSTRUCT
 
-
-void GenoOper_fS::prepareParams()
-{
-	minValues = {
-			{INGESTION, Model::getMinPart().ingest},
-			{FRICTION,  Model::getMinPart().friction},
-			{STIFFNESS, 0.1},
-			{ROT_X,     -M_PI},
-			{ROT_Y,     -M_PI},
-			{ROT_Z,     -M_PI},
-			{RX,        -M_PI},
-			{RY,        -M_PI},
-			{RZ,        -M_PI},
-			{SIZE,      0.01},
-			{SIZE_X,    Model::getMinPart().scale.x},
-			{SIZE_Y,    Model::getMinPart().scale.y},
-			{SIZE_Z,    Model::getMinPart().scale.z}
-	};
-
-	maxValues = {
-			{INGESTION, Model::getMaxPart().ingest},
-			{FRICTION,  Model::getMaxPart().friction},
-			{STIFFNESS, 0.5},
-			{ROT_X,     M_PI},
-			{ROT_Y,     M_PI},
-			{ROT_Z,     M_PI},
-			{RX,        M_PI},
-			{RY,        M_PI},
-			{RZ,        M_PI},
-			{SIZE,      100.0},
-			{SIZE_X,    Model::getMaxPart().scale.x},
-			{SIZE_Y,    Model::getMaxPart().scale.y},
-			{SIZE_Z,    Model::getMaxPart().scale.z}
-	};
-}
-
 GenoOper_fS::GenoOper_fS()
 {
-	prepareParams();
 	par.setParamTab(genooper_fS_paramtab);
 	par.select(this);
@@ -332,5 +295,5 @@
 bool GenoOper_fS::addPart(fS_Genotype &geno, const vector <Part::Shape> &availablePartShapes, bool mutateSize)
 {
-	geno.getState();
+	geno.getState(false);
 	Node *node = geno.chooseNode();
 	char partType = SHAPE_TO_GENE.at(availablePartShapes[rndUint(availablePartShapes.size())]);
@@ -375,5 +338,5 @@
 	if (mutateSize)
 	{
-		geno.getState();
+		geno.getState(false);
 		mutateSizeParam(newNode, SIZE_X, true);
 		mutateSizeParam(newNode, SIZE_Y, true);
@@ -387,5 +350,6 @@
 	Node *randomNode, *selectedChild;
 	// Choose a parent with children
-	for (int i = 0; i < mutationTries; i++)
+	// It may be difficult to choose a eligible node, so the number of tries should be high
+	for (int i = 0; i < 10 * mutationTries; i++)
 	{
 		randomNode = geno.chooseNode();
@@ -425,5 +389,5 @@
 #endif
 
-		geno.getState();
+		geno.getState(false);
 		double sizeMultiplier = randomNode->getParam(SIZE) * randomNode->state->s;
 		double relativeVolume = randomNode->calculateVolume() / pow(sizeMultiplier, 3.0);
@@ -485,5 +449,4 @@
 	bool isRadius = isRadiusOfBase || key == SIZE_X;
 	if (ensureCircleSection && isRadius)
-	if (ensureCircleSection && isRadius)
 	{
 		if (randomNode->partType == Part::Shape::SHAPE_ELLIPSOID)
@@ -493,6 +456,7 @@
 	}
 	// Add modified default value for param
-	randomNode->params[key] = mutateCreep('f', randomNode->defaultValues.at(key), minValues.at(key), maxValues.at(key), true);
-	return true;
+	randomNode->params[key] = randomNode->defaultValues.at(key);
+	geno.getState(false);
+	return mutateParamValue(randomNode, key);
 }
 
@@ -508,6 +472,14 @@
 			auto it = randomNode->params.begin();
 			advance(it, rndUint(paramCount));
-			randomNode->params.erase(it->first);
-			return true;
+			string key = it->first;
+			double value = it->second;
+
+			randomNode->params.erase(key);
+			if(geno.checkValidityOfPartSizes() == 0)
+				return true;
+			else
+			{
+				randomNode->params[key] = value;
+			}
 		}
 	}
@@ -515,7 +487,19 @@
 }
 
+
+bool GenoOper_fS::mutateParamValue(Node *node, string key)
+{
+	// Do not allow invalid changes in part size
+	if (std::find(SIZE_PARAMS.begin(), SIZE_PARAMS.end(), key) == SIZE_PARAMS.end())
+	{
+		node->params[key] = GenoOperators::mutateCreep('f', node->getParam(key), Node::minValues.at(key), Node::maxValues.at(key), true);
+		return true;
+	} else
+		return mutateSizeParam(node, key, ensureCircleSection);
+}
+
 bool GenoOper_fS::changeParam(fS_Genotype &geno)
 {
-	geno.getState();
+	geno.getState(false);
 	for (int i = 0; i < mutationTries; i++)
 	{
@@ -526,12 +510,5 @@
 			auto it = randomNode->params.begin();
 			advance(it, rndUint(paramCount));
-
-			// Do not allow invalid changes in part size
-			if (std::find(SIZE_PARAMS.begin(), SIZE_PARAMS.end(), it->first) == SIZE_PARAMS.end())
-			{
-				it->second = GenoOperators::mutateCreep('f', it->second, minValues.at(it->first), maxValues.at(it->first), true);
-				return true;
-			} else
-				return mutateSizeParam(randomNode, it->first, ensureCircleSection);
+			return mutateParamValue(randomNode, it->first);
 		}
 	}
@@ -543,4 +520,6 @@
 	Node *randomNode = geno.chooseNode();
 	char randomModifier = MODIFIERS[rndUint(MODIFIERS.length())];
+	int oldValue = randomNode->modifiers[randomModifier];
+
 	randomNode->modifiers[randomModifier] += rndUint(2) == 0 ? 1 : -1;
 
@@ -548,5 +527,5 @@
 	if (isSizeMod && geno.checkValidityOfPartSizes() != 0)
 	{
-		randomNode->modifiers[randomModifier]++;
+		randomNode->modifiers[randomModifier] = oldValue;
 		return false;
 	}
@@ -715,6 +694,6 @@
 	}
 
-	double min = std::max(minValues.at(key), valueAtMinVolume);
-	double max = std::min(maxValues.at(key), valueAtMaxVolume);
+	double min = std::max(Node::minValues.at(key), valueAtMinVolume);
+	double max = std::min(Node::maxValues.at(key), valueAtMaxVolume);
 
 	node->params[key] = GenoOperators::mutateCreep('f', node->getParam(key), min, max, true);
Index: cpp/frams/genetics/fS/fS_oper.h
===================================================================
--- cpp/frams/genetics/fS/fS_oper.h	(revision 1012)
+++ cpp/frams/genetics/fS/fS_oper.h	(revision 1017)
@@ -42,7 +42,4 @@
 	paInt useElli, useCub,  useCyl;
 	paInt strongAddPart;
-
-	std::map<string, double> minValues;
-	std::map<string, double> maxValues;
 
 	GenoOper_fS();
@@ -118,4 +115,13 @@
 
 	/**
+	 * Changes the value of specified parameter.
+	 * The state of the node must be previously calculated
+	 * @param node - the node on which parameter is modified
+	 * @param key - the key of parameter
+	 * @return
+	 */
+	bool mutateParamValue(Node *node, string key);
+
+	/**
 	 * Performs change modifier mutation on genotype
 	 * @return true if mutation succeeded, false otherwise
@@ -144,6 +150,4 @@
 	 */
 	bool mutateSizeParam(Node *node, string key, bool ensureCircleSection);
-
-	void prepareParams();
 };
 
Index: cpp/frams/genetics/fS/part_distance_estimator.h
===================================================================
--- cpp/frams/genetics/fS/part_distance_estimator.h	(revision 1012)
+++ cpp/frams/genetics/fS/part_distance_estimator.h	(revision 1017)
@@ -5,4 +5,6 @@
 #ifndef _PART_DISTANCE_ESTIMATOR_H_
 #define _PART_DISTANCE_ESTIMATOR_H_
+
+#include "frams/model/geometry/meshbuilder.h"
 
 class fS_Utils
@@ -44,146 +46,45 @@
 class PartDistanceEstimator
 {
-	/**
-	 * Used in finding the proper distance between the parts
-	 * distance between spheres / sphere radius
-	 * That default value can be changed in certain cases
-	 * */
-	static constexpr float SPHERE_RELATIVE_DISTANCE = 0.5;
-	/**
-	 * Used in finding the proper distance between the parts
-	 * The maximal allowed value for
-	 * maximal radius of the node / sphere radius
-	 */
-	static const int MAX_DIAMETER_QUOTIENT = 30;
-	/**
-	 * The tolerance of the value of distance between parts
-	 */
-	static constexpr double SPHERE_DISTANCE_TOLERANCE = 0.96;
 
-	static constexpr double SPHERE_SIZE_QUOTIENT = 0.25;
+public:
+	static constexpr double PRECISION = 0.05;
+	static constexpr double RELATIVE_DENSITY = 5.0;
 
 
-	static bool isInsidePart(Part::Shape shape, const Pt3D &partRadii, const Pt3D &center, double sphereRadius)
+	static Part *buildTemporaryPart(Part::Shape shape, const Pt3D &scale, const Pt3D &rotations)
 	{
-		if(sphereRadius >= fS_Utils::min3(partRadii))
-			return true; // TODO Special case
-		Pt3D mostRemote = Pt3D(fabs(center.x), fabs(center.y), fabs(center.z)) + Pt3D(sphereRadius);
-		bool isInEllipsis;
-
-		bool result = false;
-
-		switch (shape)
-		{
-			case Part::Shape::SHAPE_CUBOID:
-				if(mostRemote.x <= partRadii.x && mostRemote.y <= partRadii.y && mostRemote.z <= partRadii.z)
-					result = true;
-				break;
-			case Part::Shape::SHAPE_CYLINDER:
-				isInEllipsis = pow(center.y / (partRadii.y - sphereRadius), 2) + pow(center.z / (partRadii.z - sphereRadius), 2) <= 1.0;
-				if (mostRemote.x <= partRadii.x && isInEllipsis)
-					result = true;
-				break;
-			case Part::Shape::SHAPE_ELLIPSOID:
-				if(pow(center.x / (partRadii.x - sphereRadius), 2) + pow(center.y / (partRadii.y - sphereRadius), 2) + pow(center.z / (partRadii.z - sphereRadius), 2) <= 1.0)
-					result = true;
-				break;
-			default:
-				logMessage("fS", "calculateVolume", LOG_ERROR, "Invalid part type");
-		}
-		return result;
+		Part *tmpPart = new Part(shape);
+		tmpPart->scale = scale;
+		tmpPart->setRot(rotations);
+		return tmpPart;
 	}
 
-public:
-	static vector<Pt3D> findSphereCenters(Part::Shape shape, double &sphereRadius, const Pt3D &radii, const Pt3D &rotations)
+	static vector <Pt3D> findSphereCenters(Part *part)
 	{
-		double sphereRelativeDistance = SPHERE_RELATIVE_DISTANCE;
-		double minRadius = fS_Utils::min3(radii);
-		if(minRadius <= 0)
-			throw fS_Exception("Invalid part size in PartDistanceEstimator", 0);
-		double maxRadius = fS_Utils::max3(radii);
-		sphereRadius = SPHERE_SIZE_QUOTIENT * minRadius;
-		if (MAX_DIAMETER_QUOTIENT < maxRadius / sphereRadius)
+		// Divide by maximal radius to avoid long computations
+		MeshBuilder::PartSurface surface(RELATIVE_DENSITY / fS_Utils::max3(part->scale));
+		surface.initialize(part);
+
+		vector <Pt3D> centers;
+		Pt3D point;
+		while (surface.tryGetNext(point))
 		{
-			// When max radius is much bigger than sphereRadius and there are to many spheresZ
-			sphereRelativeDistance = 1.0;   // Make the spheres adjacent to speed up the computation
-			sphereRadius = maxRadius / MAX_DIAMETER_QUOTIENT;
-		}
-		else if(shape == Part::Shape::SHAPE_ELLIPSOID)
-			sphereRadius = minRadius;
-
-		double sphereDiameter = 2 * sphereRadius;
-
-		double radiiArr[3]{radii.x, radii.y, radii.z};
-		vector<double> coordinates[3];
-		for (int dim = 0; dim < 3; dim++)
-		{
-			double spaceForSphereCenters = 2 * radiiArr[dim] - sphereDiameter;
-			if (spaceForSphereCenters > 0)
-			{
-				int lastIndex = ceil(spaceForSphereCenters / (sphereDiameter * sphereRelativeDistance));
-				for (int i = 0; i <= lastIndex; i++)
-				{
-					coordinates[dim].push_back(spaceForSphereCenters * (double(i) / lastIndex - 0.5));
-				}
-			}
-			else
-				coordinates[dim].push_back(0.0);
-		}
-
-		vector<Pt3D> centers;
-		for (int xi = 0; xi < int(coordinates[0].size()); xi++)
-		{
-			for (int yi = 0; yi < int(coordinates[1].size()); yi++)
-			{
-				for (int zi = 0; zi < int(coordinates[2].size()); zi++)
-				{
-						Pt3D center = Pt3D(coordinates[0][xi], coordinates[1][yi], coordinates[2][zi]);
-						if (isInsidePart(shape, radii, center, sphereRadius))
-						{
-							fS_Utils::rotateVector(center, rotations);
-							centers.push_back(center);
-						}
-				}
-			}
+			centers.push_back(point);
 		}
 		return centers;
 	}
 
-	static int isCollision(vector<Pt3D> &centersParent, vector<Pt3D> &centers, Pt3D &vectorBetweenParts,
-					double distanceThreshold)
+	static bool isCollision(Part *parentPart, vector <Pt3D> &centers, Pt3D &vectorBetweenParts)
 	{
-		double upperThresholdSq = distanceThreshold * distanceThreshold;
-		double lowerThresholdSq = pow(SPHERE_DISTANCE_TOLERANCE * distanceThreshold, 2);
-		double distanceSq;
-		double dx, dy, dz;
-		bool existsAdjacent = false;
-		Pt3D *tmpPoint;
-		for (int sc = 0; sc < int(centers.size()); sc++)
+		static double CBRT_3 = std::cbrt(3);
+		double maxParentReachSq = pow(CBRT_3 * fS_Utils::max3(parentPart->scale), 2);
+		for (int i = 0; i < int(centers.size()); i++)
 		{
-			Pt3D shiftedSphere = Pt3D(centers[sc]);
-			shiftedSphere += vectorBetweenParts;
-			for (int psc = 0; psc < int(centersParent.size()); psc++)
-			{
-				tmpPoint = &centersParent[psc];
-				dx = shiftedSphere.x - tmpPoint->x;
-				dy = shiftedSphere.y - tmpPoint->y;
-				dz = shiftedSphere.z - tmpPoint->z;
-				distanceSq = dx * dx + dy * dy + dz * dz;
-
-				if (distanceSq <= upperThresholdSq)
-				{
-					if (distanceSq >= lowerThresholdSq)
-						existsAdjacent = true;
-					else
-					{
-						return COLLISION;
-					}
-				}
-			}
+			Pt3D shifted = centers[i] + vectorBetweenParts;
+			double distanceToCenterSq = shifted.x * shifted.x + shifted.y * shifted.y + shifted.z * shifted.z;
+			if (distanceToCenterSq <= maxParentReachSq && GeometryUtils::isPointInsidePart(shifted, parentPart))
+				return true;
 		}
-		if (existsAdjacent)
-			return ADJACENT;
-		else
-			return DISJOINT;
+		return false;
 	}
 };
@@ -191,43 +92,34 @@
 double Node::getDistance()
 {
-	Pt3D size = calculateSize();
-	Pt3D parentSize = parent->calculateSize();	// Here we are sure that parent is not nullptr
-	double parentSphereRadius, sphereRadius;
-	vector<Pt3D> centersParent = PartDistanceEstimator::findSphereCenters(parent->partType, parentSphereRadius, parentSize, parent->getRotation());
-	vector<Pt3D> centers = PartDistanceEstimator::findSphereCenters(partType, sphereRadius, size, getRotation());
+	Pt3D size;
+	calculateSize(size);
+	Pt3D parentSize;
+	parent->calculateSize(parentSize);    // Here we are sure that parent is not nullptr
+	Part *tmpPart = PartDistanceEstimator::buildTemporaryPart(partType, size, getRotation());
+	Part *parentTmpPart = PartDistanceEstimator::buildTemporaryPart(parent->partType, parentSize, parent->getRotation());
 
-	double distanceThreshold = sphereRadius + parentSphereRadius;
+	vector <Pt3D> centers = PartDistanceEstimator::findSphereCenters(tmpPart);
+
 	double minDistance = 0.0;
 	double maxDistance = 2 * (fS_Utils::max3(parentSize) + fS_Utils::max3(size));
 	double currentDistance = fS_Utils::avg(maxDistance, minDistance);
-	int result = -1;
-	int iterationNo = 0;
-	while (result != ADJACENT)
+	int collisionDetected = false;
+	while (maxDistance - minDistance > PartDistanceEstimator::PRECISION)
 	{
-		iterationNo++;
 		Pt3D vectorBetweenParts = state->v * currentDistance;
-		result = PartDistanceEstimator::isCollision(centersParent, centers, vectorBetweenParts, distanceThreshold);
+		collisionDetected = PartDistanceEstimator::isCollision(parentTmpPart, centers, vectorBetweenParts);
 
-		if (result == DISJOINT)
+		if (collisionDetected)
+		{
+			minDistance = currentDistance;
+			currentDistance = fS_Utils::avg(maxDistance, currentDistance);
+		} else
 		{
 			maxDistance = currentDistance;
 			currentDistance = fS_Utils::avg(currentDistance, minDistance);
-		} else if (result == COLLISION)
-		{
-			minDistance = currentDistance;
-			currentDistance = fS_Utils::avg(maxDistance, currentDistance);
 		}
-
-		if(maxDistance <= 0 || iterationNo > 1000)
-			throw fS_Exception("Computing of distances between parts failed", 0);
-		if (currentDistance > maxDistance)
-		{
-			throw fS_Exception("Internal error; the maximal distance between parts exceeded.", 0);
-		}
-		if (currentDistance < minDistance)
-			throw fS_Exception("Internal error; the minimal distance between parts exceeded.", 0);
-
 	}
-
+	delete tmpPart;
+	delete parentTmpPart;
 	return currentDistance;
 }
