/* GTS - Library for the manipulation of triangulated surfaces * Copyright (C) 1999 Stéphane Popinet * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Library General Public License for more details. * * You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include #include #include #include #include "config.h" #ifdef HAVE_GETOPT_H # include #endif /* HAVE_GETOPT_H */ #ifdef HAVE_UNISTD_H # include #endif /* HAVE_UNISTD_H */ #include "gts.h" #ifndef PI # define PI 3.14159265359 #endif static void sphere (gdouble ** f, GtsCartesianGrid g, guint k, gpointer data) { gdouble x, y, z = g.z; guint i, j; for (i = 0, x = g.x; i < g.nx; i++, x += g.dx) for (j = 0, y = g.y; j < g.ny; j++, y += g.dy) f[i][j] = x*x + y*y + z*z; } /* The Clebsch diagonal surface: a smooth cubic surface admitting the * symmetry of the tetrahedron (courtesy Johannes Beigel). */ static void clebsch (gdouble ** f, GtsCartesianGrid g, guint k, gpointer data) { gdouble x, y, z = g.z; guint i, j; gdouble w2 = sqrt (2.0); for (i = 0, x = g.x; i < g.nx; i++, x += g.dx) for (j = 0, y = g.y; j < g.ny; j++, y += g.dy) { /* tetrahedral coordinates */ gdouble p = 1. - z - w2*x; gdouble q = 1. - z + w2*x; gdouble r = 1. + z + w2*y; gdouble s = 1. + z - w2*y; /* symmetric polynomials */ gdouble c1 = p + q + r - s; gdouble c2 = p*p*p + q*q*q + r*r*r - s*s*s; f[i][j] = c2 - c1*c1*c1; } } /* The Barth decic: a degree 10 surface with 345 ordinary double points. * (courtesy Johannes Beigel again). */ static void barth (gdouble ** f, GtsCartesianGrid g, guint k, gpointer data) { gdouble x, y, z = g.z; guint i, j; gdouble t = (1. + sqrt(5.))/2.; for (i = 0, x = g.x; i < g.nx; i++, x += g.dx) for (j = 0, y = g.y; j < g.ny; j++, y += g.dy) { gdouble t4 = t*t*t*t; gdouble p1 = x*x - t4*y*y; gdouble p2 = y*y - t4*z*z; gdouble p3 = z*z - t4*x*x; gdouble p4 = x*x*x*x + y*y*y*y + z*z*z*z - 2.*x*x*y*y - 2.*y*y*z*z - 2.*x*x*z*z; gdouble tmp = x*x + y*y + z*z - 1.; gdouble q1 = tmp*tmp; gdouble tmp1 = x*x + y*y + z*z - (2. - t); gdouble q2 = tmp1*tmp1; f[i][j] = 8.*p1*p2*p3*p4 + (3. + 5.*t)*q1*q2; } } /* Returns isosurface f(x, y, z) = iso with f(x, y, z) = x*x + y*y + z*z. */ int main (int argc, char * argv[]) { int c = 0; GtsCartesianGrid g; GtsSurface * surface; gdouble iso; gboolean verbose = FALSE, tetra = FALSE, dual = FALSE; GtsIsoCartesianFunc func = sphere; if (!setlocale (LC_ALL, "POSIX")) g_warning ("cannot set locale to POSIX"); /* parse options using getopt */ while (c != EOF) { #ifdef HAVE_GETOPT_LONG static struct option long_options[] = { {"dual", no_argument, NULL, 'd'}, {"tetra", no_argument, NULL, 't'}, {"barth", no_argument, NULL, 'b'}, {"clebsch", no_argument, NULL, 'c'}, {"help", no_argument, NULL, 'h'}, {"verbose", no_argument, NULL, 'v'}, { NULL } }; int option_index = 0; switch ((c = getopt_long (argc, argv, "hvbctd", long_options, &option_index))) { #else /* not HAVE_GETOPT_LONG */ switch ((c = getopt (argc, argv, "hvbctd"))) { #endif /* not HAVE_GETOPT_LONG */ case 'd': /* dual */ dual = TRUE; break; case 't': /* tetra */ tetra = TRUE; break; case 'b': /* Barth function */ func = barth; break; case 'c': /* Clebsch function */ func = clebsch; break; case 'v': /* verbose */ verbose = TRUE; break; case 'h': /* help */ fprintf (stderr, "Usage: iso [OPTION] NX NY NZ VAL > file.gts\n" "Compute the isosurface of various functions (default is a sphere).\n" "The size of the cartesian mesh is given by NX, NY and NZ and the isosurface value by VAL.\n" "\n" " -b --barth use the Barth decic function\n" " -c --clebsch use the Clebsch function\n" " -t --tetra use marching tetrahedra (default is marching cubes)\n" " -d --dual use \"dual\" marching tetrahedra\n" " -v --verbose print statistics about the surface\n" " -h --help display this help and exit\n" "\n" "Reports bugs to %s\n", GTS_MAINTAINER); return 0; /* success */ break; case '?': /* wrong options */ fprintf (stderr, "Try `iso --help' for more information.\n"); return 1; /* failure */ } } if (optind >= argc) { /* missing NX */ fprintf (stderr, "iso: missing NX\n" "Try `iso --help' for more information.\n"); return 1; /* failure */ } g.nx = atoi (argv[optind]); if (optind + 1 >= argc) { /* missing NY */ fprintf (stderr, "iso: missing NY\n" "Try `iso --help' for more information.\n"); return 1; /* failure */ } g.ny = atoi (argv[optind + 1]); if (optind + 2 >= argc) { /* missing NZ */ fprintf (stderr, "iso: missing NZ\n" "Try `iso --help' for more information.\n"); return 1; /* failure */ } g.nz = atoi (argv[optind + 2]); if (optind + 3 >= argc) { /* missing VAL */ fprintf (stderr, "iso: missing VAL\n" "Try `iso --help' for more information.\n"); return 1; /* failure */ } iso = atof (argv[optind + 3]); /* interval is [-10:10][-10:10][-10:10] */ g.x = -10.0; g.dx = 20./(gdouble) (g.nx - 1); g.y = -10.0; g.dy = 20./(gdouble) (g.ny - 1); g.z = -10.0; g.dz = 20./(gdouble) (g.nz - 1); surface = gts_surface_new (gts_surface_class (), gts_face_class (), gts_edge_class (), gts_vertex_class ()); if (tetra) gts_isosurface_tetra (surface, g, func, NULL, iso); else if (dual) gts_isosurface_tetra_bcl (surface, g, func, NULL, iso); else gts_isosurface_cartesian (surface, g, func, NULL, iso); if (verbose) gts_surface_print_stats (surface, stderr); gts_surface_write (surface, stdout); return 0; }