source: GTP/trunk/Lib/Geom/shared/GTGeometry/src/libs/vmi/src/mesh.cpp @ 2090

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <string.h>
#include <math.h>

#include "../include/mesh.h"
#include "../include/area.h"
#include "../include/global.h"

using namespace VMI;

#define OUTPUT_SEVERAL_GROUPS

///////////////////////////////////////////////////////////////////////////////
// list of integers
void VMI::printItemList(int *list, int n) {
    int i;
    
    for (i=0; i<n; i++)
        printf("%d ",list[i]);

    printf("\n");
}

int VMI::findItem(int *list, int n, int item) {
    int i;

    for (i=0; i<n; i++)
        if (list[i] == item) return TRUE;
        return FALSE;
}

void VMI::addItem(int *list, int *n, int item) {

    if (!findItem(list, *n, item))
        list[(*n)++] = item;
}

void VMI::delItem(int *list, int *n, int item) {
    int i, j;
    
    for (i=0; i<*n; i++) {
        if (list[i] == item) {
            // delete it
            for (j =i + 1 ; j<*n; j++)
                list[j - 1] = list[j];
            (*n)--;
            i--; // Delete all ocurrencies of an item
        }
    }
}
///////////////////////////////////////////////////////////////////////////////

int VMI::findEdge(Edge *e, GLuint num, GLuint _u, GLuint _v) {
    GLuint i, u, v;
    int found = -1; 

    for (i=0; i<num; i++) {

        if (e[i].enable == TRUE) {

            u = e[i].u;
            v = e[i].v;

            if (((u == _u) && (v == _v)) ||
                ((u == _v) && (v == _u))) {
                found = i;
                break;
            }
        }
    }

    return found;
}

// Adds a new edge at the end of the mesh edge list
Edge *VMI::addEdge(Edge *list, int *n, int u, int v, int *pos) {
          Edge *newList;

                // Reallocate memory for the new edge
    newList = (Edge *)realloc(list, ((*n) + 1) * sizeof(Edge));

                // New vertex indices
                newList[*n].u = u;
    newList[*n].v = v;

                // This edge has no triangle adjancency
                newList[*n].numTriangles = 0;
                newList[*n].triangles = NULL;
                
                // Enable the new edge
                newList[*n].enable = TRUE;

                // New edge position
                *pos = *n;

                // New number of edges
          (*n)++;

                return newList;
}

GLdouble VMI::computeTriangleArea(Vertex *vertices, Triangle *t) {
    GLuint i;
    GLdouble v0[3], v1[3], v2[3], a;
    
    // Compute triangle area
    i = t->indices[0];
    v0[0] = vertices[i].x;
    v0[1] = vertices[i].y;
    v0[2] = vertices[i].z;

    i = t->indices[1];
    v1[0] = vertices[i].x;
    v1[1] = vertices[i].y;
    v1[2] = vertices[i].z;

    i = t->indices[2];
    v2[0] = vertices[i].x;
    v2[1] = vertices[i].y;
    v2[2] = vertices[i].z;

    a = triangleArea(v0, v1, v2);
    
    //printf("Triangle area: %f\n", a)
    return a;  
}

void VMI::computeTriangleNormal(Vertex *vertices, Triangle *t) {
    GLuint i;
    GLfloat x0, y0, z0,
            x1, y1, z1,
            x2, y2, z2;
    GLfloat x10, y10, z10,
            x12, y12, z12,
            cpx, cpy, cpz, r;

    i = t->indices[0];
    x0 = vertices[i].x;
    y0 = vertices[i].y;
    z0 = vertices[i].z;

    i = t->indices[1];
    x1 = vertices[i].x;
    y1 = vertices[i].y;
    z1 = vertices[i].z;

    i = t->indices[2];
    x2 = vertices[i].x;
    y2 = vertices[i].y;
    z2 = vertices[i].z;

    /* Compute edge vectors */
    x10 = x1 - x0;
    y10 = y1 - y0;
    z10 = z1 - z0;
    x12 = x1 - x2;
    y12 = y1 - y2;
    z12 = z1 - z2;

    /* Compute the cross product */
    cpx = (z10 * y12) - (y10 * z12);
    cpy = (x10 * z12) - (z10 * x12);
    cpz = (y10 * x12) - (x10 * y12);

    /* Normalize the result to get the unit-length facet normal */
    r = sqrt(cpx * cpx + cpy * cpy + cpz * cpz);

    t->normal[0] = cpx / r;
    t->normal[1] = cpy / r;
    t->normal[2] = cpz / r;
}

Mesh *VMI::initMesh(GLMmodel* pmodel) {

    GLuint i, j, curGroup, v1, v2, v3, n, m, t;
    int e;
    Mesh *mesh;
    GLMgroup *group;

    mesh = (Mesh *)malloc (sizeof(Mesh));
    
    if (mesh == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }

    mesh->vertices = (Vertex *)malloc (sizeof(Vertex) * pmodel->numvertices);
    
    if (mesh->vertices == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }

    mesh->triangles = (Triangle *)malloc (sizeof(Triangle) * pmodel->numtriangles);
    
    if (mesh->triangles == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }
    
    printf("Adding vertices...");
    for (i=1; i<=pmodel->numvertices; i++) { // Vertices start at 1
        mesh->vertices[i - 1].x = pmodel->vertices[3 * i + 0];
        mesh->vertices[i - 1].y = pmodel->vertices[3 * i + 1];
        mesh->vertices[i - 1].z = pmodel->vertices[3 * i + 2];
        mesh->vertices[i - 1].numTriangles = 0;
        mesh->vertices[i - 1].triangles = NULL;
        mesh->vertices[i - 1].enable = GL_TRUE;
    }
    printf("Ok\n");
    mesh->numVertices = pmodel->numvertices;
    mesh->currentNumVertices = pmodel->numvertices;

    printf("Adding triangles...");

#ifdef OUTPUT_SEVERAL_GROUPS
    group = pmodel->groups;
    curGroup=0; // current group
    while(group) {
        
        for (i=0; i<group->numtriangles; i++) {
            
            t = group->triangles[i];

            mesh->triangles[t].id = t;
            mesh->triangles[t].group = curGroup; // set group
            mesh->triangles[t].indices[0] = pmodel->triangles[t].vindices[0] - 1;
            mesh->triangles[t].indices[1] = pmodel->triangles[t].vindices[1] - 1;
            mesh->triangles[t].indices[2] = pmodel->triangles[t].vindices[2] - 1;
            mesh->triangles[t].area = computeTriangleArea(mesh->vertices, &mesh->triangles[t]);
            //printf("\n%d a: %f",i , mesh->triangles[i].area);
            computeTriangleNormal(mesh->vertices, &mesh->triangles[t]);
            
            mesh->triangles[t].saliency = 0.0;
            
            mesh->triangles[t].enable = GL_TRUE;
            
            for (j=0; j<3; j++) {
                // Adding triangle i adjacent to 3 vertices
                v1 = mesh->triangles[t].indices[j];
                
                // Reallocate memory for the new adjacent triangle
                mesh->vertices[v1].triangles = 
                                        (GLuint *)realloc(mesh->vertices[v1].triangles, (mesh->vertices[v1].numTriangles + 1) * sizeof(GLuint));
                addItem((int *)mesh->vertices[v1].triangles, (int *)&mesh->vertices[v1].numTriangles, t);
            }
        }
        curGroup++;
        group = group->next;
    }
#else // One single group
    for (i=0; i<pmodel->numtriangles; i++) {
        mesh->triangles[i].id = i;
        mesh->triangles[i].indices[0] = pmodel->triangles[i].vindices[0] - 1;
        mesh->triangles[i].indices[1] = pmodel->triangles[i].vindices[1] - 1;
        mesh->triangles[i].indices[2] = pmodel->triangles[i].vindices[2] - 1;
        mesh->triangles[i].area = computeTriangleArea(mesh->vertices, &mesh->triangles[i]);
        //printf("\n%d a: %f",i , mesh->triangles[i].area);
        computeTriangleNormal(mesh->vertices, &mesh->triangles[i]);

        mesh->triangles[i].saliency = 0.0;

        mesh->triangles[i].enable = GL_TRUE;
        
        for (j=0; j<3; j++) {
            // Adding triangle i adjacent to 3 vertices
            v1 = mesh->triangles[i].indices[j];
   
            // Reallocate memory for the new adjacent triangle
            mesh->vertices[v1].triangles = 
                            (GLuint *)realloc(mesh->vertices[v1].triangles, (mesh->vertices[v1].numTriangles + 1) * sizeof(GLuint));
            addItem((int *)mesh->vertices[v1].triangles, (int *)&mesh->vertices[v1].numTriangles, i);
        }
    }
    printf("\n");
#endif
    printf("Ok\n");

    mesh->numTriangles = pmodel->numtriangles;
    mesh->currentNumTriangles = pmodel->numtriangles;
    printf("Number of triangles: %d\n", mesh->numTriangles);

    mesh->edges = (Edge *)malloc (sizeof(Edge) * mesh->numTriangles * 3);  // E = 3 T / 2
    
    if (mesh->edges == NULL) {
        fprintf(stderr, "Error allocating memory\n");
        exit(1);
    }
    
    printf("Adding edges...");
    n = 0;

    for (i=0; i<mesh->numTriangles; i++) {
        t = 0;
        v1 = mesh->triangles[i].indices[0];
        v2 = mesh->triangles[i].indices[1];
        v3 = mesh->triangles[i].indices[2];
        
        /////////////////////////////////////////////////////////////////////////////////
                // edge (v1,v2)
        if ((e = findEdge(mesh->edges, n, v1, v2)) == -1) {
            mesh->edges[n].u = v1;
            mesh->edges[n].v = v2;
                                                mesh->edges[n].triangles    = NULL;
                        mesh->edges[n].numTriangles = 0;
            mesh->edges[n].enable = GL_TRUE;
            m = n;
            n++;
        } else m = e;
                // Reallocate memory for the new adjacent triangle
        mesh->edges[m].triangles = 
                        (GLuint *)realloc(mesh->edges[m].triangles, (mesh->edges[m].numTriangles + 1) * sizeof(GLuint));
        // Adding triangle i adjacent to edge m
        addItem((int *)mesh->edges[m].triangles, (int *)&mesh->edges[m].numTriangles, i);
        //printf("n:%d i:%d\n", m, i);
        
        // Adding edge m adjacent to triangle i
        addItem((int *)mesh->triangles[i].edges, (int *)&t, m);

        /////////////////////////////////////////////////////////////////////////////////
        // edge (v2,v3)
        if (( e = findEdge(mesh->edges, n, v2, v3)) == -1) {
            mesh->edges[n].u = v2;
            mesh->edges[n].v = v3;
                                                mesh->edges[n].triangles    = NULL;
                        mesh->edges[n].numTriangles = 0;
            mesh->edges[n].enable = GL_TRUE;
            m = n;
            n++;
        } else m = e;
                // Reallocate memory for the new adjacent triangle
        mesh->edges[m].triangles = 
                        (GLuint *)realloc(mesh->edges[m].triangles, (mesh->edges[m].numTriangles + 1) * sizeof(GLuint));
        // Adding triangle i adjacent to edge m
        addItem((int *)mesh->edges[m].triangles, (int *)&mesh->edges[m].numTriangles, i);
        //printf("n:%d i:%d\n", m, i);
        
        // Adding edge m adjacent to triangle i
        addItem((int *)mesh->triangles[i].edges, (int *)&t, m);

        /////////////////////////////////////////////////////////////////////////////////
        // edge (v3,v1)
        if ((e = findEdge(mesh->edges, n, v3, v1)) == -1) {
            mesh->edges[n].u = v3;
            mesh->edges[n].v = v1;
                                                mesh->edges[n].triangles    = NULL;
                        mesh->edges[n].numTriangles = 0;
            mesh->edges[n].enable = GL_TRUE;
            m = n;
            n++;
        } else m = e;
                // Reallocate memory for the new adjacent triangle
        mesh->edges[m].triangles = 
                        (GLuint *)realloc(mesh->edges[m].triangles, (mesh->edges[m].numTriangles + 1) * sizeof(GLuint));
        // Adding triangle i adjacent to edge m
        addItem((int *)mesh->edges[m].triangles, (int *)&mesh->edges[m].numTriangles, i);
        //printf("n:%d i:%d\n", m, i);
        
        // Adding edge m adjacent to triangle i
        addItem((int *)mesh->triangles[i].edges, (int *)&t, m);
    }
    printf("Ok\n");
    mesh->numEdges = n;

    return mesh;
}

void VMI::printMesh(Mesh *mesh) {
    GLuint i, j;

    printf("Vertices (%d)\n", mesh->numVertices);
    printf("------------------------\n");
    
    for (i=0; i<mesh->numVertices; i++) {
        printf("v%d (%f,%f,%f)\n", i, mesh->vertices[i].x, mesh->vertices[i].y,
                                     mesh->vertices[i].z);
        printf("  t(%d)[", mesh->vertices[i].numTriangles);
        if (mesh->vertices[i].triangles != NULL) {
            printf("%d", mesh->vertices[i].triangles[0]);
            for (j=1; j<mesh->vertices[i].numTriangles; j++) {
                printf(",%d", mesh->vertices[i].triangles[j]);
            }
        }
        printf("]\n");
    }

    printf("Edges (%d)\n", mesh->numEdges);
    printf("------------------------\n");
    for (i=0; i<mesh->numEdges; i++) {
        printf("e%d (%d,%d) %d\n", i, mesh->edges[i].u, mesh->edges[i].v, mesh->edges[i].enable);
        printf("  t(%d)[", mesh->edges[i].numTriangles);
        if (mesh->edges[i].triangles != NULL) {
            printf("%d", mesh->edges[i].triangles[0]);
            for (j=1; j<mesh->edges[i].numTriangles; j++) {
                printf(",%d", mesh->edges[i].triangles[j]);
            }
                                }
        printf("]\n");
    }

    printf("Triangles (%d)\n", mesh->numTriangles);
    printf("------------------------\n");
    for (i=0; i<mesh->numTriangles; i++) {
        printf("t%d (%d,%d,%d) %d\n", i, mesh->triangles[i].indices[0], mesh->triangles[i].indices[1],
                                     mesh->triangles[i].indices[2],  mesh->triangles[i].enable);
        printf("  e(3)[");
        printf("%d", mesh->triangles[i].edges[0]);
        for (j=1; j<3; j++) {
            printf(",%d", mesh->triangles[i].edges[j]);
        }
        printf("]\n");

        printf("  n(3)[");
        printf("%f", mesh->triangles[i].normal[0]);
        for (j=1; j<3; j++) {
            printf(",%f", mesh->triangles[i].normal[j]);
        }
        printf("]\n");
       
        printf("  area: %f\n", mesh->triangles[i].area);
    }
}

void VMI::deleteMesh(Mesh *mesh) {
        GLuint i;

        if (NULL != mesh) {
                for (i=0; i<mesh->numVertices; i++) {
                        if (mesh->vertices[i].triangles != NULL)
                                free(mesh->vertices[i].triangles);

                }
                if (mesh->vertices != NULL) free(mesh->vertices);
                mesh->numVertices = 0;

                for (i=0; i<mesh->numEdges; i++) {
                        if (mesh->edges[i].triangles != NULL)
                                free(mesh->edges[i].triangles);
                }
                if (mesh->edges != NULL) free(mesh->edges);
                mesh->numEdges = 0;

                if (mesh->triangles != NULL) free(mesh->triangles);
                mesh->numTriangles = 0;

                free(mesh);

                mesh = NULL;
        }
}

int *VMI::trianglesAdjToEdge(Mesh *mesh, int e, int *n) {
    int i, n0, n1, n2, num = 0, t, 
        u = mesh->edges[e].u, 
        v = mesh->edges[e].v;
    int *triangles = NULL;

        if ((mesh->vertices[u].numTriangles + 
         mesh->vertices[v].numTriangles) > 0)
                 triangles = (int *)malloc((mesh->vertices[u].numTriangles + 
                                    mesh->vertices[v].numTriangles) * sizeof(int));

    for (i=0; i<(int)mesh->vertices[u].numTriangles; i++) {

        t = mesh->vertices[u].triangles[i];
        
        if (mesh->triangles[t].enable == TRUE) {
            
            n0 = mesh->triangles[t].indices[0];
            n1 = mesh->triangles[t].indices[1];
            n2 = mesh->triangles[t].indices[2];
            
            if ((n0 == v) || (n1 == v) || (n2 == v) ||
                (n0 == u) || (n1 == u) || (n2 == u))
                
                addItem(triangles, &num, t);
        }
    }

    for (i=0; i<(int)mesh->vertices[v].numTriangles; i++) {

        t = mesh->vertices[v].triangles[i];
        
        if (mesh->triangles[t].enable == TRUE) {
            
            n0 = mesh->triangles[t].indices[0];
            n1 = mesh->triangles[t].indices[1];
            n2 = mesh->triangles[t].indices[2];
            
            if ((n0 == v) || (n1 == v) || (n2 == v) ||
                (n0 == u) || (n1 == u) || (n2 == u))
                
                addItem(triangles, &num, t);
        }
    }
    *n = num;

    return triangles;
}

int *VMI::verticesAdjToVertex(Mesh *mesh, int v, int *n) {
    int *vertices = NULL;
    int i, t, n0, n1, n2;
    int num = 0;

        if (mesh->vertices[v].numTriangles > 0)
                vertices = (int*)malloc(sizeof(int) * mesh->vertices[v].numTriangles * 2);

    for (i=0; i<(int)mesh->vertices[v].numTriangles; i++) {
        t = mesh->vertices[v].triangles[i];

        n0 = mesh->triangles[t].indices[0];
        n1 = mesh->triangles[t].indices[1];
        n2 = mesh->triangles[t].indices[2];

        if (n0 != v) addItem(vertices, &num, n0);
        if (n1 != v) addItem(vertices, &num, n1);
        if (n2 != v) addItem(vertices, &num, n2);
    }
    *n = num;

    return vertices;
}

int *VMI::edgesAdjToVertex(Mesh *mesh, int v, int *n) {
    int *edges = NULL;
    int i, t, e;
#ifndef USE_EDGE_ADJACENCY
    int n0, n1, n2;
#else
    int j;
#endif
    int num = 0;

        if (mesh->vertices[v].numTriangles > 0) 
                edges = (int*)malloc(sizeof(int) * mesh->vertices[v].numTriangles * 2);

    for (i=0; i<(int)mesh->vertices[v].numTriangles; i++) {

        t = mesh->vertices[v].triangles[i];

#ifdef USE_EDGE_ADJACENCY
        for (j=0; j<3; j++) {
            e = mesh->triangles[t].edges[j];
            
            if ((int)mesh->edges[e].u == v || (int)mesh->edges[e].v == v)
                
                addItem(edges, &num, e);
        }
#else
        n0 = mesh->triangles[t].indices[0];
        n1 = mesh->triangles[t].indices[1];
        n2 = mesh->triangles[t].indices[2];
        
        if (((n0 == v) || (n1 == v)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n0, n1)) != -1)) addItem(edges, &num, e);
 
        if (((n1 == v) || (n2 == v)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n1, n2)) != -1)) addItem(edges, &num, e);
        
        if (((n0 == v) || (n2 == v)) &&
            ((e = findEdge(mesh->edges, mesh->numEdges, n0, n2)) != -1)) addItem(edges, &num, e);

#endif
    }
    *n = num;

    return edges;
}

int *VMI::edgesAdjToVertices(Mesh *mesh, int *vertices, int numVertices, int *n) {
    int *edges = NULL;
    int i, j, num = 0;
    int *list = NULL, numEdges;

    if (numVertices > 0) {
        // Add the first
        list = edgesAdjToVertex(mesh, vertices[0], &numEdges);

        // Allocate memory
        edges = (int *)malloc(numEdges * sizeof(int));

        memcpy(edges, list, sizeof(int) * numEdges);

        free(list);

        num = numEdges;

        // Add the rest
        for (i=1; i<numVertices; i++) {

            list = edgesAdjToVertex(mesh, vertices[i], &numEdges);

            // Reallocate memory
                        if (numEdges > 0) {

                edges = (int *)realloc(edges, (num + numEdges + 1) * sizeof(int));
                
                for (j=0; j<numEdges; j++)
                    addItem(edges, &num, list[j]);
                        }
            free(list);
        }
    }
    *n = num;

    return edges;
}

GLboolean VMI::findVertex(GLfloat *vertices, GLuint num, GLfloat x, GLfloat y, GLfloat z, int *pos) {
    GLuint i;
    GLfloat _x, _y, _z;
    GLboolean found = GL_FALSE; 
    
    *pos = -1;

    for (i=1; i<num; i++) { // i=0

        _x = vertices[3 * i + 0];
        _y = vertices[3 * i + 1];
        _z = vertices[3 * i + 2];

        if (((_x == x) && (_y == y) && (_z == z))) {
            found = GL_TRUE;
            *pos = (int)i;
            break;
        }
    }

    return found;
}

// Adds a new vertex at the end of the mesh vertex list
Vertex *VMI::addVertex(Vertex *list, int *n, float x, float y, float z, int *pos) {
          Vertex *newList;

                // Reallocate memory for the new vertex
    newList = (Vertex *)realloc(list, ((*n) + 1) * sizeof(Vertex));

                // New Vertex coordinates
                newList[*n].x = x;
    newList[*n].y = y;
    newList[*n].z = z;

                // This vertex has no triangle adjancency
                newList[*n].numTriangles = 0;
                newList[*n].triangles = NULL;
        
    // Enable the new vertex
                newList[*n].enable = TRUE;

                // New vertex position
                *pos = *n;

                // New number of vertices
          (*n)++;
                
                return newList;
}



void VMI::updateModel(GLMmodel* pmodel, Mesh *mesh, GLuint numVertices, GLuint numTriangles) {
    GLuint i , v1, v2 ,v3, numV = 1, numT = 0;
    int pos;
    GLfloat x, y, z;
    
    if (pmodel->vertices != NULL) free(pmodel->vertices);
    pmodel->vertices = (GLfloat *)calloc ((numVertices + 1) * 3,  sizeof(GLfloat));

    if (pmodel->triangles != NULL) free(pmodel->triangles);
    pmodel->triangles = (GLMtriangle *)malloc (numTriangles * sizeof(GLMtriangle));

    for (i=0; i<mesh->numTriangles; i++) {
        if (mesh->triangles[i].enable) {
            v1 = mesh->triangles[i].indices[0];
            v2 = mesh->triangles[i].indices[1];
            v3 = mesh->triangles[i].indices[2];
            
            x = mesh->vertices[v1].x;
            y = mesh->vertices[v1].y;
            z = mesh->vertices[v1].z;
            
            if (!findVertex(pmodel->vertices, numV, x, y ,z, &pos)) {
                pmodel->triangles[numT].vindices[0] = numV;
                
                pmodel->vertices[3 * numV + 0] = x;
                pmodel->vertices[3 * numV + 1] = y;
                pmodel->vertices[3 * numV + 2] = z;
                numV++;
            } else pmodel->triangles[numT].vindices[0] = (GLuint)pos;
            
            x = mesh->vertices[v2].x;
            y = mesh->vertices[v2].y;
            z = mesh->vertices[v2].z;
            
            if (!findVertex(pmodel->vertices, numV, x, y ,z, &pos)) {
                pmodel->triangles[numT].vindices[1] = numV;
                
                pmodel->vertices[3 * numV + 0] = x;
                pmodel->vertices[3 * numV + 1] = y;
                pmodel->vertices[3 * numV + 2] = z;
                numV++;
            } else pmodel->triangles[numT].vindices[1] = (GLuint)pos;
            
            x = mesh->vertices[v3].x;
            y = mesh->vertices[v3].y;
            z = mesh->vertices[v3].z;
            
            if (!findVertex(pmodel->vertices, numV, x, y ,z, &pos)) {
                pmodel->triangles[numT].vindices[2] = numV;
                
                pmodel->vertices[3 * numV + 0] = x;
                pmodel->vertices[3 * numV + 1] = y;
                pmodel->vertices[3 * numV + 2] = z;
                numV++;
            } else pmodel->triangles[numT].vindices[2] = (GLuint)pos;
            numT++;
        }
    }
    pmodel->numvertices  = numV - 1;
    pmodel->numtriangles = numT;
}

void VMI::saveModel(char *filename, GLMmodel* pmodel, Mesh *mesh) {
    GLuint i, j, g, numT = 0;
    GLMgroup *group;

#ifdef OUTPUT_SEVERAL_GROUPS
    // The output model maintains groups

    // Empty all groups
    group = pmodel->groups;
    while(group) {
        // Free memory for triangles
        free(group->triangles);

        // Reset group
        group->triangles = NULL;
        group->numtriangles = 0;

        group = group->next;
    }
    // Fill all groups
    for (i=0; i<mesh->numTriangles; i++) {

        if (mesh->triangles[i].enable) {
            g = mesh->triangles[i].group;

            // Get group
            group = pmodel->groups;
            for (j=0; j <g; j++)
                group = group->next;

            // Reallocate memory for triangles
            group->triangles = (GLuint *)realloc(group->triangles, (group->numtriangles + 1) * sizeof(GLuint));
            // Add triangle into group
            addItem((int *)group->triangles, (int *)&group->numtriangles, numT);

            numT++;
        }
    }
#else
    // The output model contains only one group
    pmodel->numgroups = 1;

    if (pmodel->groups->name != NULL) free(pmodel->groups->name);
    pmodel->groups->name = strdup("default");
    pmodel->groups->material = 0;
    pmodel->groups->next = NULL;
    pmodel->groups->numtriangles = mesh->currentNumTriangles;

    if (pmodel->groups->triangles != NULL) free(pmodel->groups->triangles);
    pmodel->groups->triangles = (GLuint*)malloc(sizeof(GLuint) * mesh->currentNumTriangles);

    for (i=0; i<mesh->currentNumTriangles; i++)
        pmodel->groups->triangles[i] = i;
#endif

    updateModel(pmodel, mesh, mesh->currentNumTriangles * 3, mesh->currentNumTriangles);

    printf("glmWriteOBJ...");
    
    glmWriteOBJ(pmodel, filename, GLM_NONE);
    printf("Ok\n");
}