lib/Generate/trimesh.c

   1 /*
   2     This file is part of the FElt finite element analysis package.
   3     Copyright (C) 1993 Jason I. Gobat and Darren C. Atkinson
   4
   5     This program is free software; you can redistribute it and/or modify
   6     it under the terms of the GNU General Public License as published by
   7     the Free Software Foundation; either version 2 of the License, or
   8     (at your option) any later version.
   9
  10     This program is distributed in the hope that it will be useful,
  11     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13     GNU General Public License for more details.
  14
  15     You should have received a copy of the GNU General Public License
  16     along with this program; if not, write to the Free Software
  17     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18 */
  19
  20 /****************************************************************************
  21  *
  22  * File:        triangle.c
  23  *
  24  ***************************************************************************/
  25
  26 # include <math.h>
  27 # include <stdio.h>
  28 # include "allocate.h"
  29 # include "error.h"
  30 # include "fe.h"
  31 # include "objects.h"
  32 # include "mesh.h"
  33 # include "triangle.h"
  34
  35 typedef double dbl_pair [2];
  36 typedef int    triple_int [3];
  37
  38 static double PolygonArea (vcl, n)
  39    double  (*vcl) [2];
  40    int       n;
  41 {
  42    int       i;
  43    double    A;
  44
  45    A = 0.0;
  46    for (i = 0 ; i < n - 1 ; i++)
  47       A += (vcl [i][0]*vcl [i+1][1] - vcl [i+1][0]*vcl [i][1]);
  48
  49    A += (vcl [n-1][0]*vcl [0][1] - vcl [0][0]*vcl [n-1][1]);
  50
  51    return 0.5*A;
  52 }
  53
  54 /****************************************************************************
  55  *
  56  * Function:    GenerateTriMesh
  57  *
  58  * Description: a procedure to interface to Triangle and generate a mesh
  59  *              of triangular elements.
  60  *
  61  ****************************************************************************/
  62
  63 unsigned GenerateTriMesh (trimesh,element,node,numelts,numnodes,bnode,belement)
  64    TriMesh      trimesh;
  65    Element      **element;
  66    Node         **node;
  67    unsigned     *numelts;
  68    unsigned     *numnodes;
  69    unsigned     bnode;
  70    unsigned     belement;
  71 {
  72    unsigned     i,j;
  73    unsigned     count;
  74    int          n, m;
  75    int          ne, nn;
  76    int          status;
  77    struct triangulateio   in, out;
  78    int          npoints;
  79    double       x, y;
  80    int          start;
  81    double       A, Atot;
  82    double       Acrit;
  83    char         opts [32];
  84
  85    if (trimesh -> numcurves <= 0) {
  86       error ("must have at least a boundary curve to generate a TriMesh");
  87       return 1;
  88    }
  89
  90    if (trimesh -> definition -> numnodes != 3) {
  91       error ("TriMesh generation requires three node elements");
  92       return 1;
  93    }
  94
  95
  96    npoints = 0;
  97
  98    Atot = 0.0;
  99
 100    for (i = 0 ; i < trimesh -> numcurves ; i++) {
 101       if (trimesh -> curves [i] -> numvc < 3) {
 102          error ("each curve must have at least 3 points");
 103          return 1;
 104       }
 105
 106       npoints += trimesh -> curves [i] -> numvc;
 107
 108       A = PolygonArea(trimesh -> curves [i] -> vcl,
 109                       trimesh -> curves [i] -> numvc);
 110
 111       if (i == 0 && A < 0) {
 112          error ("main boundary points must be in CCW order");
 113          return 1;
 114       }
 115       else if (i > 0 && A > 0) {
 116          error ("hole boundary points must be in CW order");
 117          return 1;
 118       }
 119
 120       Atot += A;
 121    }
 122
 123    Acrit = trimesh -> alpha * Atot / trimesh -> target;
 124
 125    if (npoints <= 0) {
 126       error ("nothing to generate");
 127       return 1;
 128    }
 129
 130    in.numberofpoints = npoints;
 131    in.numberofpointattributes = 0;
 132    in.pointlist = (double *) malloc(in.numberofpoints * 2 * sizeof(double));
 133    in.pointattributelist = NULL;
 134    in.pointmarkerlist = (int *) malloc(sizeof(int) * in.numberofpoints);
 135    for (i = 0 ; i < in.numberofpoints ; i++)
 136        in.pointmarkerlist [i] = 1;
 137
 138    in.numberofholes = trimesh -> numcurves - 1;
 139    if (in.numberofholes)
 140       in.holelist = (double *) malloc(2 * in.numberofholes * sizeof(double));
 141    else
 142       in.holelist = NULL;
 143
 144    in.numberofregions = 0;
 145    in.numberofsegments = npoints;
 146    in.segmentlist = (int *) malloc(2 * in.numberofsegments * sizeof(int));
 147    in.segmentmarkerlist = (int *) malloc(sizeof(int) * in.numberofsegments);
 148    for (i = 0 ; i < in.numberofsegments ; i++)
 149        in.segmentmarkerlist [i] = 1;
 150
 151    n = 0;
 152    m = 0;
 153    for (i = 0 ; i < trimesh -> numcurves ; i++) {
 154       start = m;
 155       for (j = 0 ; j < trimesh -> curves [i] -> numvc ; j++) {
 156          in.pointlist [n]     = trimesh -> curves [i] -> vcl [j][0];
 157          in.pointlist [n + 1] = trimesh -> curves [i] -> vcl [j][1];
 158
 159          in.segmentlist [n]     = m;
 160          in.segmentlist [n + 1] = m + 1;
 161
 162          n += 2;
 163          m += 1;
 164       }
 165
 166       in.segmentlist [n - 1] = start;
 167    }
 168
 169    for (i = 0 ; i < in.numberofholes ; i++) {
 170       x = 0.0;
 171       y = 0.0;
 172
 173       for (j = 0 ; j < trimesh -> curves [i+1] -> numvc ; j++) {
 174          x += trimesh -> curves [i+1] -> vcl [j][0];
 175          y += trimesh -> curves [i+1] -> vcl [j][1];
 176       }
 177
 178       x /= trimesh -> curves [i+1] -> numvc;
 179       y /= trimesh -> curves [i+1] -> numvc;
 180
 181       in.holelist [2*i]     = x;
 182       in.holelist [2*i + 1] = y;
 183    }
 184
 185    out.pointlist                  = (double *) NULL;
 186    out.pointattributelist         = (double *) NULL;
 187    out.pointmarkerlist            = (int *) NULL;
 188    out.trianglelist               = (int *) NULL;
 189    out.numberoftriangleattributes = 0;
 190    out.triangleattributelist      = (double *) NULL;
 191    out.segmentmarkerlist          = (int *) NULL;
 192    out.segmentlist                = (int *) NULL;
 193
 194    sprintf (opts, "Qpqza%f", Acrit);
 195
 196    triangulate(opts, &in, &out, NULL);
 197
 198    ne = out.numberoftriangles;
 199    nn = out.numberofpoints;
 200    if (ne <= 0 || nn <= 0) {
 201       error ("nothing to generate");
 202       return 1;
 203    }
 204
 205         /*
 206          * allocate some memory to hold everything that we will generate
 207          */
 208
 209    if (!(*node = Allocate(Node, nn)))
 210       Fatal ("allocation error in TriMesh generation");
 211
 212    UnitOffset (*node);
 213
 214    for (i = 1 ; i <= nn ; i++) {
 215       if (!((*node) [i] = CreateNode (0)))
 216          Fatal ("allocation error in TriMesh generation");
 217    }
 218
 219    if (!(*element = Allocate(Element, ne)))
 220       Fatal ("allocation error in TriMesh generation");
 221
 222    UnitOffset (*element);
 223
 224    for (i = 1 ; i <= ne ; i++) {
 225       if (!((*element) [i] = CreateElement (0, trimesh -> definition)))
 226          Fatal ("allocation error in TriMesh generation");
 227    }
 228
 229         /*
 230          * generate all the nodes
 231          */
 232
 233    for (i = 1 ; i <= nn ; i++) {
 234
 235       (*node) [i] -> number = i + bnode;
 236       (*node) [i] -> x = out.pointlist [2*i - 2];
 237       (*node) [i] -> y = out.pointlist [2*i - 1];
 238       (*node) [i] -> z = 0;
 239
 240    }
 241
 242         /*
 243          * attach all the elements to the nodes
 244          */
 245
 246    for (i = 1 ; i <= ne ; i++) {
 247
 248       (*element) [i] -> number = i + belement;
 249       (*element) [i] -> node [1] = (*node) [out.trianglelist [3*(i - 1) + 0] + 1];
 250       (*element) [i] -> node [2] = (*node) [out.trianglelist [3*(i - 1) + 1] + 1];
 251       (*element) [i] -> node [3] = (*node) [out.trianglelist [3*(i - 1) + 2] + 1];
 252
 253    }
 254
 255    *numnodes = nn;
 256    *numelts  = ne;
 257
 258    Deallocate(in.pointlist);
 259    Deallocate(in.pointmarkerlist);
 260    Deallocate(in.holelist);
 261    Deallocate(in.segmentlist);
 262    Deallocate(in.segmentmarkerlist);
 263    Deallocate(out.pointlist);
 264    Deallocate(out.pointattributelist);
 265    Deallocate(out.pointmarkerlist);
 266    Deallocate(out.trianglelist);
 267    Deallocate(out.triangleattributelist);
 268    Deallocate(out.segmentmarkerlist);
 269    Deallocate(out.segmentlist);
 270
 271    return 0;
 272 }