CMesh.cpp

#include "CMesh.h"

Elem::Elem()
{
    trailUpper = NULL;
    trailLower = NULL;
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
//                     P U B L I C   M E M B E R S                           //
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

//-----------------------------------------------------------------------/
// basic constructor - initialize mesh data
//-----------------------------------------------------------------------/
CMesh::CMesh()
{
    initializeData();
}

//-----------------------------------------------------------------------/
// basic destructor - clear node and element maps
//-----------------------------------------------------------------------/
CMesh::~CMesh()
{
    clearData();
}

//-----------------------------------------------------------------------/
// get basic private member data
//-----------------------------------------------------------------------/
uInt CMesh::getNumBodyQuad() {return numBodyQuad;}
uInt CMesh::getNumDummQuad() {return numDummQuad;}
uInt CMesh::getNumWakeQuad() {return numWakeQuad;}
uInt CMesh::getNumBodyTria() {return numBodyTria;}
uInt CMesh::getNumDummTria() {return numDummTria;}
uInt CMesh::getNumWakeTria() {return numWakeTria;}
uInt CMesh::getNumNode()     {return numNode    ;}
uInt CMesh::getMaxNodeID()   {return maxNodeId  ;}
uInt CMesh::getMaxElemID()   {return maxElemId  ;}
double* CMesh::getMaxNode()  {return maxNode    ;}
double* CMesh::getMinNode()  {return minNode    ;}
double* CMesh::getMaxBody()  {return maxBody    ;}
double* CMesh::getMinBody()  {return minBody    ;}
double* CMesh::getMaxWake()  {return maxWake    ;}
double* CMesh::getMinWake()  {return minWake    ;}
bool CMesh::checkIfNodeExists(const uInt& nodeID) {return meshNodes.find(nodeID)!=meshNodes.end();}
bool CMesh::checkIfElemExists(const uInt& elemID) {return meshElems.find(elemID)!=meshElems.end();}

//-----------------------------------------------------------------------/
// get coordinates of node
//-----------------------------------------------------------------------/
int CMesh::getNodeCoords(
        const uInt& nodeID,
        double coordinates[3])
{
    // get iterator to node with given ID
    mapNodeConstIt nodeIt = meshNodes.find(nodeID);

    // check that it exists
    if (nodeIt == meshNodes.end())
        return -1;

    // extract coordinates
    for (uInt i=0; i<3; i++)
        coordinates[i] = nodeIt->second.coordinates[i];

    return 0;
}

//-----------------------------------------------------------------------/
// extract node data
// data is used for writing ApameSolver input file and for OpenGL drawing
//-----------------------------------------------------------------------/
void CMesh::getNodeData(vecUInt& nodeIDs, vecVecDbl& coordinates)
{
    // resize vectors to number of nodes
    nodeIDs.resize(numNode);
    coordinates.resize(numNode);

    mapNodeConstIt nodeIt;
    uInt nodeIndex;
    for (nodeIt=meshNodes.begin(), nodeIndex=0; nodeIt!=meshNodes.end(); nodeIt++, nodeIndex++){
        // extract node ID
        nodeIDs[nodeIndex] = nodeIt->second.id;

        // extract coordinates
        coordinates[nodeIndex].resize(3);
        for (uInt coordIndex=0; coordIndex<3; coordIndex++)
            coordinates[nodeIndex][coordIndex] = nodeIt->second.coordinates[coordIndex];
    }
}

//-----------------------------------------------------------------------/
// extract element data
// data is used for filling edit element dialog
//-----------------------------------------------------------------------/
int CMesh::getElemData(
        const uInt& elemID,
        uInt& elemType,
        vecUInt& nodeIDs,
        vecUInt& neighborTrailIDs)
{
    // get iterator to element with given ID
    mapElemConstIt elemIt = meshElems.find(elemID);

    // check that it exists
    if (elemIt == meshElems.end())
        return -1;

    // extract element type
    elemType = elemIt->second.type;

    // extract node IDs
    nodeIDs.resize(elemIt->second.nodes.size());
    for (uInt i=0; i<nodeIDs.size(); i++)
        nodeIDs[i] = elemIt->second.nodes[i]->id;

    // check type
    if (elemType==1 || elemType==2){
        // extract neighbors
        neighborTrailIDs.resize(elemIt->second.neighbors.size());
        uInt i = 0;
        for (setPElemConstIt neighIt=elemIt->second.neighbors.begin(); neighIt!=elemIt->second.neighbors.end(); neighIt++)
            neighborTrailIDs[i++] = (*neighIt)->id;
    }
    else if (elemType==10 || elemType==11){
        // extract trailing elements
        if (elemIt->second.trailUpper)
            neighborTrailIDs.push_back(elemIt->second.trailUpper->id);
        if (elemIt->second.trailLower)
            neighborTrailIDs.push_back(elemIt->second.trailLower->id);
    }
    else
        neighborTrailIDs.clear();

    return 0;
}

//-----------------------------------------------------------------------/
// get coordinates of element
//-----------------------------------------------------------------------/
int CMesh::getElemCoords(
        const uInt& elemID,
        vecVecDbl& coordinates)
{
    // get iterator to element with given ID
    mapElemConstIt elemIt = meshElems.find(elemID);

    // check that it exists
    if (elemIt == meshElems.end())
        return -1;

    // extract coordinates
    coordinates.resize(elemIt->second.nodes.size());
    for (uInt nodeIndex=0; nodeIndex<coordinates.size(); nodeIndex++){
        coordinates[nodeIndex].resize(3);
        for (uInt i=0; i<3; i++)
            coordinates[nodeIndex][i] = elemIt->second.nodes[nodeIndex]->coordinates[i];
    }

    return 0;
}

//-----------------------------------------------------------------------/
// extract element data
//-----------------------------------------------------------------------/
int CMesh::getElemsData(
        vecUInt& elemIDs,
        vecUInt& elemTypes,
        vecVecUInt& nodeIDs)
{
    // get number of elements
    uInt elemNum = meshElems.size();

    // resize vectors to number of elements
    elemIDs.resize(elemNum);
    elemTypes.resize(elemNum);
    nodeIDs.resize(elemNum);

    // for each element
    mapElemConstIt elemIt;
    uInt elemIndex;
    for (elemIt=meshElems.begin(), elemIndex=0; elemIt!=meshElems.end(); elemIt++, elemIndex++){
        // extract element ID
        elemIDs[elemIndex] = elemIt->second.id;

        // extract element type
        elemTypes[elemIndex] = elemIt->second.type;

        // extract element nodes
        nodeIDs[elemIndex].resize(elemIt->second.nodes.size());
        for (uInt i=0; i<nodeIDs[elemIndex].size(); i++)
            nodeIDs[elemIndex][i] = elemIt->second.nodes[i]->id;
    }

    return 0;
}

//-----------------------------------------------------------------------/
// extract element data for writing ApameSolver input file
//-----------------------------------------------------------------------/
int CMesh::getElemIDsTypes(vecUInt& elemIDs, vecUInt& elemTypes)
{
    // get number of elements
    uInt elemNum = meshElems.size();

    // resize vectors to number of elements
    elemIDs.resize(elemNum);
    elemTypes.resize(elemNum);

    // for each element
    mapElemConstIt elemIt;
    uInt elemIndex;
    for (elemIt=meshElems.begin(), elemIndex=0; elemIt!=meshElems.end(); elemIt++, elemIndex++){
        // extract element ID
        elemIDs[elemIndex] = elemIt->second.id;

        // extract element type
        elemTypes[elemIndex] = elemIt->second.type;
    }

    return 0;
}

//-----------------------------------------------------------------------/
// extract element data for writing ApameSolver input file
//-----------------------------------------------------------------------/
int CMesh::getElemSolverData(vecUInt& elemIDs, vecUInt& elemTypes, vecVecUInt& elemNodeIDs, vecVecUInt& neighborIDs, vecVecUInt& trailIDs)
{
    // get number of elements
    uInt elemNum = meshElems.size();

    // resize vectors to number of elements
    elemIDs.resize(elemNum);
    elemTypes.resize(elemNum);
    elemNodeIDs.resize(elemNum);
    neighborIDs.resize(elemNum);
    trailIDs.resize(elemNum);

    // for each element
    bool neighWarning = false;
    mapElemConstIt elemIt;
    uInt elemIndex;
    for (elemIt=meshElems.begin(), elemIndex=0; elemIt!=meshElems.end(); elemIt++, elemIndex++){
        // extract element ID
        elemIDs[elemIndex] = elemIt->second.id;

        // extract element type
        elemTypes[elemIndex] = elemIt->second.type;

        // extract node IDs
        elemNodeIDs[elemIndex].resize(elemIt->second.nodes.size());
        for (uInt nodeIndex=0; nodeIndex<elemIt->second.nodes.size(); nodeIndex++)
            elemNodeIDs[elemIndex][nodeIndex] = elemIt->second.nodes[nodeIndex]->id;

        // extract neighbor IDs
        if (elemIt->second.type==1 || elemIt->second.type==2){
            if (elemIt->second.neighbors.empty())
                neighWarning = true;
            neighborIDs[elemIndex].resize(elemIt->second.neighbors.size());

            // go through neighbors and add by ID
            uInt neighborIndex = 0;
            for (setPElemConstIt neighIt=elemIt->second.neighbors.begin(); neighIt!=elemIt->second.neighbors.end(); neighIt++)
                neighborIDs[elemIndex][neighborIndex++] = (*neighIt)->id;
        }

        // extract trailing elements for wake elements
        if (elemIt->second.type==10 || elemIt->second.type==11){
            if (!elemIt->second.trailUpper || !elemIt->second.trailLower)
                return -1;
            trailIDs[elemIndex].resize(2);
            trailIDs[elemIndex][0] = elemIt->second.trailUpper->id;
            trailIDs[elemIndex][1] = elemIt->second.trailLower->id;
        }
    }

    if (neighWarning)
        return 1;

    return 0;
}

//-----------------------------------------------------------------------/
// extract element data for drawing
//-----------------------------------------------------------------------/
void CMesh::getElemIDsTypesCoords(vecUInt& ids, vecUInt& types, vecVecVecDbl& coordinates)
{
    // resize vectors to number of elements
    uInt elemNum = meshElems.size();
    ids.resize(elemNum);
    types.resize(elemNum);
    coordinates.resize(elemNum);

    mapElemConstIt elemIt;
    uInt elemIndex;
    for (elemIt=meshElems.begin(), elemIndex=0; elemIt!=meshElems.end(); elemIt++, elemIndex++){
        // extract element ID
        ids[elemIndex] = elemIt->second.id;

        // extract element type
        types[elemIndex] = elemIt->second.type;

        // extract constituing node coordinates
        coordinates[elemIndex].resize(elemIt->second.nodes.size());
        for (uInt nodeIndex=0; nodeIndex<elemIt->second.nodes.size(); nodeIndex++){
            coordinates[elemIndex][nodeIndex].resize(3);
            for (uInt coordIndex=0; coordIndex<3; coordIndex++)
                coordinates[elemIndex][nodeIndex][coordIndex] = elemIt->second.nodes[nodeIndex]->coordinates[coordIndex];
        }
    }
}

//-----------------------------------------------------------------------/
// extract neighbors for neighbors display mode
//-----------------------------------------------------------------------/
void CMesh::getNeighbors(
        const uInt& monitoredElement,
        vecVecUInt& neighboring)
{
    neighboring.resize(4);
    setUInt usedIDs; // set used not to store already used elements
    // try to find monitored element
    mapElemConstIt elemIt = meshElems.find(monitoredElement);
    if (elemIt != meshElems.end()){
        // element exists so store it as first order
        neighboring[0].push_back(monitoredElement);
        usedIDs.insert(monitoredElement);
        
        // find its neighbors (second order)
        for (setPElemConstIt neighIt=elemIt->second.neighbors.begin(); neighIt!=elemIt->second.neighbors.end(); neighIt++){
            neighboring[1].push_back((*neighIt)->id);
            usedIDs.insert((*neighIt)->id);
        }
        
        // find second neighbors (third order)
        for (uInt i=0; i<neighboring[1].size(); i++){
            mapElemConstIt secElemIt = meshElems.find(neighboring[1][i]);
            if (secElemIt != meshElems.end()){
                for (setPElemConstIt neighIt=secElemIt->second.neighbors.begin(); neighIt!=secElemIt->second.neighbors.end(); neighIt++){
                    if (usedIDs.find((*neighIt)->id) == usedIDs.end())
                        neighboring[2].push_back((*neighIt)->id);
                }
            }
        }
        
        // find third neighbors (fourth order)
        for (uInt i=0; i<neighboring[2].size(); i++){
            mapElemConstIt thirdElemIt = meshElems.find(neighboring[2][i]);
            if (thirdElemIt != meshElems.end()){
                for (setPElemConstIt neighIt=thirdElemIt->second.neighbors.begin(); neighIt!=thirdElemIt->second.neighbors.end(); neighIt++){
                    if (usedIDs.find((*neighIt)->id) == usedIDs.end())
                        neighboring[3].push_back((*neighIt)->id);
                }
            }
        }
        
        // add wakes to second order elements
        for (setPElemConstIt neighIt=elemIt->second.wakes.begin(); neighIt!=elemIt->second.wakes.end(); neighIt++)
            neighboring[1].push_back((*neighIt)->id);
        
        // add upper trail to second order elements
        if (elemIt->second.trailUpper)
            neighboring[1].push_back(elemIt->second.trailUpper->id);
        
        // add lower trail to third order elements
        if (elemIt->second.trailLower)
            neighboring[2].push_back(elemIt->second.trailLower->id);
    }
}

//-----------------------------------------------------------------------/
// extract mesh data for IDW
//-----------------------------------------------------------------------/
void CMesh::getIDWData(
        vecUInt& nodeIDs,
        vecVecDbl& coordinates,
        vecUInt& elemIDs,
        vecVecUInt& elemNodeIDs)
{
    // resize node vectors to number of nodes
    nodeIDs.resize(numNode);
    coordinates.resize(numNode);

    // for each node
    mapNodeConstIt nodeIt;
    uInt nodeIndex;
    for (nodeIt=meshNodes.begin(), nodeIndex=0; nodeIt!=meshNodes.end(); nodeIt++, nodeIndex++){
        // extract node ID
        nodeIDs[nodeIndex] = nodeIt->second.id;

        // extract coordinates
        coordinates[nodeIndex].resize(3);
        for (uInt coordIndex=0; coordIndex<3; coordIndex++)
            coordinates[nodeIndex][coordIndex] = nodeIt->second.coordinates[coordIndex];
    }

    // resize element vectors to number of elements
    elemIDs.resize(meshElems.size());
    elemNodeIDs.resize(meshElems.size());

    // for each element
    mapElemConstIt elemIt;
    uInt elemIndex;
    for (elemIt=meshElems.begin(), elemIndex=0; elemIt!=meshElems.end(); elemIt++, elemIndex++){
        // extract element ID
        elemIDs[elemIndex] = elemIt->second.id;

        // extract node IDs
        elemNodeIDs[elemIndex].resize(elemIt->second.nodes.size());
        for (uInt nodeIndex=0; nodeIndex<elemIt->second.nodes.size(); nodeIndex++)
            elemNodeIDs[elemIndex][nodeIndex] = elemIt->second.nodes[nodeIndex]->id;
    }
}

//-----------------------------------------------------------------------/
// insert node with defined coordinates
// returned value:
//  node ID
//-----------------------------------------------------------------------/
int CMesh::insertNodeWithCoords(const double coordinates[3])
{
    // insert node with private function
    Node* newNodeP = insertNodeWithCoordsPrivate(coordinates);

    // return node ID
    return newNodeP->id;
}

//-----------------------------------------------------------------------/
// insert node defined with ID and coordinates
// returned error codes:
//   0 = OK
//  -1 = given ID already exists
//  -2 = given ID is equal 0
//-----------------------------------------------------------------------/
int CMesh::insertNodeWithIDCoords(const uInt& id, const double coordinates[3])
{
    // check if node ID already exist
    if (meshNodes.find(id) != meshNodes.end())
        return -1;

    // check if node ID is equal to zero
    if (id == 0)
        return -2;

    // add new node with given ID to node map
    Node* newNodeP = &meshNodes[id] ;

    // add node ID
    newNodeP->id = id;
    // add coordinates
    for (uInt i=0; i<3; i++)
        newNodeP->coordinates[i] = coordinates[i];

    // check if node ID is greatest yet
    if (id > maxNodeId)
        maxNodeId = id;

    // update bounds
    updateBoundsWithNode(newNodeP);

    // increase node counter
    numNode++;

    return 0;
}

//-----------------------------------------------------------------------/
// insert element with type and coordinates of constituing nodes
//
// returned codes:
//  >0 = element ID
//  -1 = number of coordinates per node is not 3
//  -2 = zero surface element
//  -3 = incorrect element type
//  -4 = number of nodes is not 3 or 4
//-----------------------------------------------------------------------/
int CMesh::insertElemWithNodeCoords(const vecVecDbl& coordinatesAll, const uInt& type)
{
    Elem* newElemP = NULL;
    return insertElemWithNodeCoordsPrivate(coordinatesAll, type, newElemP);
}

//-----------------------------------------------------------------------/
// insert element with type and constituing nodes ID's
//  >0 = element ID
//  -1 = incorrect number of nodes given
//  -2 = zero surface element
//  -3 = incorrect element type
//  -4 = duplicate node ID(s) given
//  -5 = node(s) does not exist
//  -6 = type does not conform to number of nodes
//-----------------------------------------------------------------------/
int CMesh::insertElemWithNodeIDs(const uInt& type, const vecUInt& nodeIDs)
{
    // check if incorrect type is given
    if (type != 1  &&
        type != 2  &&
        type != 10 &&
        type != 11 &&
        type != 20 &&
        type != 21 )
        return -3;

    // check if proper node number is given
    if (nodeIDs.size()!=3 && nodeIDs.size()!=4)
        return -1;

    // check if type does not conform to number of nodes
    if (nodeIDs.size()==3 && (type==1 || type==10 || type==20))
        return -6;

    // check if duplicate node ID's are given
    if (nodeIDs.size() == 4){
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[3] ||
            nodeIDs[3] == nodeIDs[0] ||
            nodeIDs[0] == nodeIDs[2] ||
            nodeIDs[1] == nodeIDs[3] )
            return -4;
    }
    else if (nodeIDs.size() == 3){
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[0] )
            return -4;
    }

    // check if nodes exist and assemble their pointers
    vecPNode nodes(nodeIDs.size());
    for (uInt i=0; i<nodeIDs.size(); i++){
        mapNodeIt it = meshNodes.find(nodeIDs[i]);
        if (it == meshNodes.end())
            return -5;
        else
            nodes[i] = &it->second;
    }

    // check normal
    if (calcNormal(nodes) == -1)
        return -2;

    // get new element ID
    uInt id = ++maxElemId;

    // add new element with given ID to element map
    Elem* newElemP = &meshElems[id];

    // add element ID
    newElemP->id = id;
    // add element type
    newElemP->type = type;
    // add constituing nodes (pointers)
    newElemP->nodes = nodes;
    // add element associativity to nodes
    for (uInt i=0; i<nodeIDs.size(); i++)
        newElemP->nodes[i]->elements.insert(newElemP);

    // update bounds with element
    updateBoundsWithElem(newElemP);

    // add to counter of elements
    if      (type == 1 ) numBodyQuad++;
    else if (type == 2 ) numBodyTria++;
    else if (type == 10) numWakeQuad++;
    else if (type == 11) numWakeTria++;
    else if (type == 20) numDummQuad++;
    else if (type == 21) numDummTria++;

    return id;
}

//-----------------------------------------------------------------------/
// insert element with ID and constituing nodes ID's
//   0 = OK
//  -1 = element ID already exists
//  -2 = zero element ID
//  -3 = incorrect number of nodes given
//  -4 = duplicate node ID(s) given
//  -5 = node(s) does not exist
//  -6 = zero surface element
//-----------------------------------------------------------------------/
int CMesh::insertElemWithIDNodeIDs(const uInt& id, const vecUInt& nodeIDs)
{
    // check if element ID already exists
    if (meshElems.find(id) != meshElems.end())
        return -1;

    // check if element ID equals zero
    if (id == 0)
        return -2;

    // check if duplicate node ID's are given
    uInt type;
    if (nodeIDs.size() == 4){
        type = 1;
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[3] ||
            nodeIDs[3] == nodeIDs[0] ||
            nodeIDs[0] == nodeIDs[2] ||
            nodeIDs[1] == nodeIDs[3] )
            return -4;
    }
    else if (nodeIDs.size() == 3){
        type = 2;
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[0] )
            return -4;
    }
    else
        return -3;

    // check if nodes exist and assemble their pointers
    vecPNode nodes(nodeIDs.size());
    for (uInt i=0; i<nodeIDs.size(); i++){
        mapNodeIt it = meshNodes.find(nodeIDs[i]);
        if (it == meshNodes.end())
            return -5;
        else
            nodes[i] = &it->second;
    }

    // check normal
    if (calcNormal(nodes) == -1)
        return -6;

    // update maximum element ID
    if (id > maxElemId)
        maxElemId = id;

    // add new element with given ID to element map
    Elem* newElemP = &meshElems[id];

    // add element ID
    newElemP->id = id;
    // add element type
    newElemP->type = type;
    // add constituing nodes (pointers)
    newElemP->nodes = nodes;
    // add element associativity to node
    for (uInt i=0; i<nodeIDs.size(); i++)
        newElemP->nodes[i]->elements.insert(newElemP);

    // update bounds with element
    updateBoundsWithElem(newElemP);

    // add to counter of elements
    if      (type == 1) numBodyQuad++;
    else if (type == 2) numBodyTria++;

    return 0;
}

//-----------------------------------------------------------------------/
// insert element with ID and constituing nodes ID's
//   0 = OK
//  -1 = element ID already exists
//  -2 = zero element ID
//  -3 = incorrect number of nodes given
//  -4 = duplicate node ID(s) given
//  -5 = node(s) does not exist
//  -6 = zero surface element
//  -7 = type does not match number of nodes
//  -8 = incorrect type given
//-----------------------------------------------------------------------/
int CMesh::insertElemWithIDTypeNodeIDs(
        const uInt& id,
        const uInt& type,
        const vecUInt& nodeIDs)
{
    // check if element ID already exists
    if (meshElems.find(id) != meshElems.end())
        return -1;

    // check if element ID equals zero
    if (id == 0)
        return -2;

    // check for proper element type
    if (type!=1 && type!=2 && type!=10 && type!=11 && type!=20 && type!=21)
        return -8;

    // check if proper node number is given
    if (nodeIDs.size()!=3 && nodeIDs.size()!=4)
        return -3;

    // check if duplicate node ID's are given
    if (nodeIDs.size() == 4){
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[3] ||
            nodeIDs[3] == nodeIDs[0] ||
            nodeIDs[0] == nodeIDs[2] ||
            nodeIDs[1] == nodeIDs[3] )
            return -4;
    }
    else if (nodeIDs.size() == 3){
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[0] )
            return -4;
    }

    // check if type matches number of nodes
    if (nodeIDs.size()==4 && (type==2 || type==11 || type==21))
        return -7;

    // check if nodes exist and assemble their pointers
    vecPNode nodes(nodeIDs.size());
    for (uInt i=0; i<nodeIDs.size(); i++){
        mapNodeIt it = meshNodes.find(nodeIDs[i]);
        if (it == meshNodes.end())
            return -5;
        else
            nodes[i] = &it->second;
    }

    // check normal
    if (calcNormal(nodes) == -1)
        return -6;

    // update maximum element ID
    if (id > maxElemId)
        maxElemId = id;

    // add new element with given ID to element map
    Elem* newElemP = &meshElems[id];

    // add element ID
    newElemP->id = id;
    // add element type
    newElemP->type = type;
    // add constituing nodes (pointers)
    newElemP->nodes = nodes;
    // add element associativity to node
    for (uInt i=0; i<nodeIDs.size(); i++)
        newElemP->nodes[i]->elements.insert(newElemP);

    // update bounds with element
    updateBoundsWithElem(newElemP);

    // add to counter of elements
    if      (type == 1 ) numBodyQuad++;
    else if (type == 2 ) numBodyTria++;
    else if (type == 10) numWakeQuad++;
    else if (type == 11) numWakeTria++;
    else if (type == 20) numDummQuad++;
    else if (type == 21) numDummTria++;

    return 0;
}

//-----------------------------------------------------------------------/
// insert neighboring information to element
//   0 = OK
//  -1 = element does not exist
//  -2 = neighbor(s) does not exist
//  -3 = same ID in element and neighbor
//  -4 = trying to assign neighbors to wake element
//  -5 = cannot assign neighbors to wake element
//-----------------------------------------------------------------------/
int CMesh::insertNeighborsToElem(const uInt& id, const vecUInt& neighbors)
{
    // check if element ID exists
    mapElemIt elemIt = meshElems.find(id);
    if (elemIt == meshElems.end())
        return -1;

    // extract element
    Elem* element = &elemIt->second;

    // check if wake type
    if (element->type==10 || element->type==11)
        return -5;

    // check neighbors
    vecPElem pNeighbors(neighbors.size());
    for (uInt i=0; i<neighbors.size(); i++){
        // check that exists
        mapElemIt neighIt = meshElems.find(neighbors[i]);
        if (neighIt == meshElems.end())
            return -2;
        pNeighbors[i] = &neighIt->second;
        // check that neighbor is not same as watched element
        if (pNeighbors[i]->id == element->id)
            return -3;
    }

    // add neighboring associativity
    for (uInt i=0; i<pNeighbors.size(); i++){
        // neighbors to element
        element->neighbors.insert(pNeighbors[i]);
        // element to neighbors
        pNeighbors[i]->neighbors.insert(element);
    }

    return 0;
}

//-----------------------------------------------------------------------/
// insert trailing elements information to wake elements
//   0 = OK
//  -1 = element does not exist
//  -2 = trail(s) does not exist
//  -3 = same ID in element and trail
//  -4 = non-wake element
//  -5 = more that two trails given
//-----------------------------------------------------------------------/
int CMesh::insertTrailsToElem(const uInt& id, const vecUInt& trailIDs)
{
    // check if element ID exists
    mapElemIt elemIt = meshElems.find(id);
    if (elemIt == meshElems.end())
        return -1;

    // extract element
    Elem* element = &elemIt->second;

    // check if not wake type
    if (element->type!=10 && element->type!=11)
        return -4;

    // check number of trails
    if (trailIDs.size() > 2)
        return -5;

    // check trails
    vecPElem pTrails(trailIDs.size());
    for (uInt i=0; i<trailIDs.size(); i++){
        // check that exists
        mapElemIt trailIt = meshElems.find(trailIDs[i]);
        if (trailIt == meshElems.end())
            return -2;
        pTrails[i] = &trailIt->second;
        // check that trail is not same as watched element
        if (pTrails[i]->id == element->id)
            return -3;
    }

    // add trailing associativity to current element (wake)
    if (pTrails.size() == 2){
        element->trailUpper = pTrails[0];
        element->trailLower = pTrails[1];
    }
    else if (pTrails.size() == 1)
        element->trailUpper = pTrails[0];

    // add wake associativity to trailing elements
    for (uInt i=0; i<pTrails.size(); i++)
        pTrails[i]->wakes.insert(element);

    return 0;
}

//-----------------------------------------------------------------------/
// edit node
//-----------------------------------------------------------------------/
int CMesh::editNode(
        const uInt& nodeID,
        const uInt& newNodeID,
        const double newCoordinates[3])
{
    // check if node exists and get its iterator
    mapNodeIt nodeIt = meshNodes.find(nodeID);
    if (nodeIt == meshNodes.end())
        return -1;

    // extract node
    Node* node = &(nodeIt->second);

    // check if change of node ID is to happen
    if (nodeID == newNodeID){
        // change only coordinates
        for (uInt i=0; i<3; i++)
            node->coordinates[i] = newCoordinates[i];
    }
    else{
        // check if new node exists
        if (meshNodes.find(newNodeID) != meshNodes.end())
            return -2;

        // get node elements
        setPElem elements = node->elements;

        // erase node from map
        meshNodes.erase(node->id);

        // insert new node
        Node* newNodeP;
        newNodeP = &meshNodes[newNodeID];

        // assign node ID
        newNodeP->id = newNodeID;

        // assign coordinates
        for (uInt i=0; i<3; i++)
            newNodeP->coordinates[i] = newCoordinates[i];

        // assign elements
        newNodeP->elements = elements;

        // update node references in elements
        for (setPElemIt elemIt=newNodeP->elements.begin(); elemIt!=newNodeP->elements.end(); elemIt++){
            for (uInt i=0; i<(*elemIt)->nodes.size(); i++){
                if ((*elemIt)->nodes[i] == node)
                    (*elemIt)->nodes[i] = newNodeP;
            }
        }

        // update bounds
        updateNodeBounds();

        // update maximum node ID
        updateMaxNodeID();
    }
    return 0;
}

//-----------------------------------------------------------------------/
// edit element
//-----------------------------------------------------------------------/
int CMesh::editElement(
        const uInt& elemID,
        const uInt& newElemID,
        const uInt& type,
        const vecUInt& nodeIDs,
        const vecUInt& neighborTrailIDs)
{
    // check if element exists
    if (meshElems.find(elemID) == meshElems.end())
        return -1;

    // check if element with new ID exists
    if (newElemID!=elemID && meshElems.find(newElemID)!=meshElems.end())
        return -2;

    // check element type
    if (type!=1 && type!=2 && type!=10 && type!=11 && type!=20 && type!=21)
        return -3;

    // check if duplicate node IDs are given
    if (nodeIDs.size() == 4){
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[3] ||
            nodeIDs[3] == nodeIDs[0] ||
            nodeIDs[0] == nodeIDs[2] ||
            nodeIDs[1] == nodeIDs[3] )
            return -4;
    }
    else if (nodeIDs.size() == 3){
        if (nodeIDs[0] == nodeIDs[1] ||
            nodeIDs[1] == nodeIDs[2] ||
            nodeIDs[2] == nodeIDs[0] )
            return -4;
    }
    else
        return -5;

    // check if nodes exist and assemble their pointers
    vecPNode nodes(nodeIDs.size());
    for (uInt i=0; i<nodeIDs.size(); i++){
        mapNodeIt it = meshNodes.find(nodeIDs[i]);
        if (it == meshNodes.end())
            return -6;
        else
            nodes[i] = &it->second;
    }

    // check normal
    if (calcNormal(nodes))
        return -7;

    // check if duplicate element IDs are given
    if (neighborTrailIDs.size() == 4){
        if (neighborTrailIDs[0] == neighborTrailIDs[1] ||
            neighborTrailIDs[1] == neighborTrailIDs[2] ||
            neighborTrailIDs[2] == neighborTrailIDs[3] ||
            neighborTrailIDs[3] == neighborTrailIDs[0] ||
            neighborTrailIDs[0] == neighborTrailIDs[2] ||
            neighborTrailIDs[1] == neighborTrailIDs[3] )
            return -8;
    }
    else if (neighborTrailIDs.size() == 3){
        if (neighborTrailIDs[0] == neighborTrailIDs[1] ||
            neighborTrailIDs[1] == neighborTrailIDs[2] ||
            neighborTrailIDs[2] == neighborTrailIDs[0] )
            return -8;
    }
    else if (neighborTrailIDs.size() == 2){
        if (neighborTrailIDs[0] == neighborTrailIDs[1])
            return -8;
    }

    // check that neighboring/trailing elements exist
    for (uInt i=0; i<neighborTrailIDs.size(); i++){
        if (meshElems.find(neighborTrailIDs[i]) == meshElems.end())
            return -9;
    }

    // delete old element
    deleteElement(meshElems.find(elemID));

    // insert new element
    forceInsertElemWithIDNodes(newElemID, type, nodes);

    // add neighboring/trailing information
    if      (type==1  || type==2 ) insertNeighborsToElem(newElemID, neighborTrailIDs);
    else if (type==10 || type==11) insertTrailsToElem   (newElemID, neighborTrailIDs);

    return 0;
}

//-----------------------------------------------------------------------/
// inverse normals of given elements
//-----------------------------------------------------------------------/
void CMesh::reverseNormal(const setUInt& elemSelection)
{
    for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
        mapElemIt elemIt = meshElems.find(*it);
        if (elemIt == meshElems.end())
            continue;

        // extract element
        Elem* element = &(elemIt->second);

        // substitute first and third node
        Node* tempNode = element->nodes[0];
        element->nodes[0] = element->nodes[2];
        element->nodes[2] = tempNode;
    }
}

//-----------------------------------------------------------------------/
// merge nodes
//-----------------------------------------------------------------------/
void CMesh::mergeNodes(
        const double& tolerance,
        const string& type,
        const setUInt& nodeSelection)
{
    // assemble node pointers vector
    vecPNode nodePointers;
    uInt nodeIndex = 0;
    if (nodeSelection.empty()){
        nodePointers.resize(numNode);
        for (mapNodeIt nodeIt=meshNodes.begin(); nodeIt!=meshNodes.end(); nodeIt++)
            nodePointers[nodeIndex++] = &nodeIt->second;
    }
    else{
        for (setConstUIntIt it=nodeSelection.begin(); it!=nodeSelection.end(); it++){
            mapNodeIt nodeIt = meshNodes.find(*it);
            if (nodeIt != meshNodes.end())
                nodePointers.push_back(&nodeIt->second);
        }
    }

    // assemble distance matrix
    vecVecDbl distanceMatrix(nodePointers.size());
    for (uInt i=0; i<nodePointers.size(); i++){
        distanceMatrix[i].resize(nodePointers.size());
        for (uInt j=i+1; j<nodePointers.size(); j++){ // notice i+1 in order to skip symmetric part of matrix
            distanceMatrix[i][j] = sqrt( pow(nodePointers[j]->coordinates[0] - nodePointers[i]->coordinates[0], 2) +
                                         pow(nodePointers[j]->coordinates[1] - nodePointers[i]->coordinates[1], 2) +
                                         pow(nodePointers[j]->coordinates[2] - nodePointers[i]->coordinates[2], 2) );
        }
    }

    // perform merging across node matrix
    for (uInt i=0; i<nodePointers.size(); i++){
        for (uInt j=i+1; j<nodePointers.size(); j++){
            // check if already merged
            if (nodePointers[i] == nodePointers[j])
                continue;

            // check tolerance
            if (distanceMatrix[i][j] > tolerance)
                continue;

            // check merge type
            if      (type == "higher"){
                // merge to higher node ID
                if (nodePointers[i]->id > nodePointers[j]->id) mergeNodePair(nodePointers[i], nodePointers[j]);
                else                                           mergeNodePair(nodePointers[j], nodePointers[i]);
            }
            else if (type == "lower"){
                // merge to lower node ID
                if (nodePointers[i]->id > nodePointers[j]->id) mergeNodePair(nodePointers[j], nodePointers[i]);
                else                                           mergeNodePair(nodePointers[i], nodePointers[j]);
            }
            else // merge to middle
                mergeNodePair(nodePointers[i], nodePointers[j], true);
        }
    }
}

//-----------------------------------------------------------------------/
// merge pair of nodes
//-----------------------------------------------------------------------/
void CMesh::mergeNodePair(
        Node* mastr,
        Node* slave,
        bool toMiddle)
{
    // edit master node coordinates if merge to middle requested
    if (toMiddle){
        for (uInt i=0; i<3; i++)
            mastr->coordinates[i] = (mastr->coordinates[i] + slave->coordinates[i])/2;
    }

    // for each element of slave node
    for (setPElemIt elemIt=slave->elements.begin(); elemIt!=slave->elements.end(); elemIt++){
        // extract element
        Elem* element = *elemIt;

        // find reference to slave node and replace it with master node
        for (uInt nodeIndex=0; nodeIndex<element->nodes.size(); nodeIndex++){
            if (element->nodes[nodeIndex]->id == slave->id){
                // remove slave node reference from that element
                element->nodes.erase(element->nodes.begin()+nodeIndex);

                // add master node reference to that element
                element->nodes.insert(element->nodes.begin()+nodeIndex, mastr);

                // add that element to master node
                mastr->elements.insert(element);

                // no need to search for more nodes as element can have only one (slave) node
                break;
            }
        }
    }

    // remove element references from slave node
    slave->elements.clear();

    // delete slave node
    deleteNode(meshNodes.find(slave->id));

    // update bounds
    updateNodeBounds();

    // update maximum node ID
    updateMaxNodeID();

    // declare slave as master since they are now merged (to be able to merge more than two nodes into one later on)
    slave = mastr;
}

//-----------------------------------------------------------------------/
// scale (selection) by ratio and origin
//-----------------------------------------------------------------------/
void CMesh::scale(
        const double ratio[3],
        const double origin[3],
        const setUInt& nodeSelection,
        const setUInt& elemSelection)
{
    if (!nodeSelection.empty()){
        // scale node selection
        for (setConstUIntIt it=nodeSelection.begin(); it!=nodeSelection.end(); it++){
            mapNodeIt nodeIt = meshNodes.find(*it);
            if (nodeIt == meshNodes.end())
                continue;

            Node* node = &(nodeIt->second);
            for (uInt i=0; i<3; i++)
                node->coordinates[i] = origin[i]-(origin[i]-node->coordinates[i])*ratio[i];
        }
    }
    else if (!elemSelection.empty()){
        // scale element selection
        // extract nodes of all elements into one set (thus avoiding duplicates)
        setPNode nodes;
        for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
            mapElemIt elemIt = meshElems.find(*it);
            if (elemIt == meshElems.end())
                continue;

            for (uInt nodeIndex=0; nodeIndex<elemIt->second.nodes.size(); nodeIndex++)
                nodes.insert(elemIt->second.nodes[nodeIndex]);
        }

        // scale them
        for (setPNodeIt nodeIt=nodes.begin(); nodeIt!=nodes.end(); nodeIt++){
            for (uInt i=0; i<3; i++)
                (*nodeIt)->coordinates[i] = origin[i]-(origin[i]-(*nodeIt)->coordinates[i])*ratio[i];
        }
    }
    else{
        // scale all nodes
        for (mapNodeIt nodeIt=meshNodes.begin(); nodeIt!=meshNodes.end(); nodeIt++){
            Node* node = &nodeIt->second;
            for (uInt i=0; i<3; i++)
                node->coordinates[i] = origin[i]-(origin[i]-node->coordinates[i])*ratio[i];
        }
    }

    updateBounds();
}

//-----------------------------------------------------------------------/
// translate (selection)
//-----------------------------------------------------------------------/
void CMesh::translate(
        const double distance[3],
        const setUInt& nodeSelection,
        const setUInt& elemSelection)
{
    if (!nodeSelection.empty()){
        // translate node selection
        for (setConstUIntIt it=nodeSelection.begin(); it!=nodeSelection.end(); it++){
            mapNodeIt nodeIt = meshNodes.find(*it);
            if (nodeIt == meshNodes.end())
                continue;

            Node* node = &(nodeIt->second);
            for (uInt i=0; i<3; i++)
                node->coordinates[i] += distance[i];
        }
    }
    else if (!elemSelection.empty()){
        // translate element selection
        // extract nodes of all elements into one set (thus avoiding duplicates)
        setPNode nodes;
        for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
            mapElemIt elemIt = meshElems.find(*it);
            if (elemIt == meshElems.end())
                continue;

            for (uInt nodeIndex=0; nodeIndex<elemIt->second.nodes.size(); nodeIndex++)
                nodes.insert(elemIt->second.nodes[nodeIndex]);
        }

        // translate them
        for (setPNodeIt nodeIt=nodes.begin(); nodeIt!=nodes.end(); nodeIt++){
            for (uInt i=0; i<3; i++)
                (*nodeIt)->coordinates[i] += distance[i];
        }
    }
    else{
        // translate all nodes
        for (mapNodeIt nodeIt=meshNodes.begin(); nodeIt!=meshNodes.end(); nodeIt++){
            Node* node = &nodeIt->second;
            for (uInt i=0; i<3; i++)
                node->coordinates[i] += distance[i];
        }
    }

    updateBounds();
}

//-----------------------------------------------------------------------/
// rotate (selection)
//-----------------------------------------------------------------------/
void CMesh::rotate(
        const double angle,
        const double origin[3],
        const double tip[3],
        const setUInt& nodeSelection,
        const setUInt& elemSelection)
{
    if (!nodeSelection.empty()){
        // rotate node selection
        for (setConstUIntIt it=nodeSelection.begin(); it!=nodeSelection.end(); it++){
            mapNodeIt nodeIt = meshNodes.find(*it);
            if (nodeIt == meshNodes.end())
                continue;

            rotateNode(&(nodeIt->second), angle, origin, tip);
        }
    }
    else if (!elemSelection.empty()){
        // rotate element selection
        // extract nodes of all elements into one set (thus avoiding duplicates)
        setPNode nodes;
        for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
            mapElemIt elemIt = meshElems.find(*it);
            if (elemIt == meshElems.end())
                continue;

            for (uInt nodeIndex=0; nodeIndex<elemIt->second.nodes.size(); nodeIndex++)
                nodes.insert(elemIt->second.nodes[nodeIndex]);
        }

        // rotate them
        for (setPNodeIt nodeIt=nodes.begin(); nodeIt!=nodes.end(); nodeIt++)
            rotateNode(*nodeIt, angle, origin, tip);
    }
    else{
        // rotate all nodes
        for (mapNodeIt nodeIt=meshNodes.begin(); nodeIt!=meshNodes.end(); nodeIt++)
            rotateNode(&nodeIt->second, angle, origin, tip);
    }

    updateBounds();
}

//-----------------------------------------------------------------------/
// rotate node
// function found on web attributed to Ronald Goldman
// http://paulbourke.net/geometry/rotate/
//-----------------------------------------------------------------------/
void CMesh::rotateNode(
        Node* node,
        const double angle,
        const double origin[3],
        const double tip[3])
{
    // calculate rotation axis direction (unit vector)
    double axis[3];
    for (uInt i=0; i<3; i++)
        axis[i] = tip[i] - origin[i];
    // normalize
    double axisMod = sqrt(axis[0]*axis[0] +
                          axis[1]*axis[1] +
                          axis[2]*axis[2] );
    if (axisMod == 0.0)
        return;
    axis[0] /= axisMod;
    axis[1] /= axisMod;
    axis[2] /= axisMod;

    // extract node coordinates (with origin subtracted)
    double coordinates[3];
    for (uInt i=0; i<3; i++)
        coordinates[i] = node->coordinates[i] - origin[i];

    double cosAng = cos(angle);
    double sinAng = sin(angle);
    double minCosAng = 1-cosAng;

    vecDbl modCoords(3, 0.0);
    modCoords[0] += (minCosAng*axis[0]*axis[0] +         cosAng) * coordinates[0];
    modCoords[0] += (minCosAng*axis[0]*axis[1] - axis[2]*sinAng) * coordinates[1];
    modCoords[0] += (minCosAng*axis[0]*axis[2] + axis[1]*sinAng) * coordinates[2];

    modCoords[1] += (minCosAng*axis[0]*axis[1] + axis[2]*sinAng) * coordinates[0];
    modCoords[1] += (minCosAng*axis[1]*axis[1] +         cosAng) * coordinates[1];
    modCoords[1] += (minCosAng*axis[1]*axis[2] - axis[0]*sinAng) * coordinates[2];

    modCoords[2] += (minCosAng*axis[0]*axis[2] - axis[1]*sinAng) * coordinates[0];
    modCoords[2] += (minCosAng*axis[1]*axis[2] + axis[0]*sinAng) * coordinates[1];
    modCoords[2] += (minCosAng*axis[2]*axis[2] +         cosAng) * coordinates[2];

    modCoords[0] += origin[0];
    modCoords[1] += origin[1];
    modCoords[2] += origin[2];

    for (uInt i=0; i<3; i++)
        node->coordinates[i] = modCoords[i];
}

//-----------------------------------------------------------------------/
// mirror (element selection)
//-----------------------------------------------------------------------/
void CMesh::mirror(
        const double origin[3],
        const double tip[3],
        const setUInt& elemSelection)
{
    // get normal unit vector
    double normal[3];
    for (uInt i=0; i<3; i++)
        normal[i] = tip[i]-origin[i];
    // get modulus
    double normalMod = sqrt(normal[0]*normal[0] +
                            normal[1]*normal[1] +
                            normal[2]*normal[2] );
    if (normalMod == 0.0)
        return;
    // normalize
    normal[0] /= normalMod;
    normal[1] /= normalMod;
    normal[2] /= normalMod;

    // store ID offsets as maximum ID numbers
    uInt nodeIDOffset = maxNodeId;
    uInt elemIDOffset = maxElemId;

    // assemble set of node and element pointers
    setPNode nodes;
    setPElem elems;
    if (elemSelection.empty()){
        // add all nodes
        for (mapNodeIt nodeIt=meshNodes.begin(); nodeIt!=meshNodes.end(); nodeIt++)
            nodes.insert(&nodeIt->second);
        // add all elements
        for (mapElemIt elemIt=meshElems.begin(); elemIt!=meshElems.end(); elemIt++)
            elems.insert(&elemIt->second);
    }
    else{
        // add nodes of elements in element selection
        for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
            // try to find this element in map
            mapElemIt elemIt = meshElems.find(*it);
            if (elemIt == meshElems.end())
                continue;
            // add element pointer
            elems.insert(&elemIt->second);
            // add node pointers
            for (uInt nodeIndex=0; nodeIndex<elemIt->second.nodes.size(); nodeIndex++)
                nodes.insert(elemIt->second.nodes[nodeIndex]);
        }
    }

    // mirror all nodes first
    for (setPNodeIt nodeIt=nodes.begin(); nodeIt!=nodes.end(); nodeIt++){
        // get node ID
        uInt nodeID = (*nodeIt)->id+nodeIDOffset;
        // calculate mirrored node coordinates
        double v[3], coordinates[3];
        for (uInt i=0; i<3; i++)
            v[i] = (*nodeIt)->coordinates[i]-origin[i];
        double dblDist = (v[0]*normal[0] + v[1]*normal[1] + v[2]*normal[2])*2;
        for (uInt i=0; i<3; i++)
            coordinates[i] = (*nodeIt)->coordinates[i]-dblDist*normal[i];
        // insert node
        forceInsertNodeWithIDCoords(nodeID, coordinates, numNode-1);
    }

    // mirror elements
    for (setPElemIt elemIt=elems.begin(); elemIt!=elems.end(); elemIt++){
        // get element ID
        uInt elemID = (*elemIt)->id+elemIDOffset;
        // get element type
        uInt type = (*elemIt)->type;
        // get element nodes
        vecUInt nodeIDs((*elemIt)->nodes.size());
        nodeIDs[0] = (*elemIt)->nodes[2]->id+nodeIDOffset; // switched first and third (to keep normal orientation)
        nodeIDs[1] = (*elemIt)->nodes[1]->id+nodeIDOffset;
        nodeIDs[2] = (*elemIt)->nodes[0]->id+nodeIDOffset; // switched first and third (to keep normal orientation)
        if (nodeIDs.size() == 4) nodeIDs[3] = (*elemIt)->nodes[3]->id+nodeIDOffset;
        // insert element
        insertElemWithIDTypeNodeIDs(elemID, type, nodeIDs);
    }

    // mirror neighboring information
    for (setPElemIt elemIt=elems.begin(); elemIt!=elems.end(); elemIt++){
        // get element ID
        uInt elemID = (*elemIt)->id+elemIDOffset;
        // get neighbor IDs
        vecUInt neighbors((*elemIt)->neighbors.size());
        uInt i = 0;
        for (setPElemIt neighIt=(*elemIt)->neighbors.begin(); neighIt!=(*elemIt)->neighbors.end(); neighIt++)
            neighbors[i++] = (*neighIt)->id+elemIDOffset;
        insertNeighborsToElem(elemID, neighbors);
    }

    // mirror wake information
    for (setPElemIt elemIt=elems.begin(); elemIt!=elems.end(); elemIt++){
        // get element ID
        uInt elemID = (*elemIt)->id+elemIDOffset;
        // get trail IDs
        vecUInt trailIDs;
        bool trailExists = false;
        if ((*elemIt)->trailUpper){
            trailExists = true;
            trailIDs.resize(1);
            trailIDs[0] = (*elemIt)->trailUpper->id+elemIDOffset;
        }
        if ((*elemIt)->trailLower){
            trailExists = true;
            trailIDs.resize(2);
            trailIDs[1] = (*elemIt)->trailLower->id+elemIDOffset;
        }
        if (trailExists)
            insertTrailsToElem(elemID, trailIDs);
    }

    updateBounds();
}

//-----------------------------------------------------------------------/
// convert element type
//   0 = OK
//  -1 = given element ID does not exist
//  -2 = incorrect requested type ("body", "dummy", "wake")
//-----------------------------------------------------------------------/
int CMesh::convertElemType(const uInt& id, const string& convertTo)
{
    // extract element with given ID (check if it exists)
    mapElemIt it = meshElems.find(id);
    if (it == meshElems.end())
        return -1;
    Elem* element = &it->second;
    uInt oldType = element->type;

    // check to which element type to convert
    uInt newType;
    if (element->nodes.size() == 3){
        if      (convertTo == "body" ) newType = 2;
        else if (convertTo == "dummy") newType = 21;
        else if (convertTo == "wake" ) newType = 11;
        else return -2;
    }
    else{
        if      (convertTo == "body" ) newType = 1;
        else if (convertTo == "dummy") newType = 20;
        else if (convertTo == "wake" ) newType = 10;
        else return -2;
    }

    // exit if types are the same
    if (oldType == newType)
        return 0;

    // change type
    element->type = newType;

    // if conversion from wake element
    if (oldType==10 || oldType==11){
        // delete references to current element in upper and lower trailing elements
        if (element->trailUpper) element->trailUpper->wakes.erase(element);
        if (element->trailLower) element->trailLower->wakes.erase(element);
        // delete references to upper and lower trailing elements in current element
        element->trailUpper = NULL;
        element->trailLower = NULL;
    }

    // if conversion to wake element
    if (newType==10 || newType==11){
        // delete references to current element in neighbors
        for (setPElemIt it=element->neighbors.begin(); it!=element->neighbors.end(); it++)
            (*it)->neighbors.erase(element);
        // delete neighbors
        element->neighbors.clear();
        // delete references to current element in wakes (references can be in upper and lower "trails")
        for (setPElemIt it=element->wakes.begin(); it!=element->wakes.end(); it++){
            if ((*it)->trailUpper == element) (*it)->trailUpper = NULL;
            if ((*it)->trailLower == element) (*it)->trailLower = NULL;
        }
        // delete wakes
        element->wakes.clear();
    }

    // if conversion to dummy element
    if (newType==20 || newType==21){
        // delete references to current element in neighbors
        for (setPElemIt it=element->neighbors.begin(); it!=element->neighbors.end(); it++)
            (*it)->neighbors.erase(element);
        // delete neighbors
        element->neighbors.clear();
    }

    // recount elements
    if      (oldType == 1 ) numBodyQuad--;
    else if (oldType == 2 ) numBodyTria--;
    else if (oldType == 10) numWakeQuad--;
    else if (oldType == 11) numWakeTria--;
    else if (oldType == 20) numDummQuad--;
    else if (oldType == 21) numDummTria--;
    if      (newType == 1 ) numBodyQuad++;
    else if (newType == 2 ) numBodyTria++;
    else if (newType == 10) numWakeQuad++;
    else if (newType == 11) numWakeTria++;
    else if (newType == 20) numDummQuad++;
    else if (newType == 21) numDummTria++;

    // update bounds
    updateElemBounds();

    return 0;
}

//-----------------------------------------------------------------------/
// updates wake trailing elements
//-----------------------------------------------------------------------/
int CMesh::updateWakeTrails(
        const uInt& upperID,
        const uInt& lowerID,
        const setUInt& elemSelection)
{
    // check if both elements exist and extract element pointers
    mapElemIt upper = meshElems.find(upperID);
    mapElemIt lower = meshElems.find(lowerID);
    if (upper==meshElems.end() || lower==meshElems.end())
        return -1;
    Elem* trailUpper = &(upper->second);
    Elem* trailLower = &(lower->second);

    // go through element selection and update trails
    bool nonWakeFound = false;
    for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
        // find element in map
        mapElemIt elemIt = meshElems.find(*it);
        if (elemIt == meshElems.end())
            continue;

        // check if wake type and insert trails
        if (elemIt->second.type==10 || elemIt->second.type==11){
            elemIt->second.trailUpper = trailUpper;
            elemIt->second.trailLower = trailLower;
        }
        else
            nonWakeFound = true;
    }

    if (nonWakeFound) return 1;
    else              return 0;
}

//-----------------------------------------------------------------------/
// extrudes line from node selection
//-----------------------------------------------------------------------/
void CMesh::extrudeLine(
        const double& extrusionLength,
        const uInt& elemNum,
        const setUInt& nodeSelection)
{
    // get node line branches
    vecDeqPNode nodeBranches;
    getNodeBranches(nodeSelection, nodeBranches);

    // get total and single wake element length
    double totalExtrusionLength = extrusionLength*(maxBody[1]-minBody[1]);
    double singleExtrusionLength = totalExtrusionLength/elemNum;

    // for each branch
    for (uInt i=0; i<nodeBranches.size(); i++){
        // for each pair of nodes from line branch
        for (deqPNodeIt nodeIt=nodeBranches[i].begin(); nodeIt!=nodeBranches[i].end()-1; nodeIt++){
            // initiate starting coordinates
            vecVecDbl coordinates(4);
            for (uInt i=0; i<4; i++)
                coordinates[i].resize(3);
            for (uInt i=0; i<3; i++){
                coordinates[0][i] = (*nodeIt  )->coordinates[i];
                deqPNodeIt next = nodeIt;
                coordinates[1][i] = (*(++next))->coordinates[i];
            }

            // perform element creation number of segments times
            for (uInt segmentIndex=0; segmentIndex<elemNum; segmentIndex++){
                // assemble wake coordinates
                coordinates[2] = coordinates[1];
                coordinates[3] = coordinates[0];
                // add segment length to X coordinates
                coordinates[2][0] += singleExtrusionLength;
                coordinates[3][0] += singleExtrusionLength;

                // insert wake element
                Elem* newElemP = NULL;
                insertElemWithNodeCoordsPrivate(coordinates, 1, newElemP);

                // store end coordinates as start coordinates for next segment
                coordinates[0] = coordinates[3];
                coordinates[1] = coordinates[2];
            }
        }
    }
}

//-----------------------------------------------------------------------/
// separates neighbors
//-----------------------------------------------------------------------/
void CMesh::separateAtSeam(const setUInt& nodeSelection)
{
    // get node line branches
    vecDeqPNode nodeBranches;
    getNodeBranches(nodeSelection, nodeBranches);

    // for each branch
    for (uInt i=0; i<nodeBranches.size(); i++){
        // for each pair of nodes from line branch
        for (deqPNodeIt nodeIt=nodeBranches[i].begin(); nodeIt!=nodeBranches[i].end()-1; nodeIt++){
            // get nodes
            Node* firstNode = *(nodeIt  );
            Node* secndNode = *(nodeIt+1);

            // find two elements that have these two nodes
            Elem* firstElem = NULL;
            Elem* secndElem = NULL;
            for (setPElemIt elemItA=firstNode->elements.begin(); elemItA!=firstNode->elements.end(); elemItA++){
                for (setPElemIt elemItB=secndNode->elements.begin(); elemItB!=secndNode->elements.end(); elemItB++){
                    if ((*elemItA)->id != (*elemItB)->id)
                        continue;
                    if (!firstElem)
                        firstElem = *elemItA;
                    else if (!secndElem){
                        secndElem = *elemItA;
                        break;
                    }
                }
                if (secndElem)
                    break;
            }

            // remove neighboring references between these two elements if second element is found
            if (secndElem){
                firstElem->neighbors.erase(secndElem);
                secndElem->neighbors.erase(firstElem);
            }
        }
    }
}

//-----------------------------------------------------------------------/
// updates neighboring information
//-----------------------------------------------------------------------/
void CMesh::updateNeighboringInfo(const double& angleThreshold, const setUInt& elemSelection)
{
    // assemble set of element pointers (only body elements)
    setPElem workingSelection;
    if (elemSelection.empty()){
        // with all elements
        for (mapElemIt it=meshElems.begin(); it!=meshElems.end(); it++){
            if (it->second.type==1 || it->second.type==2)
                workingSelection.insert(&it->second);
        }
    }
    else{
        // with element selection
        for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
            // find element in map
            mapElemIt elemIt = meshElems.find(*it);
            if (elemIt == meshElems.end())
                continue;

            // check if body type and insert into working selection
            if (elemIt->second.type==1 || elemIt->second.type==2)
                workingSelection.insert(&(elemIt->second));
        }
    }

    // clear neighboring info
    for (setPElemIt it=workingSelection.begin(); it!=workingSelection.end(); it++){
        // clear references to this element from his neighbors
        for (setPElemIt neighIt=(*it)->neighbors.begin(); neighIt!=(*it)->neighbors.end(); neighIt++){
            for (setPElemIt neighRevIt=(*neighIt)->neighbors.begin(); neighRevIt!=(*neighIt)->neighbors.end(); neighRevIt++){
                if ((*neighRevIt)->id == (*neighIt)->id)
                    (*neighIt)->neighbors.erase((*neighIt));
            }
        }
        // clear references to neighbors from this element
        (*it)->neighbors.clear();
    }

    // find neighbors for each element
    for (setPElemIt it=workingSelection.begin(); it!=workingSelection.end(); it++){
        // extract element pointer
        Elem* element = *it;

        // try each node pair
        for (uInt firstNodeIndex=0; firstNodeIndex<element->nodes.size(); firstNodeIndex++){
            // get second node index (perturbation)
            uInt secondNodeIndex;
            if   (firstNodeIndex == element->nodes.size()-1) secondNodeIndex = 0;
            else                                             secondNodeIndex = firstNodeIndex+1;

            // find neighbor that has these two nodes
            Elem* neighbor = getNeighbor(element, firstNodeIndex, secondNodeIndex);

            // check if element is in working selection
            if (workingSelection.find(neighbor) == workingSelection.end())
                continue;

            // check if angle between element normals is below threshold
            if (angleBetweenElems(neighbor, element) > angleThreshold)
                continue;

            // insert neighbor to observed element and observed element to neighbor
            neighbor->neighbors.insert(element );
            element ->neighbors.insert(neighbor);
        }
    }
}

//-----------------------------------------------------------------------/
// updates wake information
//-----------------------------------------------------------------------/
void CMesh::updateWakeInfo(
        const double& angleThreshold,
        const double& wakeLengthRatio,
        const uInt& wakeNum,
        const setUInt& elemSelection)
{
    // assemble set of element pointers (only body elements)
    setPElem workingSelection;
    if (elemSelection.empty()){
        // with all elements that
        for (mapElemIt it=meshElems.begin(); it!=meshElems.end(); it++){
            // check that element is not wake
            if (it->second.type!=10 && it->second.type!=11)
                workingSelection.insert(&it->second);
        }
    }
    else{
        // with element selection
        for (setConstUIntIt it=elemSelection.begin(); it!=elemSelection.end(); it++){
            // extract element pointer from map
            Elem* pElem = &meshElems[*it];
            // check that element is not wake
            if (pElem->type!=10 && pElem->type!=11)
                workingSelection.insert(&meshElems[*it]);
        }
    }

    // delete wakes of each element
    vecUInt elemIDs;
    vecUInt nodeIDs;
    for (setPElemIt it=workingSelection.begin(); it!=workingSelection.end(); it++){
        // for each wake element
        for (setPElemIt wakeIt=(*it)->wakes.begin(); wakeIt!=(*it)->wakes.end(); wakeIt++){
            // store element ID
            elemIDs.push_back((*wakeIt)->id);

            // store node IDs of element
            for (uInt nodeIndex=0; nodeIndex<nodeIDs.size(); nodeIndex++)
                nodeIDs.push_back((*wakeIt)->nodes[nodeIndex]->id);
        }

        // delete wake references (just in case - this should have been cleaned with clearing references from trailing elements)
        (*it)->wakes.clear();
    }
    // delete elements
    deleteElemsByID(elemIDs);
    // delete nodes
    deleteNodesByID(nodeIDs);

    // calculate wake length
    double wakeLength = (maxBody[1]-minBody[1])*wakeLengthRatio;

    // get segment length
    double segmentLength = wakeLength/wakeNum;

    // find neighbors for each element
    setUInt usedElems;
    for (setPElemIt it=workingSelection.begin(); it!=workingSelection.end(); it++){
        // extract element pointer
        Elem* element = *it;

        // check if element already used in wake creation
        if (usedElems.find(element->id) != usedElems.end())
            continue;

        // check if element normal is pointing forward
        if (calcNormal(element)[0] < 0.0)
            continue;

        // with each node pair of current element, try to find neighbor
        for (uInt firstNodeIndex=0; firstNodeIndex<element->nodes.size(); firstNodeIndex++){
            // get second node index (perturbation)
            uInt secondNodeIndex;
            if   (firstNodeIndex == element->nodes.size()-1) secondNodeIndex = 0;
            else                                             secondNodeIndex = firstNodeIndex+1;

            // find neighbor that has these two nodes
            Elem* neighbor = getNeighbor(element, firstNodeIndex, secondNodeIndex);

            // check if element is in working selection
            if (workingSelection.find(neighbor) == workingSelection.end())
                continue;

            // check if neighbor already used in wake creation
            if (usedElems.find(neighbor->id) != usedElems.end())
                continue;

            // check if angle between element normals is above threshold
            if (angleBetweenElems(neighbor, element) < angleThreshold)
                continue;

            // mark these two elements as already used in wake creation
            usedElems.insert(element ->id);
            usedElems.insert(neighbor->id);

            // get starting node coordinates (normal should point upwards for convenience)
            vecVecDbl startCoords(2);
            for (uInt i=0; i<startCoords.size(); i++)
                startCoords[i].resize(3);
            if (element->nodes[ firstNodeIndex]->coordinates[1] > element->nodes[secondNodeIndex]->coordinates[1] ){
                for (uInt coordIndex=0; coordIndex<3; coordIndex++){
                    startCoords[0][coordIndex] = element->nodes[ firstNodeIndex]->coordinates[coordIndex];
                    startCoords[1][coordIndex] = element->nodes[secondNodeIndex]->coordinates[coordIndex];
                }
            }
            else{
                for (uInt coordIndex=0; coordIndex<3; coordIndex++){
                    startCoords[1][coordIndex] = element->nodes[ firstNodeIndex]->coordinates[coordIndex];
                    startCoords[0][coordIndex] = element->nodes[secondNodeIndex]->coordinates[coordIndex];
                }
            }

            // perform wake element creation number of segments times
            vecVecDbl coordinates(4);
            coordinates[0] = startCoords[0];
            coordinates[1] = startCoords[1];
            for (uInt segmentIndex=0; segmentIndex<wakeNum; segmentIndex++){
                // assemble wake coordinates
                coordinates[2] = coordinates[1];
                coordinates[3] = coordinates[0];
                // add segment length to X coordinates
                coordinates[2][0] += segmentLength;
                coordinates[3][0] += segmentLength;

                // insert wake element
                Elem* newElemP = NULL;
                insertElemWithNodeCoordsPrivate(coordinates, 10, newElemP);

                // add this wake to both trailing elements
                (element )->wakes.insert(newElemP);
                (neighbor)->wakes.insert(newElemP);

                // add trailing elements to this wake element, but watch for orientation
                if (calcNormal(element)[2] > 0.0){
                    newElemP->trailUpper = (element );
                    newElemP->trailLower = (neighbor);
                }
                else{
                    newElemP->trailUpper = (neighbor);
                    newElemP->trailLower = (element );
                }

                // store end coordinates as start coordinates for next segment
                coordinates[0] = coordinates[3];
                coordinates[1] = coordinates[2];
            }

            break;
        }
    }
}

//-----------------------------------------------------------------------/
// clear neighboring info
//-----------------------------------------------------------------------/
void CMesh::clearNeighboringInfo()
{
    for (mapElemIt elemMapIt=meshElems.begin(); elemMapIt!=meshElems.end(); elemMapIt++)
        elemMapIt->second.neighbors.clear();
}

//-----------------------------------------------------------------------/
// clear mesh data
//-----------------------------------------------------------------------/
void CMesh::clearData()
{
    meshNodes.clear();
    meshElems.clear();

    initializeData();
}

//-----------------------------------------------------------------------/
// delete all nodes
//-----------------------------------------------------------------------/
void CMesh::deleteAllNodes()
{
    // delete node from map (check is performed that node is not referenced by elements)
    for (mapNodeIt it=meshNodes.begin(); it!=meshNodes.end();){
        // store iterator before deleting it
        mapNodeIt savedIt = it++;
        deleteNode(savedIt);
    }

    // update bounds
    updateNodeBounds();

    // update maximum node ID
    updateMaxNodeID();
}

//-----------------------------------------------------------------------/
// delete nodes by IDs
//-----------------------------------------------------------------------/
void CMesh::deleteNodesByID(const vecUInt& nodeIDs)
{
    // go through all ID's
    for (uInt nodeIndex=0; nodeIndex<nodeIDs.size(); nodeIndex++)
        deleteNode(meshNodes.find(nodeIDs[nodeIndex]));

    // update bounds
    updateNodeBounds();

    // update maximum node ID
    updateMaxNodeID();
}

//-----------------------------------------------------------------------/
// delete all elements
//-----------------------------------------------------------------------/
void CMesh::deleteAllElems()
{
    // clear element map
    meshElems.clear();

    // clear elements from nodes
    for (mapNodeIt nodeMapIt=meshNodes.begin(); nodeMapIt!=meshNodes.end(); nodeMapIt++)
        nodeMapIt->second.elements.clear();

    // initialize element data
    initializeElemData();
}

//-----------------------------------------------------------------------/
// delete elements by IDs
//-----------------------------------------------------------------------/
void CMesh::deleteElemsByID(const vecUInt& elemIDs)
{
    // go through all ID's
    for (uInt elemIndex=0; elemIndex<elemIDs.size(); elemIndex++){
        // extract element with given ID
        mapElemIt it = meshElems.find(elemIDs[elemIndex]);
        if (it != meshElems.end()){
            Elem* element = &it->second;

            // removing element references from nodes
            for (uInt i=0; i<element->nodes.size(); i++)
                element->nodes[i]->elements.erase(element);

            // removing neighboring references from elements
            for (setPElemIt elemIt=element->neighbors.begin(); elemIt!=element->neighbors.end(); elemIt++)
                (*elemIt)->neighbors.erase(element);

            // remove wake information
            for (setPElemIt elemIt=element->wakes.begin(); elemIt!=element->wakes.end(); elemIt++){
                if ((*elemIt)->trailUpper == element) (*elemIt)->trailUpper = NULL;
                if ((*elemIt)->trailLower == element) (*elemIt)->trailLower = NULL;
            }

            // delete references to current element in trailing elements
            if (element->trailUpper) element->trailUpper->wakes.erase(element);
            if (element->trailLower) element->trailLower->wakes.erase(element);

            // erase element from map
            meshElems.erase(it);

            // update counter
            if      (element->type == 1 ) numBodyQuad--;
            else if (element->type == 2 ) numBodyTria--;
            else if (element->type == 10) numWakeQuad--;
            else if (element->type == 11) numWakeTria--;
            else if (element->type == 20) numDummQuad--;
            else if (element->type == 21) numDummTria--;
        }
    }

    // update bounds
    updateElemBounds();

    // update maximum element ID
    updateMaxElemID();
}

///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
//                    P R I V A T E   M E M B E R S                          //
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////

//-----------------------------------------------------------------------/
// initialize mesh data
//-----------------------------------------------------------------------/
void CMesh::initializeData()
{
    // initialize node data
    initializeNodeData();

    // initialize element data
    initializeElemData();
}

//-----------------------------------------------------------------------/
// initialize node data
//-----------------------------------------------------------------------/
void CMesh::initializeNodeData()
{
    // initialize general node data
    maxNodeId = 0;
    numNode = 0;
    // initialize node bounds
    for (uInt i=0; i<3; i++) minNode[i] = MAX_BOUND;
    for (uInt i=0; i<3; i++) maxNode[i] = MIN_BOUND;
}

//-----------------------------------------------------------------------/
// initialize element data
//-----------------------------------------------------------------------/
void CMesh::initializeElemData()
{
    // initialize general element data
    maxElemId = 0;
    numBodyQuad = 0;
    numBodyTria = 0;
    numDummQuad = 0;
    numDummTria = 0;
    numWakeQuad = 0;
    numWakeTria = 0;
    // initialize body bounds
    for (uInt i=0; i<3; i++) minBody[i] = MAX_BOUND;
    for (uInt i=0; i<3; i++) maxBody[i] = MIN_BOUND;
    // initialize wake bounds
    for (uInt i=0; i<3; i++) minWake[i] = MAX_BOUND;
    for (uInt i=0; i<3; i++) maxWake[i] = MIN_BOUND;
}

//-----------------------------------------------------------------------/
// get neighbor of the element (pointers)
// returned value:
//  element pointer
//-----------------------------------------------------------------------/
Elem* CMesh::getNeighbor(
        const Elem* element,
        const uInt& firstNodeIndex,
        const uInt& secondNodeIndex)
{
    Elem* neighbor = NULL;

    // check element validity
    if (!element)
        return neighbor;

    // check indexes
    if (firstNodeIndex  >= element->nodes.size() ||
        secondNodeIndex >= element->nodes.size() )
        return neighbor;

    // loop through elements of first node
    for (setPElemIt elemItI = element->nodes[firstNodeIndex]->elements.begin();
         elemItI != element->nodes[firstNodeIndex]->elements.end();
         elemItI++){
        // loop through elements of second node
        for (setPElemIt elemItJ = element->nodes[secondNodeIndex]->elements.begin();
             elemItJ != element->nodes[secondNodeIndex]->elements.end();
             elemItJ++){
            // check that these nodes do not point to different elements
            if (*elemItI != *elemItJ)
                continue;

            // check that potential neighbor is not in fact the observed element
            if (*elemItI == element)
                continue;

            // this is the neighbor so return it
            return *elemItI;
        }
    }

    return neighbor;
}

//-----------------------------------------------------------------------/
// insert node with defined coordinates
// returned value:
//  node pointer
//-----------------------------------------------------------------------/
Node* CMesh::insertNodeWithCoordsPrivate(const double coordinates[3])
{
    // get new node ID
    uInt id = ++maxNodeId;

    // add new node to node map
    Node* newNodeP = &meshNodes[id] ;

    // add node ID
    newNodeP->id = id;
    // add coordinates
    for (uInt i=0; i<3; i++)
        newNodeP->coordinates[i] = coordinates[i];

    // update bounds
    updateBoundsWithNode(newNodeP);

    // increase node counter
    numNode++;

    return newNodeP;
}

//-----------------------------------------------------------------------/
// insert element with type and coordinates of constituing nodes
//
// returned codes:
//  >0 = element ID
//  -1 = number of coordinates per node is not 3
//  -2 = zero surface element
//  -3 = incorrect element type
//  -4 = number of nodes is not 3 or 4
//-----------------------------------------------------------------------/
int CMesh::insertElemWithNodeCoordsPrivate(const vecVecDbl& coordinatesAll, const uInt& type, Elem*& newElemP)
{
    // check if incorrect type is given
    if (type != 1  &&
        type != 2  &&
        type != 10 &&
        type != 11 &&
        type != 20 &&
        type != 21 )
        return -3;

    // check number of nodes
    if (coordinatesAll.size()!=4 && coordinatesAll.size()!=3)
        return -4;

    // check number of coordinates for each node
    for (uInt i=0; i<coordinatesAll.size(); i++){
        if (coordinatesAll[i].size() != 3)
            return -1;
    }

    // check normal
    if (calcNormal(coordinatesAll) != 0)
        return -2;

    // get new element ID
    uInt id = ++maxElemId;

    // add new element with given ID to element map
    newElemP = &meshElems[id];

    // add element ID
    newElemP->id = id;
    // add element type
    newElemP->type = type;
    // add constituing nodes (pointers)
    newElemP->nodes.resize(coordinatesAll.size());
    for (uInt i=0; i<coordinatesAll.size(); i++){
        // create node and add to element constituing node pointers
        newElemP->nodes[i] = insertNodeWithCoordsPrivate(&coordinatesAll[i][0]);
        // add element associativity to node
        newElemP->nodes[i]->elements.insert(newElemP);
    }

    // update bounds with element
    updateBoundsWithElem(newElemP);

    // add to counter of elements
    if      (type == 1 ) numBodyQuad++;
    else if (type == 2 ) numBodyTria++;
    else if (type == 10) numWakeQuad++;
    else if (type == 11) numWakeTria++;
    else if (type == 20) numDummQuad++;
    else if (type == 21) numDummTria++;

    return id;
}

//-----------------------------------------------------------------------/
// insert node defined with ID and coordinates without checking input data
//  id          - node ID
//  coordinates - node coordinates
//  position    - (optional for speed-up) position of node in node map
//-----------------------------------------------------------------------/
void CMesh::forceInsertNodeWithIDCoords(const uInt& id, const double coordinates[3], const int& position)
{
    // create new node pointer
    Node* newNodeP;
    // check if position is given
    if (position == -1)
        // no position given so just insert node to node map
        newNodeP = &meshNodes[id];
    else{
        // position is given so define new node
        Node newNode;
        // and add it to node map with position hint
        mapNodeIt p  = meshNodes.lower_bound(position);
        mapNodeIt it = meshNodes.insert(p, pair<uInt,Node>(id, newNode));
        newNodeP = &it->second;
    }

    // add node ID
    newNodeP->id = id;
    // add coordinates
    for (uInt i=0; i<3; i++)
        newNodeP->coordinates[i] = coordinates[i];

    // check if node ID is greatest yet
    if (id > maxNodeId)
        maxNodeId = id;

    // update bounds
    updateBoundsWithNode(newNodeP);

    // increase node counter
    numNode++;
}

//-----------------------------------------------------------------------/
// insert element with ID, type and constituing nodes ID's without
// checking input data
//  id          - element ID
//  type        - element type
//  nodes       - constituing node pointers
//  position    - (optional for speed-up) position of element in element map
//-----------------------------------------------------------------------/
void CMesh::forceInsertElemWithIDNodes(const uInt& id, const uInt& type, const vecPNode& nodes, const int& position)
{
    // create new element pointer
    Elem* newElemP;
    // check if position is given
    if (position == -1)
        // no position given so just insert element to element map
        newElemP = &meshElems[id];
    else{
        // position is given so define new element
        Elem newElem;
        // and add it to element map with position hint
        map<uInt,Elem>::iterator p  = meshElems.lower_bound(position);
        map<uInt,Elem>::iterator it = meshElems.insert(p, pair<uInt,Elem>(id, newElem));
        newElemP = &it->second;
    }

    // add element ID
    newElemP->id = id;
    // add element type
    newElemP->type = type;
    // add constituing nodes (pointers)
    newElemP->nodes.resize(nodes.size());
    for (uInt i=0; i<nodes.size(); i++)
        newElemP->nodes[i] = &meshNodes[nodes[i]->id];
    // add element associativity to nodes
    for (uInt i=0; i<newElemP->nodes.size(); i++)
        newElemP->nodes[i]->elements.insert(newElemP);

    // check if element ID is greatest yet
    if (id > maxElemId)
        maxElemId = id;

    // update bounds with element
    updateBoundsWithElem(newElemP);

    // add to counter of elements
    if      (type == 1 ) numBodyQuad++;
    else if (type == 2 ) numBodyTria++;
    else if (type == 10) numWakeQuad++;
    else if (type == 11) numWakeTria++;
    else if (type == 20) numDummQuad++;
    else if (type == 21) numDummTria++;
}

//-----------------------------------------------------------------------/
// delete node by given valid iterator to its position in node map
//-----------------------------------------------------------------------/
void CMesh::deleteNode(mapNodeIt nodeIt)
{
    if (nodeIt!=meshNodes.end() && nodeIt->second.elements.empty()){
        meshNodes.erase(nodeIt);
        numNode--;
    }
}

//-----------------------------------------------------------------------/
// delete element by given valid iterator to its position in element map
//-----------------------------------------------------------------------/
void CMesh::deleteElement(mapElemIt elemIt)
{
    if (elemIt!=meshElems.end()){
        if      (elemIt->second.type == 1 ) numBodyQuad--;
        else if (elemIt->second.type == 2 ) numBodyTria--;
        else if (elemIt->second.type == 10) numWakeQuad--;
        else if (elemIt->second.type == 11) numWakeTria--;
        else if (elemIt->second.type == 20) numDummQuad--;
        else if (elemIt->second.type == 21) numDummTria--;
        meshElems.erase(elemIt);
    }
}

//-----------------------------------------------------------------------/
// this function updates all bounds
//-----------------------------------------------------------------------/
void CMesh::updateBounds()
{
    updateNodeBounds();
    updateElemBounds();
}

//-----------------------------------------------------------------------/
// this function updates node bounds
//-----------------------------------------------------------------------/
void CMesh::updateNodeBounds()
{
    // initialize node bounds
    for (uInt i=0; i<3; i++) minNode[i] = MAX_BOUND;
    for (uInt i=0; i<3; i++) maxNode[i] = MIN_BOUND;

    for (mapNodeIt it=meshNodes.begin(); it!=meshNodes.end(); it++)
        updateBoundsWithNode(&it->second);
}

//-----------------------------------------------------------------------/
// this function updates element bounds
//-----------------------------------------------------------------------/
void CMesh::updateElemBounds()
{
    // initialize body bounds
    for (uInt i=0; i<3; i++) minBody[i] = MAX_BOUND;
    for (uInt i=0; i<3; i++) maxBody[i] = MIN_BOUND;
    // initialize wake bounds
    for (uInt i=0; i<3; i++) minWake[i] = MAX_BOUND;
    for (uInt i=0; i<3; i++) maxWake[i] = MIN_BOUND;

    for (mapElemIt it=meshElems.begin(); it!=meshElems.end(); it++)
        updateBoundsWithElem(&it->second);
}

//-----------------------------------------------------------------------/
// this function updates node bounds with given node
//-----------------------------------------------------------------------/
void CMesh::updateBoundsWithNode(const Node* node)
{
    for (uInt xyzIndex=0; xyzIndex<3; xyzIndex++){
        if (node->coordinates[xyzIndex] < minNode[xyzIndex]) minNode[xyzIndex] = node->coordinates[xyzIndex];
        if (node->coordinates[xyzIndex] > maxNode[xyzIndex]) maxNode[xyzIndex] = node->coordinates[xyzIndex];
    }
}

//-----------------------------------------------------------------------/
// this function updates body and node bounds with given element
//-----------------------------------------------------------------------/
void CMesh::updateBoundsWithElem(const Elem* element)
{
    // update node bounds
    for (uInt i=0; i<element->nodes.size(); i++){
        for (uInt xyzIndex=0; xyzIndex<3; xyzIndex++){
            if (element->nodes[i]->coordinates[xyzIndex] < minNode[xyzIndex]) minNode[xyzIndex] = element->nodes[i]->coordinates[xyzIndex];
            if (element->nodes[i]->coordinates[xyzIndex] > maxNode[xyzIndex]) maxNode[xyzIndex] = element->nodes[i]->coordinates[xyzIndex];
        }
    }

    // check element type
    if (element->type == 1  ||
        element->type == 2  ||
        element->type == 20 ||
        element->type == 21 ){
        // update body bounds
        for (uInt i=0; i<element->nodes.size(); i++){
            for (uInt xyzIndex=0; xyzIndex<3; xyzIndex++){
                if (element->nodes[i]->coordinates[xyzIndex] < minBody[xyzIndex]) minBody[xyzIndex] = element->nodes[i]->coordinates[xyzIndex];
                if (element->nodes[i]->coordinates[xyzIndex] > maxBody[xyzIndex]) maxBody[xyzIndex] = element->nodes[i]->coordinates[xyzIndex];
            }
        }
    }
    else if (element->type == 10 ||
             element->type == 11 ){
        // update wake bounds
        for (uInt i=0; i<element->nodes.size(); i++){
            for (uInt xyzIndex=0; xyzIndex<3; xyzIndex++){
                if (element->nodes[i]->coordinates[xyzIndex] < minWake[xyzIndex]) minWake[xyzIndex] = element->nodes[i]->coordinates[xyzIndex];
                if (element->nodes[i]->coordinates[xyzIndex] > maxWake[xyzIndex]) maxWake[xyzIndex] = element->nodes[i]->coordinates[xyzIndex];
            }
        }
    }
}

//-----------------------------------------------------------------------/
// this function updates maximum node ID
//-----------------------------------------------------------------------/
void CMesh::updateMaxNodeID()
{
    maxNodeId = 0;
    for (mapNodeIt it=meshNodes.begin(); it!=meshNodes.end(); it++){
        if (it->second.id > maxNodeId)
            maxNodeId = it->second.id;
    }
}

//-----------------------------------------------------------------------/
// this function updates maximum element ID
//-----------------------------------------------------------------------/
void CMesh::updateMaxElemID()
{
    maxElemId = 0;
    for (mapElemIt it=meshElems.begin(); it!=meshElems.end(); it++){
        if (it->second.id > maxElemId)
            maxElemId = it->second.id;
    }
}

//-----------------------------------------------------------------------/
// this function calculates normal or checks if zero surface element is
// produced with given coordinates
// error codes:
//   0 - OK
//  -1 - zero surface element
//-----------------------------------------------------------------------/
int CMesh::calcNormal(const vecVecDbl& coordinates, bool onlyCheck, double* normal)
{
    // calculating diagonal vectors
    double a[3], b[3];
                                       for (uInt i=0; i<3; i++) a[i] = coordinates[2][i]-coordinates[0][i];
    if      (coordinates.size() == 4){ for (uInt i=0; i<3; i++) b[i] = coordinates[3][i]-coordinates[1][i]; }
    else if (coordinates.size() == 3){ for (uInt i=0; i<3; i++) b[i] = coordinates[2][i]-coordinates[1][i]; }

    if (onlyCheck){
        // calculating vector normal to diagonals (vector product)
        double normalVec[3];
        normalVec[0] = a[1]*b[2] - a[2]*b[1];
        normalVec[1] = a[2]*b[0] - a[0]*b[2];
        normalVec[2] = a[0]*b[1] - a[1]*b[0];

        // calculating modulus of that vector
        double normalMod;
        if ((normalMod = sqrt(pow(normalVec[0], 2) +
                              pow(normalVec[1], 2) +
                              pow(normalVec[2], 2) )) == 0.0)
            return -1;
    }
    else{
        // calculating vector normal to diagonals (vector product)
        normal[0] = a[1]*b[2] - a[2]*b[1];
        normal[1] = a[2]*b[0] - a[0]*b[2];
        normal[2] = a[0]*b[1] - a[1]*b[0];

        // calculating modulus of that vector
        double normalMod;
        if ((normalMod = sqrt(pow(normal[0], 2) +
                              pow(normal[1], 2) +
                              pow(normal[2], 2) )) == 0.0)
            return -1;

        // normalize normal
        for (uInt i=0; i<3; i++)
            normal[i] /= normalMod;
    }

    return 0;
}

//-----------------------------------------------------------------------/
// calculates normal or checks if zero surface element is
// produced with given node pointers
// error codes:
//   0 - OK
//  -1 - zero surface element
//-----------------------------------------------------------------------/
int CMesh::calcNormal(const vecPNode& nodes, bool onlyCheck, double* normal)
{
    // calculating diagonal vectors
    double a[3], b[3];
                                 for (uInt i=0; i<3; i++) a[i] = nodes[2]->coordinates[i]-nodes[0]->coordinates[i];
    if      (nodes.size() == 4){ for (uInt i=0; i<3; i++) b[i] = nodes[3]->coordinates[i]-nodes[1]->coordinates[i]; }
    else if (nodes.size() == 3){ for (uInt i=0; i<3; i++) b[i] = nodes[2]->coordinates[i]-nodes[1]->coordinates[i]; }

    if (onlyCheck){
        // calculating vector normal to diagonals (vector product)
        double normalVector[3];
        normalVector[0] = a[1]*b[2] - a[2]*b[1];
        normalVector[1] = a[2]*b[0] - a[0]*b[2];
        normalVector[2] = a[0]*b[1] - a[1]*b[0];

        // calculating modulus of that vector
        double normalMod;
        if ((normalMod = sqrt(pow(normalVector[0], 2) +
                              pow(normalVector[1], 2) +
                              pow(normalVector[2], 2) )) == 0.0)
            return -1;

        return 0;
    }
    else{
        // calculating vector normal to diagonals (vector product)
        normal[0] = a[1]*b[2] - a[2]*b[1];
        normal[1] = a[2]*b[0] - a[0]*b[2];
        normal[2] = a[0]*b[1] - a[1]*b[0];

        // calculating modulus of that vector
        double normalMod;
        if ((normalMod = sqrt(pow(normal[0], 2) +
                              pow(normal[1], 2) +
                              pow(normal[2], 2) )) == 0.0)
            return -1;

        // normalize normal
        for (uInt i=0; i<3; i++)
            normal[i] /= normalMod;

        return 0;
    }
}

//-----------------------------------------------------------------------/
// calculates element normal
// check for element validity is not performed (0,0,0 returned)
//-----------------------------------------------------------------------/
vecDbl CMesh::calcNormal(const Elem* element)
{
    // declare and initialize normal vector
    vecDbl normal(3);
    normal[0] = 0.0;
    normal[1] = 0.0;
    normal[2] = 0.0;

    // declare diagonal vectors
    double a[3], b[3];

    // calculate diagonal vectors
    if      (element->nodes.size() == 3){
        for (uInt i=0; i<3; i++) b[i] = element->nodes[2]->coordinates[i]-element->nodes[1]->coordinates[i]; }
    else if (element->nodes.size() == 4){
        for (uInt i=0; i<3; i++) b[i] = element->nodes[3]->coordinates[i]-element->nodes[1]->coordinates[i]; }
    else return normal;
        for (uInt i=0; i<3; i++) a[i] = element->nodes[2]->coordinates[i]-element->nodes[0]->coordinates[i];

    // calculating vector normal to diagonals (vector product)
    normal[0] = a[1]*b[2] - a[2]*b[1];
    normal[1] = a[2]*b[0] - a[0]*b[2];
    normal[2] = a[0]*b[1] - a[1]*b[0];

    // calculating modulus of that vector
    double normalMod;
    if ((normalMod = sqrt(pow(normal[0], 2) +
                          pow(normal[1], 2) +
                          pow(normal[2], 2) )) > 0.0){
        // normalize normal
        for (uInt i=0; i<3; i++)
            normal[i] /= normalMod;
    }

    return normal;
}

//-----------------------------------------------------------------------/
// calculates angle between two element normals
//-----------------------------------------------------------------------/
double CMesh::angleBetweenElems(Elem* elem_1, Elem* elem_2)
{
    // assemble normal vectors
    vecDbl vector_1, vector_2;
    vector_1 = calcNormal(elem_1);
    vector_2 = calcNormal(elem_2);

    // calculate cosine
    double cosine = vector_1[0]*vector_2[0] +
                    vector_1[1]*vector_2[1] +
                    vector_1[2]*vector_2[2] ;

    // check if due to numerical errors we got out of unit circle bounds
    if      (cosine >= 1.0)  return 0.0;
    else if (cosine <= -1.0) return pi;
    else                     return acos(cosine);
}

//-----------------------------------------------------------------------/
// get continuous node branches from node selection
//-----------------------------------------------------------------------/
void CMesh::getNodeBranches(
        const setUInt& nodeSelection,
        vecDeqPNode& nodeBranches)
{
    // extract node pointers from node selection
    setPNode nodes;
    for (setConstUIntIt it=nodeSelection.begin(); it!=nodeSelection.end(); it++){
        // check if node with given ID exists
        mapNodeIt nodeIt = meshNodes.find(*it);
        if (nodeIt != meshNodes.end())
            nodes.insert(&(nodeIt->second));
    }

    // exit if no nodes exist
    if (nodes.empty())
        return;

    // repeat process until all nodes are used
    setUInt usedNodes;
    while (usedNodes.size() != nodes.size()){
        // find start node as first unused node
        Node* startNode = NULL;
        for (setPNodeIt nodeIt=nodes.begin(); nodeIt!=nodes.end(); nodeIt++){
            if (usedNodes.find((*nodeIt)->id) == usedNodes.end()){
                startNode = *nodeIt;
                break;
            }
        }

        // check just in case if startNode was not found
        if (!startNode)
            break;

        // initialize node branch
        deqPNode nodeBranch;
        nodeBranch.push_back(startNode);

        // add starting node to used nodes set
        usedNodes.insert(startNode->id);

        // find adjacent nodes (each for growing in its own direction) by
        // searching elements of starting node for nodes that are also in selection
        Node* nodeBack  = NULL;
        Node* nodeFront = NULL;
        // find adjacent nodes to starting node
        for (setPElemIt elemIt=startNode->elements.begin(); elemIt!=startNode->elements.end(); elemIt++){
            for (uInt nodeIndex=0; nodeIndex<(*elemIt)->nodes.size(); nodeIndex++){
                // check if node exists in node selection
                if (nodes.find((*elemIt)->nodes[nodeIndex]) == nodes.end())
                    continue;
                // check that this node is not the starting node
                if ((*elemIt)->nodes[nodeIndex] == startNode)
                    continue;
                if (!nodeBack )
                    nodeBack  = (*elemIt)->nodes[nodeIndex];
                else if (!nodeFront && (*elemIt)->nodes[nodeIndex]!=nodeBack)
                    nodeFront = (*elemIt)->nodes[nodeIndex];
            }
            // quit searching if both are found
            if (nodeFront)
                break;
        }

        // start growing this branch first at the back then at the front
        if (nodeBack ) growBranch(nodeBranch, nodeBack,  nodes, usedNodes, "back" );
        if (nodeFront) growBranch(nodeBranch, nodeFront, nodes, usedNodes, "front");

        // if first and last node share element, then put first node to the end
        if (nodeBranch.front() != nodeBranch.back() &&
            nodeBranch.size()  >  2                 ){
            Node* nodeA = nodeBranch.front();
            Node* nodeB = nodeBranch.back();
            bool elemFound = false;
            for (setPElemIt elemAIt=nodeA->elements.begin(); elemAIt!=nodeA->elements.end(); elemAIt++){
                for (setPElemIt elemBIt=nodeB->elements.begin(); elemBIt!=nodeB->elements.end(); elemBIt++){
                    if ((*elemAIt)->id == (*elemBIt)->id){
                        nodeBranch.push_back(nodeBranch.front());
                        elemFound = true;
                        break;
                    }
                }
                if (elemFound)
                    break;
            }
        }

        nodeBranches.push_back(nodeBranch);
    }
}

//-----------------------------------------------------------------------/
// grow branch of node line
//-----------------------------------------------------------------------/
void CMesh::growBranch(
        deqPNode& nodeBranch,
        Node* startNode,
        const setPNode& nodes,
        setUInt& usedNodes,
        const string& direction)
{
    // add starting node to branch
    if (usedNodes.find(startNode->id) == usedNodes.end()){
        if (direction == "back") nodeBranch.push_back (startNode);
        else                     nodeBranch.push_front(startNode);
    }

    // add starting node ID to used nodes set
    usedNodes.insert(startNode->id);

    // grow branch by finding adjacent node to first one (initially being starting node)
    Node* firstNode = startNode;
    for(;;){
        // extract elements of node (because node (iterator) might get invalidated)
        setPElem elements = firstNode->elements;

        // search elements for the second node
        bool nodeFound = false;
        for (setPElemIt elemIt=elements.begin(); elemIt!=elements.end(); elemIt++){
            // search nodes of element for second node
            for (uInt nodeIndex=0; nodeIndex<(*elemIt)->nodes.size(); nodeIndex++){
                // extract node
                Node* secondNode = (*elemIt)->nodes[nodeIndex];

                // check that this node
                if (nodes.find(secondNode)         != nodes.end()     && // exists in node selection
                    secondNode                     != firstNode       && // is not the first node
                    usedNodes.find(secondNode->id) == usedNodes.end() ){ // is not already used
                    nodeFound = true;
                    // add node to branch
                    if (direction == "back") nodeBranch.push_back (secondNode);
                    else                     nodeBranch.push_front(secondNode);
                    usedNodes.insert(secondNode->id);
                    // first node is now second node
                    firstNode = secondNode;
                    break;
                }
            }
            if (nodeFound)
                break;
        }
        if (!nodeFound)
            break;
    }
}