<p><b>dwj07@fsu.edu</b> 2011-09-01 11:45:41 -0600 (Thu, 01 Sep 2011)</p><p><br>
Adding a tools branch, with mpas_draw and points-mpas included.<br>
<br>
both are c++ programs, to be used with mpas.<br>
<br>
points-mpas:<br>
This program reads a point set, and a triangulation of that point set to make<br>
a full spherical grid for use with mpas.<br>
A README is included with more detail<br>
<br>
mpas_draw:<br>
This program reads in an mpas grid/input/output/restart file and allows you<br>
to visualize the fields included. It is particularly useful when exploring<br>
fields defined on edges.<br>
A README is included with more deail<br>
</p><hr noshade><pre><font color="gray">Added: branches/tools/mpas_draw/Makefile
===================================================================
--- branches/tools/mpas_draw/Makefile (rev 0)
+++ branches/tools/mpas_draw/Makefile 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,29 @@
+CXX=g++
+#CXXFLAGS= -O3
+CXXFLAGS= -O3 -m64
+NETCDF=/opt/netcdf/3.6.3/64
+
+DBGFLAGS= -g -m64
+
+#OFFSETFLAG=-D_64BITOFFSET
+
+PLATFORM=_MACOS
+PLATFORM=_LINUX
+
+ifeq ($(PLATFORM),_LINUX)
+ CXXLIBS = -I$(NETCDF)/include -L$(NETCDF)/lib -lGL -lglut -lnetcdf_c++ -lnetcdf -lGLU -lstdc++
+endif
+
+ifeq ($(PLATFORM),_MACOS)
+ CXXLIBS = -L$(NETCDF)/lib -lnetcdf -lnetcdf_c++ -I$(NETCDF)/include -lstdc++ -framework OpenGL -framework GLUT
+endif
+
+all:
+ $(CXX) mpas_draw.cpp vec_utils.cpp netcdf_utils.cpp $(CXXLIBS) $(CXXFLAGS) $(OFFSETFLAG) -D$(PLATFORM) -o MpasDraw.x
+
+debug:
+ $(CXX) mpas_draw.cpp vec_utils.cpp netcdf_utils.cpp $(CXXLIBS) $(DBGFLAGS) $(OFFSETFLAG) -D$(PLATFORM) -o MpasDraw.x
+
+clean:
+ rm MpasDraw.x
+
Added: branches/tools/mpas_draw/README
===================================================================
--- branches/tools/mpas_draw/README (rev 0)
+++ branches/tools/mpas_draw/README 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,70 @@
+MpasDraw.x:
+ This program is used to visualize MPAS grid, and output files using OpenGL.
+
+ Requirements:
+ You must have opengl and netcdf libraries installed on your computer.
+ Please edit the makefile with the specific paths.
+
+ To Compile:
+ Edit the make file for the appropriate platform, and simply type
+ make
+ The program should build for you, and the executable will be named MpasDraw.x
+
+ To Run:
+ ./MpasDraw.x [file-path]
+
+ [file-path] is an optional argument, that is the path to the mpas grid or output file you are trying to
+ visualize. Optinally, you can input this once the program is running.
+
+ Usage:
+ Once inside the program, you can use the following keys to control the plot.
+
+ Up/Down/Left/Right Arrow Keys:
+ These keys are used to rotate the sphere, to allow for viewing at different angles.
+
+ s d f e , .: (Translation)
+ The lowercase s and f keys control the horizonal position of the image.
+ The lowercase d and e keys control the vertical position of the image.
+ The , and . keys control the depth of the image.
+
+ c: (connectivity)
+ The lowercase c key cycles between cells that are draw, the options are the Delaunay triangulation,
+ the Voronoi diagram, and the edges. Edges are drawn as rectangles made up of two cell centers, and two
+ cell vertices.
+
+ w: (wireframe)
+ The lowercase w key toggles the wireframe on and off.
+
+ b: (colorbar)
+ The lowercase b key toggles the color bar scale.
+ MpasDraw starts with the color bar being scaled on the current slice that is being plotted.
+ Pressing b cycles between this and coloring based on all time and vert levels.
+ Coloring on a single slice is faster, but doesn't show global changes.
+ Coloring on all time takes longer to build the scaling range (if number of time steps is large)
+ but the color bar is static. Also, this colorbar only has to be computed once for the entire run of MpasDraw
+
+ r: (reset)
+ The lowercase r key resets the plot to the default parameters (for color, timestep, vertical level, position, rotation, and color bar).
+
+ v: (variables)
+ The lowercase v key lists available variables that can be chosen from
+ to color the current mesh with. After pressing v, the variables are listed in the terminal.
+ The variable id (given as a number in the left most column) can be entered to
+ change the color of the current plot.
+
+ l, L: (level)
+ The lowercase l key cycles through available vertical levels.
+ The uppercase L key cycles through available vertical levels in reverse.
+
+ t, T: (time)
+ The lowercase t key cycles through available time steps.
+ The uppercase T key cycles through available time steps in reverse.
+
+ [0-9]: (Color Range)
+ Change the percent of available color space.
+ 1-9: 10*Number Pressed% of Color range (eg. 10% = 1, 20% = 2...)
+ 0: 100% of Color range
+
+ ESC:
+ The escape key closes the program.
+
Added: branches/tools/mpas_draw/constants.h
===================================================================
--- branches/tools/mpas_draw/constants.h (rev 0)
+++ branches/tools/mpas_draw/constants.h 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,111 @@
+/*
+ * constants.h
+ * vordraw
+ *
+ * Created by Geoffrey Womeldorff on 3/30/09.
+ * Copyright 2009 FSU. All rights reserved.
+ *
+ */
+
+#define kWindowWidth 400
+#define kWindowHeight 400
+#define KEY_ESCAPE 27
+
+#define KEY_A 65
+#define KEY_B 66
+#define KEY_C 67
+#define KEY_D 68
+#define KEY_E 69
+#define KEY_F 70
+#define KEY_G 71
+#define KEY_H 72
+#define KEY_I 73
+#define KEY_J 74
+#define KEY_K 75
+#define KEY_L 76
+#define KEY_M 77
+#define KEY_N 78
+#define KEY_O 79
+#define KEY_P 80
+#define KEY_Q 81
+#define KEY_R 82
+#define KEY_S 83
+#define KEY_T 84
+#define KEY_U 85
+#define KEY_V 86
+#define KEY_W 87
+#define KEY_X 88
+#define KEY_Y 89
+#define KEY_Z 90
+
+#define KEY_a 97
+#define KEY_b 98
+#define KEY_c 99
+#define KEY_d 100
+#define KEY_e 101
+#define KEY_f 102
+#define KEY_g 103
+#define KEY_h 104
+#define KEY_i 105
+#define KEY_j 106
+#define KEY_k 107
+#define KEY_l 108
+#define KEY_m 109
+#define KEY_n 110
+#define KEY_o 111
+#define KEY_p 112
+#define KEY_q 113
+#define KEY_r 114
+#define KEY_s 115
+#define KEY_t 116
+#define KEY_u 117
+#define KEY_v 118
+#define KEY_w 119
+#define KEY_x 120
+#define KEY_y 121
+#define KEY_z 122
+
+#define KEY_0 48
+#define KEY_1 49
+#define KEY_2 50
+#define KEY_3 51
+#define KEY_4 52
+#define KEY_5 53
+#define KEY_6 54
+#define KEY_7 55
+#define KEY_8 56
+#define KEY_9 57
+
+#define KEY_COMMA 44
+#define KEY_PERIOD 46
+
+#define COLORS 16
+#define PENTA_R 0.1
+#define PENTA_G 0.1
+#define PENTA_B 0.5
+#define HEXA_R 0.0
+#define HEXA_G 1.0
+#define HEXA_B 0.0
+#define HEPTA_R 0.9
+#define HEPTA_G 0.0
+#define HEPTA_B 0.1
+#define OTHER_R 1.0
+#define OTHER_G 0.32549019608
+#define OTHER_B 0.2
+#define LIL_R 0.25882352941
+#define LIL_G 0.75294117647
+#define LIL_B 0.98431372549
+
+#define LINE_R 0.0
+#define LINE_G 0.0
+#define LINE_B 0.0
+
+#define REGL_R 0.0
+#define REGL_G 0.0
+#define REGL_B 0.0
+
+#define REG_R 1.0
+#define REG_G 0.0
+#define REG_B 0.0
+
+#define PI_DP 3.141592653589793
Added: branches/tools/mpas_draw/mpas_draw.cpp
===================================================================
--- branches/tools/mpas_draw/mpas_draw.cpp (rev 0)
+++ branches/tools/mpas_draw/mpas_draw.cpp 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,2103 @@
+#include <cstdlib>
+#include <ctime>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <vector>
+
+#include "constants.h"
+#include "vec_utils.h"
+#include "netcdf_utils.h"
+
+using namespace std;
+
+#ifdef _MACOS
+ #include <GLUT/glut.h>
+#elif _LINUX
+ #include <GL/glut.h>
+#endif
+
+int main ( int argc, char *argv[] );
+void display ( );
+void mouse ( int btn, int state, int x, int y );
+void gl_init ( );
+void myReshape ( int w, int h );
+void spin_image ( );
+void timestamp ( );
+
+void build_connectivity();
+void setup_ranges();
+void build_range(int id);
+
+void rescale_cells();
+void rescale_vertices();
+
+void draw_cells();
+void draw_triangles();
+void draw_edges();
+
+void draw_cell_lines();
+void draw_triangle_lines();
+void draw_edge_lines();
+
+void build_regions();
+void draw_regions();
+
+void color_mesh();
+void color_cells();
+void color_triangles();
+void color_edges();
+
+void arrowKeys( int a_keys, int x, int y );
+void keyPressed( unsigned char key, int x, int y );
+void translateView ( double updown, double leftright);
+void polarView( double distance, double twist, double elevation, double azimuth );
+void hsv_to_rgb(float h, float s, float v, float& r, float& g, float& b);
+
+double getLat(double x, double y, double z);
+double getLon(double x, double y, double z);
+double getX(double lat, double lon);
+double getY(double lat, double lon);
+double getZ(double lat, double lon);
+
+//
+// Global data.
+//
+string filename;
+static GLint axis = 1;
+GLint window;
+int drawing = 0;
+int draw_lines = 1;
+bool on_sphere;
+int color_bar = 0;
+
+double line_factor = 1.002;
+double region_line_factor = 1.006;
+double region_center_factor = 1.008;
+double range_factor = 0.80;
+
+double projUpDown = 0.0;
+double projLeftRight = 0.0;
+double projDistance = 3.0;
+double projTwist = 0.0;
+double projElevation = 0.0;
+double projAzimuth = 0.0;
+
+bool spinning = false;
+static GLfloat theta[3] = { 0.0, 0.0, 0.0 };
+double theta_speed = 0.020;
+
+int ntime;
+int nvertlevels;
+int ncells;
+int nvertices;
+int nedges;
+int maxedges;
+int pixel_height;
+int pixel_width;
+
+int edge_field, cell_field, vertex_field;
+int cur_level = 0;
+int cur_time = 0;
+
+double *xcell;
+double *ycell;
+double *zcell;
+double *xvertex;
+double *yvertex;
+double *zvertex;
+
+int *verticesoncell;
+int *verticesonedge;
+int *cellsonvertex;
+int *cellsonedge;
+int *nedgesoncell;
+
+double *cell_values;
+double *triangle_values;
+double *edge_values;
+
+vector<GLfloat> cell_cells;
+vector<GLfloat> vertex_cells;
+vector<GLfloat> edge_cells;
+vector<GLfloat> cell_lines;
+vector<GLfloat> vertex_lines;
+vector<GLfloat> edge_lines;
+vector<GLfloat> cell_colors;
+vector<GLfloat> vertex_colors;
+vector<GLfloat> edge_colors;
+
+vector< vector<double> > ranges; // 0 - min, 1 - max
+
+vector<GLfloat> region_centers;
+vector<GLfloat> region_lines;
+vector<double> region_radii;
+vector<GLfloat> region_angles;
+bool regions_built = false;
+bool region_draw = false;
+int region_line_div = 180;
+
+double xyz_center[3];
+double xyz_max[3];
+double xyz_min[3];
+double xyz_range[3];
+double xyz_scale;
+
+int main ( int argc, char *argv[] ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // MAIN is the main program for TRIANGULATION_FACES.
+ //
+ // Discussion:
+ //
+ // This program reads certain information from an MPAS NETCDF grid file,
+ // and displays the faces of the triangulation.
+ //
+ // Usage:
+ //
+ // triangulation_faces file.nc
+ //
+ // where
+ //
+ // * file.nc is an MPAS NETCDF grid file.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt, Geoff Womeldorff, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Edward Angel,
+ // Interactive Computer Graphics:
+ // A Top-Down Approach with OpenGL,
+ // Second Edition,
+ // Addison Wesley, 2000.
+ //
+ int cell;
+ int edge;
+ int i;
+ int v;
+
+ cout << "</font>
<font color="blue">";
+ timestamp ( );
+
+ cout << "</font>
<font color="blue">";
+ cout << "MPAS_DRAW:</font>
<font color="blue">";
+ cout << " C++ version</font>
<font color="blue">";
+ cout << " Read an MPAS NETCDF grid file</font>
<font color="blue">";
+ cout << " Visualize the mpas grid/output file.</font>
<font color="blue">";
+ cout << "</font>
<font color="blue">";
+ cout << " Compiled on " << __DATE__ << " at " << __TIME__ << ".</font>
<font color="blue">";
+ //
+ // If the input file was not specified, get it now.
+ //
+ if ( argc <= 1 )
+ {
+ cout << "</font>
<font color="blue">";
+ cout << "MPAS_DRAW:</font>
<font color="blue">";
+ cout << " Please enter the MPAS NETCDF grid filename.</font>
<font color="blue">";
+
+ cin >> filename;
+ }
+ else
+ {
+ filename = argv[1];
+ }
+
+ build_connectivity();
+ setup_ranges();
+ color_cells();
+ color_triangles();
+ color_edges();
+
+ //
+ // Hand things over to OpenGL.
+ //
+ glutInit ( &argc, argv );
+ glutInitDisplayMode ( GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH );
+ glutInitWindowSize ( 800, 800 );
+ glutInitWindowPosition ( 0, 0 );
+ window = glutCreateWindow ( filename.c_str ( ) );
+ glutReshapeFunc ( myReshape );
+ glutDisplayFunc ( display );
+ glutIdleFunc ( spin_image );
+ glutMouseFunc ( mouse );
+ glutKeyboardFunc( keyPressed );
+ glutSpecialFunc( arrowKeys );
+
+ //
+ // Enable hidden surface removal.
+ //
+ glEnable ( GL_DEPTH_TEST );
+
+ gl_init ( );
+
+ glutMainLoop ( );
+
+ //
+ // Things that won't actually happen because we never return from glutMainLoop:
+ //
+ //
+ // Terminate.
+ //
+ cout << "</font>
<font color="blue">";
+ cout << "TRIANGULATION_FACES::</font>
<font color="blue">";
+ cout << " Normal end of execution.</font>
<font color="blue">";
+
+ cout << "</font>
<font color="blue">";
+ timestamp ( );
+
+ return 0;
+}/*}}}*/
+void display ( ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // DISPLAY generates the graphics output.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Geoff Womeldorff, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Edward Angel,
+ // Interactive Computer Graphics:
+ // A Top-Down Approach with OpenGL,
+ // Second Edition,
+ // Addison Wesley, 2000.
+
+ glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
+
+ glMatrixMode( GL_MODELVIEW );
+ glLoadIdentity(); // moves to center of screen
+ translateView ( projUpDown, projLeftRight );
+ polarView( projDistance, projTwist, projElevation, projAzimuth );
+
+ switch(drawing){
+ case 0:
+ draw_triangles();
+ if(draw_lines)
+ draw_triangle_lines();
+ break;
+ case 1:
+ draw_cells();
+ if(draw_lines)
+ draw_cell_lines();
+ break;
+ case 2:
+ draw_edges();
+ if(draw_lines)
+ draw_edge_lines();
+ break;
+ }
+
+ switch(region_draw){
+ case true:
+ draw_regions();
+ break;
+ default:
+ break;
+ }
+
+ //
+ // Clear all the buffers.
+ //
+ glFlush ( );
+ //
+ // Switch between the two buffers for fast animation.
+ //
+ glutSwapBuffers ( );
+
+
+ return;
+}/*}}}*/
+void mouse ( int btn, int state, int x, int y ){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // MOUSE determines the response to mouse input.
+ //
+ // Discussion:
+ //
+ // The original routine assumed the user had a three button mouse, and
+ // dedicated one axis to each.
+ //
+ // Since Apple prefers the esthetics of a one button mouse, we're forced
+ // to live with that choice. This routine alternately pauses rotation,
+ // or increments the rotation axis by 1, no matter which button is pushed.
+ //
+ // Modified:
+ //
+ // 30 December 2008
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Reference:
+ //
+ // Edward Angel,
+ // Interactive Computer Graphics:
+ // A Top-Down Approach with OpenGL,
+ // Second Edition,
+ // Addison Wesley, 2000.
+ //
+ if ( ( btn == GLUT_LEFT_BUTTON && state == GLUT_DOWN ) ||
+ ( btn == GLUT_MIDDLE_BUTTON && state == GLUT_DOWN ) ||
+ ( btn == GLUT_RIGHT_BUTTON && state == GLUT_DOWN ) ) {
+ if ( spinning ) {
+ spinning = false;
+ theta_speed = 0.0;
+ } else {
+ spinning = true;
+ theta_speed = 0.020;
+ axis = axis + 1;
+ }
+ }
+
+ axis = axis % 3;
+
+ return;
+}/*}}}*/
+void gl_init ( ){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // gl_init initializes OpenGL state variables dealing with viewing and attributes.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Geoff Womeldorff, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Edward Angel,
+ // Interactive Computer Graphics:
+ // A Top-Down Approach with OpenGL,
+ // Second Edition,
+ // Addison Wesley, 2000.
+ GLfloat line_width;
+ GLfloat point_size;
+ //
+ // Set the background to WHITE.
+ //
+ glClearColor ( 1.0, 1.0, 1.0, 1.0 );
+ //
+ // Antialiasing.
+ //
+ glDepthFunc( GL_LESS );
+ glEnable(GL_DEPTH_TEST);
+ glBlendFunc ( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA );
+ glHint ( GL_LINE_SMOOTH_HINT, GL_DONT_CARE );
+ glShadeModel( GL_SMOOTH );
+
+ glEnableClientState( GL_VERTEX_ARRAY );
+ glEnableClientState( GL_COLOR_ARRAY );
+
+
+ glEnable ( GL_POINT_SMOOTH );
+ glEnable ( GL_LINE_SMOOTH );
+ //
+ // The default point size is 1.0.
+ //
+ if ( ncells <= 400 ) {
+ point_size = 16.0;
+ } else if ( ncells <= 800 ) {
+ point_size = 8.0;
+ } else if ( ncells <= 1600 ) {
+ point_size = 4.0;
+ } else if ( ncells <= 3200 ) {
+ point_size = 2.0;
+ } else {
+ point_size = 1.0;
+ }
+
+ point_size = 8.0;
+ glPointSize ( point_size );
+ //
+ // The default line width is 1.0.
+ //
+ if ( ncells <= 1600 ){
+ line_width = 1.0;
+ } else {
+ line_width = 0.1;
+ }
+ glLineWidth ( line_width );
+
+ return;
+}/*}}}*/
+void myReshape ( int w, int h ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // MYRESHAPE determines the window mapping.
+ //
+ // Modified:
+ //
+ // 30 December 2008
+ //
+ // Author:
+ //
+ // John Burkardt, Geoff Womeldorff, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Edward Angel,
+ // Interactive Computer Graphics:
+ // A Top-Down Approach with OpenGL,
+ // Second Edition,
+ // Addison Wesley, 2000.
+ //
+ //
+
+ if ( w <= h )
+ {
+ glOrtho (
+ -1.05, +1.05,
+ - 1.05 * ( GLfloat ) h / ( GLfloat ) w, + 1.05 * ( GLfloat ) h / ( GLfloat ) w,
+ -10.0, 10.0 );
+ }
+ else
+ {
+ glOrtho (
+ - 1.05 * ( GLfloat ) h / ( GLfloat ) w, + 1.05 * ( GLfloat ) h / ( GLfloat ) w,
+ - 1.05, + 1.05,
+ -10.0, 10.0 );
+ }
+
+ glViewport ( 0, 0, w, h );
+ glMatrixMode ( GL_PROJECTION );
+ glLoadIdentity ( );
+
+ gluPerspective( 45.0f, (GLfloat)w / (GLfloat)h, 0.1f, 5.0f );
+ //glMatrixMode ( GL_MODELVIEW );
+
+ return;
+}/*}}}*/
+void spin_image ( ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // SPIN_IMAGE adjusts the angle of rotation and redisplays the picture.
+ //
+ // Modified:
+ //
+ // 15 December 2008
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Reference:
+ //
+ // Edward Angel,
+ // Interactive Computer Graphics:
+ // A Top-Down Approach with OpenGL,
+ // Second Edition,
+ // Addison Wesley, 2000.
+ //
+ theta[axis] = theta[axis] + theta_speed;
+
+ if ( 360.0 < theta[axis] )
+ {
+ theta[axis] = theta[axis] - 360.0;
+ }
+ glutPostRedisplay ( );
+
+ return;
+}/*}}}*/
+void timestamp ( ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // TIMESTAMP prints the current YMDHMS date as a time stamp.
+ //
+ // Example:
+ //
+ // 31 May 2001 09:45:54 AM
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 08 July 2009
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Parameters:
+ //
+ // None
+ //
+# define TIME_SIZE 40
+
+ static char time_buffer[TIME_SIZE];
+ const struct std::tm *tm_ptr;
+ size_t len;
+ std::time_t now;
+
+ now = std::time ( NULL );
+ tm_ptr = std::localtime ( &now );
+
+ len = std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );
+
+ std::cout << time_buffer << "</font>
<font color="blue">";
+
+ return;
+# undef TIME_SIZE
+}/*}}}*/
+
+void build_connectivity(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // BUILD_CONNECTIVITY builds the connectivity arrays for Display to draw
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Geoff Womeldorff, Doug Jacobsen
+ //
+
+ int i, j;
+ int v1, v2;
+ int c1, c2, c3, dc;
+ int n_lilgons; // Smaller than 5 sides
+ int n_pentagons; // 5 Sides
+ int n_heptagons; // 7 Sides
+ int n_bigagons; // Larger than 7 sides
+
+
+ //
+ // Get sizes.
+ //
+ ntime = netcdf_mpas_read_dim(filename, "Time" );
+ nvertlevels = netcdf_mpas_read_dim(filename, "nVertLevels");
+ ncells = netcdf_mpas_read_dim (filename, "nCells");
+ nvertices = netcdf_mpas_read_dim(filename, "nVertices");
+ nedges = netcdf_mpas_read_dim(filename, "nEdges");
+ maxedges = netcdf_mpas_read_dim(filename, "maxEdges");
+ on_sphere = netcdf_mpas_read_onsphere(filename);
+
+ cell_cells.reserve((ncells+1)*maxedges*3*3);
+ cell_lines.reserve((ncells+1)*maxedges*3*2);
+ vertex_cells.reserve((nvertices+1)*3*3);
+ vertex_lines.reserve((nvertices+1)*3*3*2);
+ edge_cells.reserve((nedges+1)*4*3);
+ edge_lines.reserve((nedges+1)*4*2*3);
+
+ cout << "</font>
<font color="blue">";
+ if(!on_sphere){
+ line_factor = 0.002;
+ cout << " Points are NOT on a sphere. " << endl;
+ } else {
+ cout << " Points ARE on a sphere. " << endl;
+ }
+
+ cout << " The number of time steps NTIME = " << ntime << "</font>
<font color="blue">";
+ cout << " The number of vertical levels NVERTLEVELS = " << nvertlevels << "</font>
<font color="blue">";
+ cout << " The number of cells NCELLS = " << ncells << "</font>
<font color="blue">";
+ cout << " The number of vertices NVERTICES = " << nvertices << "</font>
<font color="blue">";
+ cout << " The number of edges NEDGES = " << nedges << "</font>
<font color="blue">";
+
+ xcell = new double[ncells];
+ ycell = new double[ncells];
+ zcell = new double[ncells];
+
+ xvertex = new double[nvertices];
+ yvertex = new double[nvertices];
+ zvertex = new double[nvertices];
+
+ verticesoncell = new int[maxedges*ncells];
+ verticesonedge = new int[2*nedges];
+ cellsonvertex = new int[3*nvertices];
+ cellsonedge = new int[2*nedges];
+ nedgesoncell = new int[ncells];
+
+ netcdf_mpas_read_xyzcell ( filename, ncells, xcell, ycell, zcell );
+ netcdf_mpas_read_xyzvertex ( filename, nvertices, xvertex, yvertex, zvertex );
+
+ rescale_cells();
+ rescale_vertices();
+
+ netcdf_mpas_read_verticesoncell ( filename, maxedges, ncells, verticesoncell );
+ netcdf_mpas_read_verticesonedge ( filename, nedges, verticesonedge );
+ netcdf_mpas_read_cellsonvertex ( filename, nvertices, cellsonvertex );
+ netcdf_mpas_read_cellsonedge ( filename, nedges, cellsonedge );
+ netcdf_mpas_read_nedgesoncell ( filename, ncells, nedgesoncell );
+
+ n_lilgons = 0;
+ n_pentagons = 0;
+ n_heptagons = 0;
+ n_bigagons = 0;
+
+ //
+ // connectivity is 1 based. Fix that.
+ //
+ for ( i = 0; i < nvertices; i++ )
+ {
+ for ( j = 0; j < 3; j++ )
+ {
+ cellsonvertex[i*3+j]--;
+ }
+ }
+
+ for( i = 0; i < ncells; i++){
+ for ( j = 0; j < maxedges; j++){
+ verticesoncell[i*maxedges + j]--;
+ }
+ if(nedgesoncell[i] < 5){
+ n_lilgons++;
+ } else if(nedgesoncell[i] == 5){
+ n_pentagons++;
+ } else if(nedgesoncell[i] == 7){
+ n_heptagons++;
+ } else if(nedgesoncell[i] > 7){
+ n_bigagons++;
+ }
+ }
+
+ cout << " Number of Lilgons (< 5 sides) " << n_lilgons << endl;
+ cout << " Number of Pentagons " << n_pentagons << endl;
+ cout << " Number of Heptagons " << n_heptagons << endl;
+ cout << " Number of Bigagons (> 7 sides) " << n_bigagons << endl;
+
+ for( i = 0; i < nedges; i++){
+ for( j = 0; j < 2; j++){
+ cellsonedge[i*2 + j]--;
+ verticesonedge[i*2 + j]--;
+ }
+ }
+
+
+ for(i = 0; i < nvertices; i++){
+ c1 = cellsonvertex[i*3];
+ c2 = cellsonvertex[i*3 + 1];
+ c3 = cellsonvertex[i*3 + 2];
+
+ if(c1 == -1){
+ if(c2 == -1){
+ dc = c3;
+ } else {
+ dc = c2;
+ }
+ } else {
+ dc = c1;
+ }
+
+ if(c1 == -1){
+ c1 = dc;
+ }
+ if(c2 == -1){
+ c2 = dc;
+ }
+ if(c3 == -1){
+ c3 = dc;
+ }
+ vertex_cells.push_back(xcell[c1]);
+ vertex_cells.push_back(ycell[c1]);
+ vertex_cells.push_back(zcell[c1]);
+
+ vertex_cells.push_back(xcell[c2]);
+ vertex_cells.push_back(ycell[c2]);
+ vertex_cells.push_back(zcell[c2]);
+
+ vertex_cells.push_back(xcell[c3]);
+ vertex_cells.push_back(ycell[c3]);
+ vertex_cells.push_back(zcell[c3]);
+
+ if(on_sphere){
+ vertex_lines.push_back(xcell[c1]*line_factor);
+ vertex_lines.push_back(ycell[c1]*line_factor);
+ vertex_lines.push_back(zcell[c1]*line_factor);
+ vertex_lines.push_back(xcell[c2]*line_factor);
+ vertex_lines.push_back(ycell[c2]*line_factor);
+ vertex_lines.push_back(zcell[c2]*line_factor);
+
+ vertex_lines.push_back(xcell[c2]*line_factor);
+ vertex_lines.push_back(ycell[c2]*line_factor);
+ vertex_lines.push_back(zcell[c2]*line_factor);
+ vertex_lines.push_back(xcell[c3]*line_factor);
+ vertex_lines.push_back(ycell[c3]*line_factor);
+ vertex_lines.push_back(zcell[c3]*line_factor);
+
+ vertex_lines.push_back(xcell[c3]*line_factor);
+ vertex_lines.push_back(ycell[c3]*line_factor);
+ vertex_lines.push_back(zcell[c3]*line_factor);
+ vertex_lines.push_back(xcell[c1]*line_factor);
+ vertex_lines.push_back(ycell[c1]*line_factor);
+ vertex_lines.push_back(zcell[c1]*line_factor);
+ } else {
+ vertex_lines.push_back(xcell[c1]);
+ vertex_lines.push_back(ycell[c1]);
+ vertex_lines.push_back(line_factor);
+ vertex_lines.push_back(xcell[c2]);
+ vertex_lines.push_back(ycell[c2]);
+ vertex_lines.push_back(line_factor);
+
+ vertex_lines.push_back(xcell[c2]);
+ vertex_lines.push_back(ycell[c2]);
+ vertex_lines.push_back(line_factor);
+ vertex_lines.push_back(xcell[c3]);
+ vertex_lines.push_back(ycell[c3]);
+ vertex_lines.push_back(line_factor);
+
+ vertex_lines.push_back(xcell[c3]);
+ vertex_lines.push_back(ycell[c3]);
+ vertex_lines.push_back(line_factor);
+ vertex_lines.push_back(xcell[c1]);
+ vertex_lines.push_back(ycell[c1]);
+ vertex_lines.push_back(line_factor);
+ }
+ }
+
+ delete(cellsonvertex);
+
+ for(i = 0; i < ncells; i++){
+ for(j = 0; j < nedgesoncell[i]; j++){
+ v1 = verticesoncell[i*maxedges + j%nedgesoncell[i]];
+ v2 = verticesoncell[i*maxedges + (j+1)%nedgesoncell[i]];
+
+ cell_cells.push_back(xcell[i]);
+ cell_cells.push_back(ycell[i]);
+ cell_cells.push_back(zcell[i]);
+ cell_cells.push_back(xvertex[v1]);
+ cell_cells.push_back(yvertex[v1]);
+ cell_cells.push_back(zvertex[v1]);
+ cell_cells.push_back(xvertex[v2]);
+ cell_cells.push_back(yvertex[v2]);
+ cell_cells.push_back(zvertex[v2]);
+
+ if(on_sphere){
+ cell_lines.push_back(xvertex[v1]*line_factor);
+ cell_lines.push_back(yvertex[v1]*line_factor);
+ cell_lines.push_back(zvertex[v1]*line_factor);
+
+ cell_lines.push_back(xvertex[v2]*line_factor);
+ cell_lines.push_back(yvertex[v2]*line_factor);
+ cell_lines.push_back(zvertex[v2]*line_factor);
+ } else {
+ cell_lines.push_back(xvertex[v1]);
+ cell_lines.push_back(yvertex[v1]);
+ cell_lines.push_back(line_factor);
+
+ cell_lines.push_back(xvertex[v2]);
+ cell_lines.push_back(yvertex[v2]);
+ cell_lines.push_back(line_factor);
+ }
+ }
+ }
+
+ delete(verticesoncell);
+
+ for(i = 0; i < nedges; i++){
+ c1 = cellsonedge[i*2];
+ c2 = cellsonedge[i*2+1];
+ v1 = verticesonedge[i*2];
+ v2 = verticesonedge[i*2+1];
+
+ if(c1 == -1){
+ c1 = c2;
+ } else if(c2 == -1){
+ c2 = c1;
+ }
+
+ if(v1 == -1){
+ v1 = v2;
+ } else if (v2 == -1){
+ v2 = v1;
+ }
+ edge_cells.push_back(xcell[c1]);
+ edge_cells.push_back(ycell[c1]);
+ edge_cells.push_back(zcell[c1]);
+
+ edge_cells.push_back(xvertex[v1]);
+ edge_cells.push_back(yvertex[v1]);
+ edge_cells.push_back(zvertex[v1]);
+
+ edge_cells.push_back(xvertex[v2]);
+ edge_cells.push_back(yvertex[v2]);
+ edge_cells.push_back(zvertex[v2]);
+
+ edge_cells.push_back(xcell[c2]);
+ edge_cells.push_back(ycell[c2]);
+ edge_cells.push_back(zcell[c2]);
+
+ edge_cells.push_back(xvertex[v2]);
+ edge_cells.push_back(yvertex[v2]);
+ edge_cells.push_back(zvertex[v2]);
+
+ edge_cells.push_back(xvertex[v1]);
+ edge_cells.push_back(yvertex[v1]);
+ edge_cells.push_back(zvertex[v1]);
+
+ if(on_sphere){
+ edge_lines.push_back(xcell[c1]*line_factor);
+ edge_lines.push_back(ycell[c1]*line_factor);
+ edge_lines.push_back(zcell[c1]*line_factor);
+ edge_lines.push_back(xvertex[v1]*line_factor);
+ edge_lines.push_back(yvertex[v1]*line_factor);
+ edge_lines.push_back(zvertex[v1]*line_factor);
+
+ edge_lines.push_back(xvertex[v1]*line_factor);
+ edge_lines.push_back(yvertex[v1]*line_factor);
+ edge_lines.push_back(zvertex[v1]*line_factor);
+ if(c1 == c2){
+ edge_lines.push_back(xvertex[v2]*line_factor);
+ edge_lines.push_back(yvertex[v2]*line_factor);
+ edge_lines.push_back(zvertex[v2]*line_factor);
+ } else {
+ edge_lines.push_back(xcell[c2]*line_factor);
+ edge_lines.push_back(ycell[c2]*line_factor);
+ edge_lines.push_back(zcell[c2]*line_factor);
+
+ edge_lines.push_back(xcell[c2]*line_factor);
+ edge_lines.push_back(ycell[c2]*line_factor);
+ edge_lines.push_back(zcell[c2]*line_factor);
+ edge_lines.push_back(xvertex[v2]*line_factor);
+ edge_lines.push_back(yvertex[v2]*line_factor);
+ edge_lines.push_back(zvertex[v2]*line_factor);
+ }
+
+ edge_lines.push_back(xvertex[v2]*line_factor);
+ edge_lines.push_back(yvertex[v2]*line_factor);
+ edge_lines.push_back(zvertex[v2]*line_factor);
+ edge_lines.push_back(xcell[c1]*line_factor);
+ edge_lines.push_back(ycell[c1]*line_factor);
+ edge_lines.push_back(zcell[c1]*line_factor);
+ } else {
+ edge_lines.push_back(xcell[c1]);
+ edge_lines.push_back(ycell[c1]);
+ edge_lines.push_back(line_factor);
+ edge_lines.push_back(xvertex[v1]);
+ edge_lines.push_back(yvertex[v1]);
+ edge_lines.push_back(line_factor);
+
+ edge_lines.push_back(xvertex[v1]);
+ edge_lines.push_back(yvertex[v1]);
+ edge_lines.push_back(line_factor);
+ if(c1 == c2){
+ edge_lines.push_back(xvertex[v2]);
+ edge_lines.push_back(yvertex[v2]);
+ edge_lines.push_back(line_factor);
+ } else {
+ edge_lines.push_back(xcell[c2]);
+ edge_lines.push_back(ycell[c2]);
+ edge_lines.push_back(line_factor);
+
+ edge_lines.push_back(xcell[c2]);
+ edge_lines.push_back(ycell[c2]);
+ edge_lines.push_back(line_factor);
+ edge_lines.push_back(xvertex[v2]);
+ edge_lines.push_back(yvertex[v2]);
+ edge_lines.push_back(line_factor);
+ }
+
+ edge_lines.push_back(xvertex[v2]);
+ edge_lines.push_back(yvertex[v2]);
+ edge_lines.push_back(line_factor);
+ edge_lines.push_back(xcell[c1]);
+ edge_lines.push_back(ycell[c1]);
+ edge_lines.push_back(line_factor);
+ }
+ }
+
+ delete(verticesonedge);
+ delete(cellsonedge);
+
+ cell_values = new double[ncells];
+ triangle_values = new double[nvertices];
+ edge_values = new double[nedges];
+
+ color_cells();
+ color_triangles();
+ color_edges();
+
+ delete [] xcell;
+ delete [] ycell;
+ delete [] zcell;
+ delete [] xvertex;
+ delete [] yvertex;
+ delete [] zvertex;
+// delete [] verticesoncell;
+// delete [] verticesonedge;
+// delete [] cellsonvertex;
+// delete [] cellsonedge;
+}/*}}}*/
+void setup_ranges(){/*{{{*/
+ int num_vars;
+
+ num_vars = netcdf_mpas_read_num_vars(filename);
+
+ ranges.resize(num_vars);
+
+ return;
+}/*}}}*/
+void build_range(int id){/*{{{*/
+ int num_items;
+ double max, min;
+ vector<double> temp_data;
+
+ if(color_bar = 1){
+ if(ranges[id].size() == 0){
+ cout << "Building min-max range for coloring." << endl;
+ num_items = netcdf_mpas_field_num_items(filename, id);
+
+ temp_data.resize(num_items);
+
+ netcdf_mpas_read_full_field(filename, id, &temp_data[0]);
+
+ r8vec_min_max(num_items, &temp_data[0], min, max);
+
+ ranges[id].push_back(min);
+ ranges[id].push_back(max);
+ cout << "Range build complete." << endl;
+ }
+ } else {
+ switch(drawing){
+ case 0:
+ num_items = nvertices;
+ break;
+ case 1:
+ num_items = ncells;
+ break;
+ case 2:
+ num_items = nedges;
+ break;
+ default:
+ return;
+ break;
+ }
+
+ temp_data.resize(num_items);
+
+ netcdf_mpas_read_field(filename, id, &temp_data.at(0), cur_time, cur_level);
+ }
+}/*}}}*/
+
+void rescale_cells(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // RESCALE_CELLS scales xcell, ycell, and zcell arrays to have a norm of 1
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+
+
+ int i;
+ int cell;
+ double norm;
+ xyz_min[0] = r8vec_min ( ncells, xcell );
+ xyz_max[0] = r8vec_max ( ncells, xcell );
+ xyz_min[1] = r8vec_min ( ncells, ycell );
+ xyz_max[1] = r8vec_max ( ncells, ycell );
+ xyz_min[2] = r8vec_min ( ncells, zcell );
+ xyz_max[2] = r8vec_max ( ncells, zcell );
+
+ xyz_range[0] = xyz_max[0] - xyz_min[0];
+ xyz_range[1] = xyz_max[1] - xyz_min[1];
+ xyz_range[2] = xyz_max[2] - xyz_min[2];
+
+ if ( xyz_range[0] == 0.0 ){
+ cout << "</font>
<font color="blue">";
+ cout << "rescale_cells(): - Fatal error!</font>
<font color="blue">";
+ cout << " The X data range is 0.</font>
<font color="blue">";
+ exit ( 1 );
+ }
+
+ if ( xyz_range[1] == 0.0 ){
+ cout << "</font>
<font color="blue">";
+ cout << "rescale_cells(): - Fatal error!</font>
<font color="blue">";
+ cout << " The Y data range is 0.</font>
<font color="blue">";
+ exit ( 1 );
+ }
+
+ if(on_sphere){
+ if ( xyz_range[2] == 0.0 ){
+ cout << "</font>
<font color="blue">";
+ cout << "rescale_cells(): - Fatal error!</font>
<font color="blue">";
+ cout << " The Z data range is 0.</font>
<font color="blue">";
+ exit ( 1 );
+ }
+ }
+
+ xyz_scale = 0.0;
+ for (i = 0; i < 3; i++ )
+ {
+ xyz_center[i] = ( xyz_min[i] + xyz_max[i] ) / 2.0;
+ xyz_scale = std::max ( xyz_scale, ( xyz_max[i] - xyz_min[i] ) / 2.0 );
+ }
+ //
+ // A sphere doesn't need this rescaling.
+ // A box does!
+ //
+ //xyz_scale = sqrt ( 3.0 ) * xyz_scale;
+ //
+ // Translate the data so it is centered.
+ // Scale the data so it fits in the unit cube.
+ //
+ for ( cell = 0; cell < ncells; cell++ )
+ {
+ xcell[cell] = ( xcell[cell] - xyz_center[0] ) / xyz_scale;
+ ycell[cell] = ( ycell[cell] - xyz_center[1] ) / xyz_scale;
+ zcell[cell] = ( zcell[cell] - xyz_center[2] ) / xyz_scale;
+
+ norm = sqrt(xcell[cell]*xcell[cell] + ycell[cell]*ycell[cell] + zcell[cell]*zcell[cell]);
+ xcell[cell] /= norm;
+ ycell[cell] /= norm;
+ zcell[cell] /= norm;
+ }
+
+}/*}}}*/
+void rescale_vertices(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // RESCALE_VERTICES scales xvertex, yvertex, and zvertex arrays to have a norm of 1
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+
+ int i;
+ int vertex;
+ double norm;
+ xyz_min[0] = r8vec_min ( nvertices, xvertex );
+ xyz_max[0] = r8vec_max ( nvertices, xvertex );
+ xyz_min[1] = r8vec_min ( nvertices, yvertex );
+ xyz_max[1] = r8vec_max ( nvertices, yvertex );
+ xyz_min[2] = r8vec_min ( nvertices, zvertex );
+ xyz_max[2] = r8vec_max ( nvertices, zvertex );
+
+ xyz_range[0] = xyz_max[0] - xyz_min[0];
+ xyz_range[1] = xyz_max[1] - xyz_min[1];
+ xyz_range[2] = xyz_max[2] - xyz_min[2];
+
+ if ( xyz_range[0] == 0.0 ){
+ cout << "</font>
<font color="blue">";
+ cout << "rescale_vertices(): - Fatal error!</font>
<font color="blue">";
+ cout << " The X data range is 0.</font>
<font color="blue">";
+ exit ( 1 );
+ }
+
+ if ( xyz_range[1] == 0.0 ){
+ cout << "</font>
<font color="blue">";
+ cout << "rescale_vertices(): - Fatal error!</font>
<font color="blue">";
+ cout << " The Y data range is 0.</font>
<font color="blue">";
+ exit ( 1 );
+ }
+ if(on_sphere){
+ if ( xyz_range[2] == 0.0 ){
+ cout << "</font>
<font color="blue">";
+ cout << "rescale_vertices(): - Fatal error!</font>
<font color="blue">";
+ cout << " The Z data range is 0.</font>
<font color="gray">";
+ exit ( 1 );
+ }
+ }
+
+ xyz_scale = 0.0;
+ for (i = 0; i < 3; i++ )
+ {
+ xyz_center[i] = ( xyz_min[i] + xyz_max[i] ) / 2.0;
+ xyz_scale = std::max ( xyz_scale, ( xyz_max[i] - xyz_min[i] ) / 2.0 );
+ }
+ //
+ // A sphere doesn't need this rescaling.
+ // A box does!
+ //
+ //xyz_scale = sqrt ( 3.0 ) * xyz_scale;
+ //
+ // Translate the data so it is centered.
+ // Scale the data so it fits in the unit cube.
+ //
+ for ( vertex = 0; vertex < nvertices; vertex++ )
+ {
+ xvertex[vertex] = ( xvertex[vertex] - xyz_center[0] ) / xyz_scale;
+ yvertex[vertex] = ( yvertex[vertex] - xyz_center[1] ) / xyz_scale;
+ zvertex[vertex] = ( zvertex[vertex] - xyz_center[2] ) / xyz_scale;
+
+ norm = sqrt(xvertex[vertex]*xvertex[vertex] + yvertex[vertex]*yvertex[vertex] + zvertex[vertex]*zvertex[vertex]);
+ xvertex[vertex] /= norm;
+ yvertex[vertex] /= norm;
+ zvertex[vertex] /= norm;
+ }
+
+}/*}}}*/
+
+void draw_triangles ( ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // DRAW_TRIANGLES draws the Delaunay triangles
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+ glColorPointer( 3, GL_FLOAT, 0, &vertex_colors[0] );
+ glVertexPointer( 3, GL_FLOAT, 0, &vertex_cells[0] );
+ glDrawArrays( GL_TRIANGLES, 0, vertex_cells.size()/3);
+
+ return;
+}/*}}}*/
+void draw_cells ( ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // DRAW_CELLS draws the Voronoi cells
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+ glColorPointer( 3, GL_FLOAT, 0, &cell_colors[0] );
+ glVertexPointer( 3, GL_FLOAT, 0, &cell_cells[0] );
+ glDrawArrays( GL_TRIANGLES, 0, cell_cells.size()/3);
+
+ return;
+}/*}}}*/
+void draw_edges ( ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // DRAW_EDGES draws quadrilaterals to represent edge values
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+ glColorPointer( 3, GL_FLOAT, 0, &edge_colors[0] );
+ glVertexPointer( 3, GL_FLOAT, 0, &edge_cells[0] );
+ glDrawArrays( GL_TRIANGLES, 0, edge_cells.size()/3);
+
+ return;
+}/*}}}*/
+
+void draw_triangle_lines(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // DRAW_TRIANGLE_LINES draws the lines for the Delaunay triangle edges
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+ glDisableClientState( GL_COLOR_ARRAY );
+ glColor3f ( LINE_R, LINE_G, LINE_B );
+ glVertexPointer( 3, GL_FLOAT, 0, &vertex_lines[0] );
+ glDrawArrays( GL_LINES, 0, vertex_lines.size()/3 );
+ glEnableClientState( GL_COLOR_ARRAY);
+
+ return;
+}/*}}}*/
+void draw_cell_lines(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // DRAW_CELL_LINES draws the lines for the Voronoi cell edges.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+ glDisableClientState( GL_COLOR_ARRAY );
+ glColor3f ( LINE_R, LINE_G, LINE_B );
+ glVertexPointer( 3, GL_FLOAT, 0, &cell_lines[0] );
+ glDrawArrays( GL_LINES, 0, cell_lines.size()/3 );
+ glEnableClientState( GL_COLOR_ARRAY);
+
+ return;
+}/*}}}*/
+void draw_edge_lines(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // DRAW_EDGE_LINES draws the lines for the edge "edges".
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+ glDisableClientState( GL_COLOR_ARRAY );
+ glColor3f ( LINE_R, LINE_G, LINE_B );
+ glVertexPointer( 3, GL_FLOAT, 0, &edge_lines[0] );
+ glDrawArrays( GL_LINES, 0, edge_lines.size()/3 );
+ glEnableClientState( GL_COLOR_ARRAY);
+
+ return;
+}/*}}}*/
+
+void build_regions(){/*{{{*/
+ ifstream reg_in("RegionCenters");
+ ifstream rad_in("RegionRadii");
+ double t;
+ double x, y, z, r;
+ double scale;
+ double planar_r;
+ int i, j;
+
+ if(!reg_in){
+ cout << "File RegionCenters not found. Skipping regions." << endl;
+ return;
+ } else {
+ while(!reg_in.eof()){
+ reg_in >> x;
+ reg_in >> y;
+ reg_in >> z;
+
+ if(reg_in.good()){
+ region_centers.push_back(x);
+ region_centers.push_back(y);
+ region_centers.push_back(z);
+ }
+ }
+ }
+
+ reg_in.close();
+
+ if(!rad_in){
+ cout << "File RegionRadii not found. Skipping regions." << endl;
+ return;
+ } else {
+ while(!rad_in.eof()){
+ rad_in >> r;
+
+ if(rad_in.good()){
+ region_radii.push_back(r);
+ }
+ }
+ }
+
+ rad_in.close();
+
+ for(i = 0; i < region_radii.size(); i++){
+ x = region_centers.at(i*3);
+ y = region_centers.at(i*3+1);
+ z = region_centers.at(i*3+2);
+
+ scale = sin(region_radii.at(i) + M_PI/2.0);
+
+ region_lines.push_back(x*scale*region_line_factor);
+ region_lines.push_back(y*scale*region_line_factor);
+ region_lines.push_back(z*scale*region_line_factor);
+
+ region_radii.at(i) = sin(region_radii.at(i));
+
+ region_angles.push_back(atan2(region_centers.at(i*3+1),region_centers.at(i*3))*180.0/M_PI);
+ region_angles.push_back((asin(region_centers.at(i*3+2))+M_PI/2.0)*180.0/M_PI);
+ }
+
+
+ for(i = 0; i < region_centers.size(); i++){
+ region_centers.at(i) = region_centers.at(i)*region_center_factor;
+ }
+
+ regions_built = true;
+}/*}}}*/
+void draw_regions(){/*{{{*/
+ double theta, phi;
+ float h, s, v;
+ float r, g, b;
+ glDisableClientState( GL_COLOR_ARRAY );
+
+ s = 0.8;
+ v = 0.0;
+
+ for(int i = 0; i < region_radii.size(); i++){
+ phi = region_angles.at(i*2);
+ theta = region_angles.at(i*2+1);
+ h = (1.0*i/region_radii.size()) * 0.8;
+
+ hsv_to_rgb(h,s,v,r,g,b);
+
+ glPushMatrix();
+ glColor3f( r, g, b);
+ glTranslated(region_lines.at(i*3), region_lines.at(i*3+1), region_lines.at(i*3+2));
+
+ glRotated(phi, 0.0, 0.0, 1.0);
+ glRotated(-theta, 0.0, 1.0, 0.0);
+
+ gluCylinder(gluNewQuadric(), region_radii.at(i)*1.015, region_radii.at(i)*1.015, 0.015, 360, 5);
+ gluDisk(gluNewQuadric(), region_radii.at(i), region_radii.at(i)*1.015, 360, 5);
+ glPopMatrix();
+ }
+
+ glColor3f( REG_R, REG_G, REG_B );
+ glVertexPointer( 3, GL_FLOAT, 0, &region_centers[0] );
+ glDrawArrays( GL_POINTS, 0, region_centers.size()/3.0 );
+
+ glEnableClientState( GL_COLOR_ARRAY );
+}/*}}}*/
+
+void color_mesh(){/*{{{*/
+ switch(drawing){
+ case 0:
+ color_triangles();
+ break;
+ case 1:
+ color_cells();
+ break;
+ case 2:
+ color_edges();
+ break;
+ default:
+ break;
+ }
+ return;
+}/*}}}*/
+void color_cells(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // COLOR_CELLS builds the color array used to display the Voronoi cells
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+
+ int i, j, k, o;
+ double max, min;
+ long *dims;
+ int num_dims = 0;
+ int cell_dim;
+ int vert_dim;
+ int time_dim;
+ int num_items;
+ float h, s, v;
+ float r, g, b;
+
+ cell_colors.clear();
+
+ s = 1.0;
+ v = 1.0;
+ if(cell_field == 0){
+ for(i = 0; i < ncells; i++){
+ o = nedgesoncell[i];
+
+ for(j = 0; j < 3*o; j++){
+ cell_colors.push_back(0.8);
+ cell_colors.push_back(0.8);
+ cell_colors.push_back(0.8);
+ }
+ }
+ } else {
+ netcdf_mpas_read_field(filename, cell_field, cell_values, cur_time, cur_level);
+
+ if(color_bar == 1){
+ min = ranges[cell_field].at(0);
+ max = ranges[cell_field].at(1);
+ } else {
+ r8vec_min_max(ncells, cell_values, min, max);
+ }
+
+ cout << "Min: " << min << " Max: " << max << endl;
+
+ for(i = 0; i < ncells; i++){
+ if((max-min) != 0.0){
+ h = (cell_values[i] - min)/(max-min) * range_factor;
+ } else {
+ h = (cell_values[i] - min)/1.0 * range_factor;
+ }
+
+ o = nedgesoncell[i];
+
+ hsv_to_rgb(h, s, v, r, g, b);
+
+ for(j = 0; j < 3*o; j++){
+ cell_colors.push_back(r);
+ cell_colors.push_back(g);
+ cell_colors.push_back(b);
+ }
+ }
+ }
+
+ return;
+}/*}}}*/
+void color_triangles(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // COLOR_TRIANGLES builds the color array used to display the Delaunay triangles
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+
+ int i, j;
+ double max, min;
+ float h, s, v;
+ float r, g, b;
+
+ vertex_colors.clear();
+
+ s = 1.0;
+ v = 1.0;
+ if(vertex_field == 0){
+ for(i = 0; i < nvertices; i++){
+ for(j = 0; j < 3; j++){
+ vertex_colors.push_back(0.8);
+ vertex_colors.push_back(0.8);
+ vertex_colors.push_back(0.8);
+ }
+ }
+ } else {
+ netcdf_mpas_read_field(filename, vertex_field, triangle_values, cur_time, cur_level);
+
+ if(color_bar == 1){
+ max = ranges[vertex_field][0];
+ min = ranges[vertex_field][1];
+ } else {
+ r8vec_min_max(nvertices, triangle_values, min, max);
+ }
+
+ cout << "Min: " << min << " Max: " << max << endl;
+
+ for(i = 0; i < nvertices; i++){
+ if(max-min != 0.0){
+ h = (triangle_values[i] - min)/(max-min)*range_factor;
+ }else{
+ h = (triangle_values[i] - min)/1.0 * range_factor;
+ }
+
+ hsv_to_rgb(h, s, v, r, g, b);
+
+ for(j = 0; j < 3; j++){
+ vertex_colors.push_back(r);
+ vertex_colors.push_back(g);
+ vertex_colors.push_back(b);
+ }
+ }
+ }
+
+ return;
+}/*}}}*/
+void color_edges(){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // COLOR_EDGES builds the color array used to display the edge quads
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+ int i, j, o;
+ double max, min;
+ float h, s, v;
+ float r, g, b;
+
+ edge_colors.clear();
+
+ s = 1.0;
+ v = 1.0;
+ if(edge_field == 0){
+ for(i = 0; i < nedges; i++){
+ for(j = 0; j < 6; j++){
+ edge_colors.push_back(0.8);
+ edge_colors.push_back(0.8);
+ edge_colors.push_back(0.8);
+ }
+ }
+ } else {
+ netcdf_mpas_read_field(filename, edge_field, edge_values, cur_time, cur_level);
+
+ if(color_bar == 1){
+ min = ranges[edge_field][0];
+ max = ranges[edge_field][1];
+ } else {
+ r8vec_min_max(nedges, edge_values, min, max);
+ }
+
+ cout << "Min: " << min << " Max: " << max << endl;
+
+ for(i = 0; i < nedges; i++){
+ if(max-min != 0.0){
+ h = (edge_values[i] - min)/(max-min)*range_factor;
+ } else {
+ h = (edge_values[i] - min)/1.0 * range_factor;
+ }
+
+ hsv_to_rgb(h, s, v, r, g, b);
+
+ for(j = 0; j < 6; j++){
+ edge_colors.push_back(r);
+ edge_colors.push_back(g);
+ edge_colors.push_back(b);
+ }
+ }
+ }
+}/*}}}*/
+
+void arrowKeys( int a_keys, int x, int y ) {/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // ARROWKEYS processes special keys, for example the arrow keys, and the f-keys
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff
+ //
+
+ switch( a_keys ) {
+ case GLUT_KEY_UP:
+ projElevation = (double)( ( (int)projElevation + 6 ) % 360 );
+ break;
+
+ case GLUT_KEY_DOWN:
+ projElevation = (double) ( ( (int)projElevation - 6 ) % 360 );
+ break;
+
+ case GLUT_KEY_RIGHT:
+ projAzimuth = (double) ( ( (int)projAzimuth + 6 ) % 360 );
+ break;
+
+ case GLUT_KEY_LEFT:
+ projAzimuth = (double) ( ( (int)projAzimuth - 6 ) % 360 );
+ break;
+
+ case GLUT_KEY_F1:
+ glutFullScreen();
+ break;
+
+ case GLUT_KEY_F2:
+ glutReshapeWindow( kWindowWidth, kWindowHeight );
+ break;
+
+ default:
+ break;
+
+ }
+
+ return;
+}/*}}}*/
+void keyPressed( unsigned char key, int x, int y ) {/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // KEYPRESSED processes non-special keys, for example the letter keys
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+
+ usleep( 100 );
+
+ // Used to get codes of keys you don't know.
+// cout << "Key pressed: " << key << endl;
+// cout << "Code is: " << (int)key << endl;
+
+ switch( key ) {
+
+ case KEY_ESCAPE:
+ glutDestroyWindow( window );
+ glDisableClientState( GL_VERTEX_ARRAY );
+ glDisableClientState( GL_COLOR_ARRAY );
+ //glDisableClientState( GL_NORMAL_ARRAY );
+ delete [] nedgesoncell;
+
+ cout << "End of normal execution. Exiting." << endl;
+
+ exit( 0 );
+ case KEY_l:
+ cur_level = (cur_level + 1) % nvertlevels;
+ cout << "Current vertical level: " << cur_level+1 << " out of " << nvertlevels << endl;
+ color_mesh();
+ break;
+ case KEY_L:
+ cur_level = (cur_level - 1) % nvertlevels;
+ if(cur_level < 0)
+ cur_level = nvertlevels - 1;
+ cout << "Current vertical level: " << cur_level+1 << " out of " << nvertlevels << endl;
+ color_mesh();
+ break;
+ case KEY_t:
+ cur_time = (cur_time + 1) % ntime;
+ cout << "Current time level: " << cur_time+1 << " out of " << ntime << endl;
+ color_mesh();
+ break;
+ case KEY_T:
+ cur_time = (cur_time - 1) % ntime;
+ if(cur_time < 0)
+ cur_time = ntime-1;
+ cout << "Current time level: " << cur_time+1 << " out of " << ntime << endl;
+ color_mesh();
+ break;
+ case KEY_s:
+ projLeftRight -= 0.1;
+ break;
+ case KEY_f:
+ projLeftRight += 0.1;
+ break;
+ case KEY_d:
+ projUpDown -= 0.1;
+ break;
+ case KEY_e:
+ projUpDown += 0.1;
+ break;
+ case KEY_b:
+ color_bar = (color_bar + 1) % 2;
+ switch(color_bar){
+ case 0:
+ cout << "Coloring based on current time/vertlevel slice" << endl;
+ break;
+ case 1:
+ cout << "Coloring base on entire time range" << endl;
+ break;
+ default:
+ break;
+ }
+ if(color_bar == 1){
+ switch(drawing){
+ case 0:
+ build_range(vertex_field);
+ cout << " Range of values: Min = " << ranges.at(vertex_field).at(0) << ", Max = " << ranges.at(vertex_field).at(1) << endl;
+ break;
+ case 1:
+ build_range(cell_field);
+ cout << " Range of values: Min = " << ranges.at(cell_field).at(0) << ", Max = " << ranges.at(cell_field).at(1) << endl;
+ break;
+ case 2:
+ build_range(edge_field);
+ cout << " Range of values: Min = " << ranges.at(edge_field).at(0) << ", Max = " << ranges.at(edge_field).at(1) << endl;
+ break;
+ default:
+ break;
+ }
+ }
+
+ break;
+ case KEY_w:
+ draw_lines = (draw_lines + 1) % 2;
+ break;
+ case KEY_c:
+ drawing = ( drawing + 1 ) % 3;
+ break;
+ case KEY_COMMA:
+ projDistance += 0.1;
+ break;
+ case KEY_PERIOD:
+ projDistance -= 0.1;
+ break;
+ case KEY_v:
+ switch(drawing){
+ case 0:
+ vertex_field = netcdf_mpas_list_nvertex_fields(filename);
+ netcdf_mpas_print_field_info(filename, vertex_field);
+ if(color_bar == 1){
+ build_range(vertex_field);
+ cout << " Range of values: Min = " << ranges.at(vertex_field).at(0) << ", Max = " << ranges.at(vertex_field).at(1) << endl;
+ }
+ break;
+ case 1:
+ cell_field = netcdf_mpas_list_ncell_fields(filename);
+ netcdf_mpas_print_field_info(filename, cell_field);
+ if(color_bar == 1){
+ build_range(cell_field);
+ cout << " Range of values: Min = " << ranges.at(cell_field).at(0) << ", Max = " << ranges.at(cell_field).at(1) << endl;
+ }
+ break;
+ case 2:
+ edge_field = netcdf_mpas_list_nedge_fields(filename);
+ netcdf_mpas_print_field_info(filename, edge_field);
+ if(color_bar == 1){
+ build_range(edge_field);
+ cout << " Range of values: Min = " << ranges.at(edge_field).at(0) << ", Max = " << ranges.at(edge_field).at(1) << endl;
+ }
+ break;
+ default:
+ break;
+ }
+ color_mesh();
+ break;
+ case KEY_R:
+ cout << " Trying regions. " << endl;
+ if(!regions_built){
+ build_regions();
+ }
+
+ if(regions_built){
+ region_draw = !region_draw;
+ }
+ break;
+ case KEY_r:
+ cout << "Resetting to default parameters." << endl;
+ cur_time = 0;
+ cur_level = 0;
+ cell_field = 0;
+ edge_field = 0;
+ vertex_field = 0;
+ color_bar = 0;
+ projDistance = 3.0;
+ projUpDown = 0.0;
+ projLeftRight = 0.0;
+ projTwist = 0.0;
+ projElevation = 0.0;
+ projAzimuth = 0.0;
+
+ color_mesh();
+ break;
+ case KEY_0:
+ range_factor = 1.0;
+ color_mesh();
+ break;
+ case KEY_1:
+ range_factor = 0.1;
+ color_mesh();
+ break;
+ case KEY_2:
+ range_factor = 0.2;
+ color_mesh();
+ break;
+ case KEY_3:
+ range_factor = 0.3;
+ color_mesh();
+ break;
+ case KEY_4:
+ range_factor = 0.4;
+ color_mesh();
+ break;
+ case KEY_5:
+ range_factor = 0.5;
+ color_mesh();
+ break;
+ case KEY_6:
+ range_factor = 0.6;
+ color_mesh();
+ break;
+ case KEY_7:
+ range_factor = 0.7;
+ color_mesh();
+ break;
+ case KEY_8:
+ range_factor = 0.8;
+ color_mesh();
+ break;
+ case KEY_9:
+ range_factor = 0.9;
+ color_mesh();
+ break;
+ default:
+ break;
+ }
+
+ return;
+}/*}}}*/
+
+void translateView ( double updown, double leftright){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // TRANSLATEVIEW Translates the drawing left, right, up, or down
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+
+ glTranslated ( leftright, updown, 0.0);
+}/*}}}*/
+void polarView( double distance, double twist, double elevation, double azimuth ) {/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // POLARVIEW Rotates the drawing (in two directions), and translates in and out
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff
+ //
+
+ glTranslated( 0.0, 0.0, -1.0 * distance );
+ glRotated( -1.0 * twist, 0.0, 0.0, 1.0 );
+ glRotated( -1.0 * elevation, 1.0, 0.0, 0.0 );
+ glRotated( azimuth, 0.0, 0.0, 1.0 );
+}/*}}}*/
+
+void hsv_to_rgb(float h, float s, float v, float& r, float& g, float& b){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // HSV_TO_RGB convertes a hue-saturation-value triplet to a red-green-blue triplet
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Geoff Womeldorff, Doug Jacobsen
+ //
+
+ int i;
+ float f;
+ float m;
+ float n;
+
+ h *= 6.0;
+
+ i = floor( h );
+ f = h - i;
+
+ if(!(i & 1)) f = 1 - f; // if i is even
+ m = v * (1 - s);
+ n = v * (1 - s * f);
+
+ switch (i) {
+ case 6:
+ case 0:
+ r = v; g = n; b = m;
+ break;
+ case 1:
+ r = n; g = v; b = m;
+ break;
+ case 2:
+ r = m; g = v; b = n;
+ break;
+ case 3:
+ r = m; g = n; b = v;
+ break;
+ case 4:
+ r = n; g = m; b = v;
+ break;
+ case 5:
+ r = v; g = m; b = n;
+ break;
+ }
+
+ return;
+}/*}}}*/
+
+double getLat(double x, double y, double z){/*{{{*/
+ return asin(z);
+}/*}}}*/
+double getLon(double x, double y, double z){/*{{{*/
+ double lon;
+
+ lon = atan2(y,x);
+
+ if(lon < 0){
+ return 2.0 * M_PI + lon;
+ } else {
+ return lon;
+ }
+}/*}}}*/
+double getX(double lat, double lon){/*{{{*/
+ return cos(lon) * cos(lat);
+}/*}}}*/
+double getY(double lat, double lon){/*{{{*/
+ return sin(lon) * cos(lat);
+}/*}}}*/
+double getZ(double lat, double lon){/*{{{*/
+ return sin(lat);
+}/*}}}*/
+
Added: branches/tools/mpas_draw/netcdf_utils.cpp
===================================================================
--- branches/tools/mpas_draw/netcdf_utils.cpp (rev 0)
+++ branches/tools/mpas_draw/netcdf_utils.cpp 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,1496 @@
+#include <string>
+#include <netcdfcpp.h>
+#include <cstdlib>
+#include <vector>
+
+using namespace std;
+
+
+bool netcdf_mpas_read_onsphere(string filename){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_ONSPHERE determines if points are on a sphere, or in a plane
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 29 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Output, int NETCDF_MPAS_READ_ONSPHERE, true or false based on the global
+ // attribute, on_a_sphere
+ //
+ NcAtt *att_id;
+ NcValues *vals;
+ bool valid;
+ string tmp_name;
+ string sph_name = "on_a_sphere";
+ //
+ // Open the file.
+ //
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ // Get Ntime, which is a NETCDF dimension.
+ //
+ valid = false;
+
+ for(int i = 0; i < ncid.num_atts() && !valid; i++){
+ tmp_name = ncid.get_att(i)->name();
+
+ if(sph_name == tmp_name){
+ valid = true;
+ }
+ }
+
+ if(valid){
+ att_id = ncid.get_att(sph_name.c_str());
+
+ vals = att_id -> values();
+
+ sph_name = "YES ";
+
+ for(int i = 0; i < vals -> num(); i++){
+ tmp_name = vals -> as_string(i);
+
+ if(tmp_name == sph_name){
+ return true;
+ }
+ }
+ } else {
+ return true;
+ }
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return false;
+
+}/*}}}*/
+//****************************************************************************80
+int netcdf_mpas_read_dim ( string filename, string dim_name ){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_DIM gets the size of the dimension with name dim_name
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 29 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, string DIM_NAME, the name of the dimension in question
+ //
+ // Output, int NETCDF_MPAS_READ_DIM, the value of the dimension.
+ //
+ int ntime;
+ long int dim_size;
+ NcDim *dim_id;
+ bool valid;
+ string tmp_name;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get Ntime, which is a NETCDF dimension.
+ //
+ valid = false;
+
+ for(int i = 0; i < ncid.num_dims() && !valid; i++){
+ tmp_name = ncid.get_dim(i)->name();
+
+ if(dim_name == tmp_name){
+ valid = true;
+ }
+ }
+
+ if(valid){
+ dim_id = ncid.get_dim(dim_name.c_str());
+ dim_size = (*dim_id).size ( );
+ } else {
+ dim_size = 1;
+ }
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ ntime = ( int ) dim_size;
+
+ return ntime;
+}/*}}}*/
+//****************************************************************************80
+int netcdf_mpas_read_num_vars(string filename){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_NUM_VARS gets the number of variables in the netcdf file
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 29 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Output, int NETCDF_MPAS_READ_NUM_VARS, the number of variables
+ //
+ long int var_size;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get Ntime, which is a NETCDF dimension.
+ //
+ var_size = ncid.num_vars();
+
+ ncid.close();
+
+ return var_size;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_xyzcell ( string filename, int ncells, double xcell[], double ycell[], double zcell[] ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_XYZCELL reads xCell, yCell, zCell.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 29 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int NCELLS, the number of nodes.
+ //
+ // Output, double XCELL[NCELLS], YCELL[NCELLS], ZCELL[NCELLS], the
+ // coordinates of the nodes.
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "xCell" );
+ (*var_id).get ( &xcell[0], ncells );
+ var_id = ncid.get_var ( "yCell" );
+ (*var_id).get ( &ycell[0], ncells );
+ var_id = ncid.get_var ( "zCell" );
+ (*var_id).get ( &zcell[0], ncells );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_xyzvertex ( string filename, int nvertices, double xvertex[], double yvertex[], double zvertex[] ){/*{{{*/
+ //****************************************************************************80
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_CELLS gets the cell center coordinates.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int NVERTICES, the number of vertices.
+ //
+ // Output, double XVERTEX[NVERTICES], YVERTEXL[NVERTICES],
+ // ZVERTEX[NVERTICES], the coordinates of the nodes.
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "xVertex" );
+ (*var_id).get ( &xvertex[0], nvertices );
+ var_id = ncid.get_var ( "yVertex" );
+ (*var_id).get ( &yvertex[0], nvertices );
+ var_id = ncid.get_var ( "zVertex" );
+ (*var_id).get ( &zvertex[0], nvertices );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_xyzedge ( string filename, int nedges, double xedge[], double yedge[], double zedge[] ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_EDGES gets the edge midpoint coordinates.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int NEDGES, the number of edges.
+ //
+ // Output, double XEDGE[NEDGES], YEDGE[NEDGES],
+ // ZEDGE[NEDGES], the coordinates of the edge midpoints.
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "xEdge" );
+ (*var_id).get ( &xedge[0], nedges );
+ var_id = ncid.get_var ( "yEdge" );
+ (*var_id).get ( &yedge[0], nedges );
+ var_id = ncid.get_var ( "zEdge" );
+ (*var_id).get ( &zedge[0], nedges );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_verticesoncell ( string filename, int maxedges, int ncells, int verticesOnCell[] ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_VERTICESONCELLS gets verticesOnCells.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int MAXEDGES, the maximum number of edges for a cell.
+ //
+ // Input, int NCELLS, the number of cells.
+ //
+ // Output, int VERTICESONCELLS[MAXEDGES*NCELLS];
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "verticesOnCell" );
+ (*var_id).get ( &verticesOnCell[0], ncells, maxedges );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_verticesonedge ( string filename, int nedges, int verticesonedge[] ){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_VERTICESONEDGE gets the verticesOnEdge information.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int NEDGES, the number of edges.
+ //
+ // Output, int VERTICESONEDGE[2*NEDGES];
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "verticesOnEdge" );
+ (*var_id).get ( &verticesonedge[0], nedges, 2 );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_cellsonvertex ( string filename, int nvertices, int cellsonvertex[] ){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_CELLSONVERTEX gets the cellsOnVertex information.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int NEDGES, the number of edges.
+ //
+ // Output, int CELLSONVERTEX[3*NVERTICES];
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "cellsOnVertex" );
+ (*var_id).get ( &cellsonvertex[0], nvertices, 3 );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_cellsonedge ( string filename, int nedges, int cellsonedge[] ){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_CELLSONEDGE gets the cellsOnEdge information.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int NEDGES, the number of edges.
+ //
+ // Output, int CELLSONEDGE[2*NEDGES];
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "cellsOnEdge" );
+ (*var_id).get ( &cellsonedge[0], nedges, 2 );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_nedgesoncell ( string filename, int ncells, int nedgesoncell[] ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_NEDGESONCELLS gets nedgesOnCells.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int NCELLS, the number of cells.
+ //
+ // Output, int NEDGESONCELLS[NCELLS];
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var ( "nEdgesOnCell" );
+ (*var_id).get ( &nedgesoncell[0], ncells );
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+int netcdf_mpas_list_ncell_fields(string filename){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_LIST_NCELL_FIELDS lists all fields that are printable based on nCells
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ NcVar *var_id;
+ NcDim *ncell_dim, *nvert_dim, *nedge_dim;
+ NcToken ncell_tok, nvert_tok, nedge_tok;
+ NcToken dim_name;
+ vector<NcToken> var_names;
+ vector<int> var_ids;
+ string name;
+ int included;
+ int has_ncells;
+ int num_vars;
+ int num_dims;
+ int i;
+ int j;
+ int valid = 0;
+ int choice = 0;
+
+ vector<string> excluded_vars;
+
+ excluded_vars.push_back("xCell");
+ excluded_vars.push_back("yCell");
+ excluded_vars.push_back("zCell");
+ excluded_vars.push_back("latCell");
+ excluded_vars.push_back("lonCell");
+ excluded_vars.push_back("indexToCellID");
+ excluded_vars.push_back("edgesOnCell");
+ excluded_vars.push_back("cellsOnCell");
+ excluded_vars.push_back("localVerticalUnitVectors");
+ excluded_vars.push_back("verticesOnCell");
+
+ num_vars = ncid.num_vars();
+ ncell_dim = ncid.get_dim("nCells");
+ ncell_tok = ncell_dim->name();
+ nvert_dim = ncid.get_dim("nVertices");
+ nvert_tok = nvert_dim->name();
+ nedge_dim = ncid.get_dim("nEdges");
+ nedge_tok = nedge_dim->name();
+
+ for(i = 0; i < num_vars; i++){
+ var_id = ncid.get_var(i);
+ num_dims = var_id->num_dims();
+ included = 1;
+ name = var_id->name();
+ for(j = 0; j < excluded_vars.size(); j++){
+ if(excluded_vars.at(j) == name){
+ included = 0;
+ }
+ }
+
+ if(included){
+ has_ncells = 0;
+
+ for(j = 0; j < num_dims; j++){
+ dim_name = (var_id->get_dim(j))->name();
+ if(ncell_tok == dim_name){
+ has_ncells = 1;
+ }
+ }
+
+ if(has_ncells){
+ var_names.push_back(var_id->name());
+ var_ids.push_back(i);
+ }
+ }
+ }
+
+ valid = 0;
+
+ if(var_names.empty()){
+ cout << endl;
+
+ cout << "No variables indexed by nCells to color with." << endl;
+ }
+
+ while(!valid && !var_names.empty()){
+ cout << endl;
+ cout << "Available variables indexed by nCells." << endl;
+ for(i = 0; i < var_names.size(); i++){
+ cout << var_ids.at(i) << "\t" << var_names.at(i) << ":\t";
+ num_dims = ncid.get_var(var_ids.at(i))->num_dims();
+ for(j = 0; j < num_dims-1; j++){
+ cout << ncid.get_var(var_ids.at(i))->get_dim(j)->name() << "*";
+ }
+ cout << ncid.get_var(var_ids.at(i))->get_dim(num_dims-1)->name() << endl;
+ }
+ cout << endl << endl;
+
+ cout << "Enter the number of the field you would like to print" << endl;
+ cout << "Choose from list above, or enter 0 for default:" << endl;
+ cin >> choice;
+
+ for(i = 0; i < var_ids.size(); i++){
+ if(var_ids.at(i) == choice){
+ valid = 1;
+ }
+ }
+
+ if(choice == 0){
+ valid = 1;
+ }
+ }
+
+ return choice;
+}/*}}}*/
+//****************************************************************************80
+int netcdf_mpas_list_nvertex_fields(string filename){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_LIST_NVERTEX_FIELDS lists all fields that are printable based on nVertices
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ NcVar *var_id;
+ NcDim *ncell_dim, *nvert_dim, *nedge_dim;
+ NcToken ncell_tok, nvert_tok, nedge_tok;
+ NcToken dim_name;
+ string name;
+ vector<NcToken> var_names;
+ vector<int> var_ids;
+ int included;
+ int has_nverts;
+ int num_vars;
+ int num_dims;
+ int i;
+ int j;
+ int choice = 0;
+ int valid = 0;
+
+ vector<string> excluded_vars;
+
+ excluded_vars.push_back("xVertex");
+ excluded_vars.push_back("yVertex");
+ excluded_vars.push_back("zVertex");
+ excluded_vars.push_back("latVertex");
+ excluded_vars.push_back("lonVertex");
+ excluded_vars.push_back("indexToVertexID");
+ excluded_vars.push_back("edgesOnVertex");
+ excluded_vars.push_back("cellsOnVertex");
+ excluded_vars.push_back("kiteAreasOnVertex");
+
+ excluded_vars.push_back("verticesOnVertex");
+ excluded_vars.push_back("kiteAreaVertex");
+
+ num_vars = ncid.num_vars();
+ ncell_dim = ncid.get_dim("nCells");
+ ncell_tok = ncell_dim->name();
+ nvert_dim = ncid.get_dim("nVertices");
+ nvert_tok = nvert_dim->name();
+ nedge_dim = ncid.get_dim("nEdges");
+ nedge_tok = nedge_dim->name();
+
+ for(i = 0; i < num_vars; i++){
+ var_id = ncid.get_var(i);
+ num_dims = var_id->num_dims();
+ included = 1;
+ name = var_id->name();
+ for(j = 0; j < excluded_vars.size(); j++){
+ if(excluded_vars.at(j) == name){
+ included = 0;
+ }
+ }
+
+ if(included){
+ has_nverts = 0;
+
+ for(j = 0; j < num_dims; j++){
+ dim_name = (var_id->get_dim(j))->name();
+ if(nvert_tok == dim_name){
+ has_nverts = 1;
+ }
+ }
+
+ if(has_nverts){
+ var_names.push_back(var_id->name());
+ var_ids.push_back(i);
+ }
+ }
+ }
+
+ valid = 0;
+
+ if(var_names.empty()){
+ cout << endl;
+ cout << "No variables indexed by nVertices to color with." << endl;
+ }
+
+ while(!valid && !var_names.empty()){
+ cout << endl;
+ cout << "Available variables indexed by nVertices." << endl;
+ for(i = 0; i < var_names.size(); i++){
+ cout << var_ids.at(i) << "\t" << var_names.at(i) << ":\t";
+ num_dims = ncid.get_var(var_ids.at(i))->num_dims();
+ for(j = 0; j < num_dims-1; j++){
+ cout << ncid.get_var(var_ids.at(i))->get_dim(j)->name() << "*";
+ }
+ cout << ncid.get_var(var_ids.at(i))->get_dim(num_dims-1)->name() << endl;
+ }
+ cout << endl << endl;
+
+ cout << "Enter the number of the field you would like to print" << endl;
+ cout << "Choose from list above, or enter 0 for default:" << endl;
+ cin >> choice;
+
+ for(i = 0; i < var_ids.size(); i++){
+ if(var_ids.at(i) == choice){
+ valid = 1;
+ }
+ }
+
+ if(choice == 0){
+ valid = 1;
+ }
+ }
+
+ return choice;
+
+}/*}}}*/
+//****************************************************************************80
+int netcdf_mpas_list_nedge_fields(string filename){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_LIST_NEDGE_FIELDS lists all fields that are printable based on nEdges
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 30 December 2010
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ NcVar *var_id;
+ NcDim *ncell_dim, *nvert_dim, *nedge_dim;
+ NcToken ncell_tok, nvert_tok, nedge_tok;
+ NcToken dim_name;
+ string name;
+ vector<NcToken> var_names;
+ vector<int> var_ids;
+ int included;
+ int has_nedges;
+ int num_vars;
+ int num_dims;
+ int i;
+ int j;
+ int choice = 0;
+ int valid = 0;
+
+ vector<string> excluded_vars;
+
+ excluded_vars.push_back("xEdge");
+ excluded_vars.push_back("yEdge");
+ excluded_vars.push_back("zEdge");
+ excluded_vars.push_back("latEdge");
+ excluded_vars.push_back("lonEdge");
+ excluded_vars.push_back("indexToEdgeID");
+ excluded_vars.push_back("cellsOnEdge");
+ excluded_vars.push_back("edgesOnEdge");
+ excluded_vars.push_back("edgeNormalVectors");
+ excluded_vars.push_back("cellTangentPlane");
+ excluded_vars.push_back("weightsOnEdge");
+ excluded_vars.push_back("verticesOnEdge");
+
+ num_vars = ncid.num_vars();
+ ncell_dim = ncid.get_dim("nCells");
+ ncell_tok = ncell_dim->name();
+ nvert_dim = ncid.get_dim("nVertices");
+ nvert_tok = nvert_dim->name();
+ nedge_dim = ncid.get_dim("nEdges");
+ nedge_tok = nedge_dim->name();
+
+ for(i = 0; i < num_vars; i++){
+ var_id = ncid.get_var(i);
+ num_dims = var_id->num_dims();
+ included = 1;
+ name = var_id->name();
+ for(j = 0; j < excluded_vars.size(); j++){
+ if(excluded_vars.at(j) == name){
+ included = 0;
+ }
+ }
+
+ if(included){
+ has_nedges = 0;
+
+ for(j = 0; j < num_dims; j++){
+ dim_name = (var_id->get_dim(j))->name();
+ if(nedge_tok == dim_name){
+ has_nedges = 1;
+ }
+ }
+
+ if(has_nedges){
+ var_names.push_back(var_id->name());
+ var_ids.push_back(i);
+ }
+ }
+ }
+
+ if(var_names.empty()){
+ cout << endl;
+ cout << "No variables indexed by nEdges to color with." << endl;
+ }
+
+ while(!valid && !var_names.empty()){
+ cout << endl;
+ cout << "Available variables indexed by nEdges." << endl;
+ for(i = 0; i < var_names.size(); i++){
+ cout << var_ids.at(i) << "\t" << var_names.at(i) << ":\t";
+ num_dims = ncid.get_var(var_ids.at(i))->num_dims();
+ for(j = 0; j < num_dims-1; j++){
+ cout << ncid.get_var(var_ids.at(i))->get_dim(j)->name() << "*";
+ }
+ cout << ncid.get_var(var_ids.at(i))->get_dim(num_dims-1)->name() << endl;
+ }
+ cout << endl << endl;
+ cout << "Enter the number of the field you would like to print" << endl;
+ cout << "Choose from list above, or enter 0 for default:" << endl;
+ cin >> choice;
+
+ for(i = 0; i < var_ids.size(); i++){
+ if(var_ids.at(i) == choice){
+ valid = 1;
+ }
+ }
+
+ if(choice == 0){
+ valid = 1;
+ }
+ }
+ return choice;
+}/*}}}*/
+//****************************************************************************80
+int netcdf_mpas_field_num_dims(string filename, int id){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_FIELD_NUM_DIMS gets the number of dimensions for a field
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int id, the variable id
+ //
+ // Output, double dims, the number of dimensions for the field
+ //
+ NcVar *var_id;
+ int num_dims;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var (id);
+ if(!var_id->is_valid()){
+ cout << "Field " << id << " doesn't exist." << endl;
+ cout << "Exiting." << endl;
+ exit(1);
+ }
+ num_dims = var_id->num_dims();
+ ncid.close ( );
+
+ return num_dims;
+}/*}}}*/
+//****************************************************************************80
+int netcdf_mpas_field_num_items(string filename, int id){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_FIELD_NUM_items gets the number of items for a field
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int id, the variable id
+ //
+ // Output, int num_items, the number of items for the field
+ //
+ NcVar *var_id;
+ int num_dims;
+ int num_items;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var (id);
+ if(!var_id->is_valid()){
+ cout << "Field " << id << " doesn't exist." << endl;
+ cout << "Exiting." << endl;
+ exit(1);
+ }
+ num_dims = var_id->num_dims();
+
+ num_items = var_id->get_dim(0)->size();
+
+ for(int i = 1; i < num_dims; i++){
+ num_items *= var_id->get_dim(i)->size();
+ }
+ ncid.close ( );
+
+ return num_items;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_print_field_info(string filename, int id){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_PRINT_FIELD_INFO prints a fields information, like dimensions and name
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // John Burkardt, Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int id, the variable id
+ //
+ NcVar *var_id;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var (id);
+
+ if(!var_id->is_valid()){
+ cout << "netcdf_mpas_print_field_info: field " << id << " doesn't exist." << endl;
+ cout << "Exiting." << endl;
+ exit(1);
+ }
+
+ cout << endl << endl;
+ cout << "Field " << id << " is " << var_id->name() << endl;
+ cout << "Dimensions: ";
+ for(int i = 0; i < var_id->num_dims()-1; i++){
+ cout << var_id->get_dim(i)->name() << "*";
+ }
+ cout << var_id->get_dim(var_id->num_dims()-1)->name() << endl;
+ cout << "Variable type: " << var_id->type() << endl;
+
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_field(string filename, int id, double values[], int cur_time, int cur_level){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_1D_FIELD gets a field based on one dimension.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int id, the variable id
+ //
+ // Output, double values[] are the values of the field from the netcdf file
+ //
+ // Output, long dims[] the dimensions of the field from the netcdf file
+ NcVar *var_id;
+ int num_dims;
+ int type;
+ long *dims;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var (id);
+ num_dims = var_id->num_dims();
+ dims = new long[num_dims];
+ for(int i = 0; i < num_dims; i++){
+ dims[i] = var_id->get_dim(i)->size();
+ }
+
+ type = var_id->type();
+ if(type == 6){ // NCDOUBLE
+ switch(num_dims){
+ case 2:
+ (*var_id).set_cur(0, cur_level);
+ (*var_id).get ( &values[0], dims[0], 1);
+ break;
+ case 3:
+ (*var_id).set_cur(cur_time, 0, cur_level);
+ (*var_id).get ( &values[0], 1, dims[1], 1);
+ break;
+ default: //1 dim
+ (*var_id).get ( &values[0], dims[0]);
+ break;
+ }
+ } else if(type == 4){ //NCINT
+ int *tmp_values;
+ int num_items;
+
+ switch(num_dims){
+ case 2:
+ num_items = dims[0];
+ tmp_values = new int[num_items];
+ (*var_id).set_cur(0, cur_level);
+ (*var_id).get ( &tmp_values[0], dims[0], 1);
+ break;
+ case 3:
+ num_items = dims[1];
+ tmp_values = new int[num_items];
+ (*var_id).set_cur(cur_time, 0, cur_level);
+ (*var_id).get ( &tmp_values[0], 1, dims[1], 1);
+ break;
+ default: //1 dim
+ num_items = dims[0];
+ tmp_values = new int[num_items];
+ (*var_id).get ( &tmp_values[0], dims[0]);
+ break;
+ }
+
+ for(int i = 0; i < num_items; i++){
+ values[i] = tmp_values[i]*1.0;
+ }
+
+ delete [] tmp_values;
+ }
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
+//****************************************************************************80
+void netcdf_mpas_read_full_field(string filename, int id, double values[]){/*{{{*/
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // NETCDF_MPAS_READ_FULL_FIELD gets a full field.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 01 January 2011
+ //
+ // Author:
+ //
+ // Doug Jacobsen
+ //
+ // Reference:
+ //
+ // Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
+ // The NETCDF User's Guide,
+ // Unidata Program Center, March 2009.
+ //
+ // Parameters:
+ //
+ // Input, string NC_FILENAME, the name of the NETCDF file to examine.
+ //
+ // Input, int id, the variable id
+ //
+ // Output, double values[] are the values of the field from the netcdf file
+ //
+ // Output, long dims[] the dimensions of the field from the netcdf file
+ NcVar *var_id;
+ int num_dims;
+ int type;
+ long *dims;
+ //
+ // Open the file.
+ #ifdef _64BITOFFSET
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly, NULL, 0, NcFile::Offset64Bits );
+ #else
+ NcFile ncid ( filename.c_str ( ), NcFile::ReadOnly );
+ #endif
+ //
+ //
+ // Get the variable values.
+ //
+ var_id = ncid.get_var (id);
+ num_dims = var_id->num_dims();
+ dims = new long[num_dims];
+ for(int i = 0; i < num_dims; i++){
+ dims[i] = var_id->get_dim(i)->size();
+ }
+
+ type = var_id->type();
+ if(type == 6){ // NCDOUBLE
+ switch(num_dims){
+ case 2:
+ (*var_id).get ( &values[0], dims[0], dims[1]);
+ break;
+ case 3:
+ (*var_id).get ( &values[0], dims[0], dims[1], dims[2]);
+ break;
+ default: //1 dim
+ (*var_id).get ( &values[0], dims[0]);
+ break;
+ }
+ } else if(type == 4){ //NCINT
+ int *tmp_values;
+ int num_items;
+
+ switch(num_dims){
+ case 2:
+ num_items = dims[0]*dims[1];
+ tmp_values = new int[num_items];
+ (*var_id).get ( &tmp_values[0], dims[0], dims[1]);
+ break;
+ case 3:
+ num_items = dims[0]*dims[1]*dims[2];
+ tmp_values = new int[num_items];
+ (*var_id).get ( &tmp_values[0], dims[0], dims[1], dims[2]);
+ break;
+ default: //1 dim
+ num_items = dims[0];
+ tmp_values = new int[num_items];
+ (*var_id).get ( &tmp_values[0], dims[0]);
+ break;
+ }
+
+ for(int i = 0; i < num_items; i++){
+ values[i] = tmp_values[i]*1.0;
+ }
+
+ delete [] tmp_values;
+ }
+ //
+ // Close the file.
+ //
+ ncid.close ( );
+
+ return;
+}/*}}}*/
Added: branches/tools/mpas_draw/netcdf_utils.h
===================================================================
--- branches/tools/mpas_draw/netcdf_utils.h (rev 0)
+++ branches/tools/mpas_draw/netcdf_utils.h 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,22 @@
+#include <string>
+#include <netcdfcpp.h>
+
+bool netcdf_mpas_read_onsphere(string filename);
+int netcdf_mpas_read_dim ( string filename, string dim_name );
+int netcdf_mpas_read_num_vars(string filename);
+void netcdf_mpas_read_xyzcell ( string filename, int ncells, double xcell[], double ycell[], double zcell[] );
+void netcdf_mpas_read_xyzvertex ( string filename, int nvertices, double xvertex[], double yvertex[], double zvertex[] );
+void netcdf_mpas_read_xyzedge ( string filename, int nedges, double xedge[], double yedge[], double zedge[] );
+void netcdf_mpas_read_verticesoncell ( string filename, int maxedges, int ncells, int verticesOnCell[] );
+void netcdf_mpas_read_verticesonedge ( string filename, int nedges, int verticesonedge[] );
+void netcdf_mpas_read_cellsonvertex ( string filename, int nvertices, int cellsonvertex[] );
+void netcdf_mpas_read_cellsonedge ( string filename, int nedges, int cellsonedge[] );
+void netcdf_mpas_read_nedgesoncell ( string filename, int ncells, int nedgesoncell[] );
+int netcdf_mpas_list_ncell_fields(string filename);
+int netcdf_mpas_list_nvertex_fields(string filename);
+int netcdf_mpas_list_nedge_fields(string filename);
+int netcdf_mpas_field_num_dims(string filename, int id);
+int netcdf_mpas_field_num_items(string filename, int id);
+void netcdf_mpas_print_field_info(string filename, int id);
+void netcdf_mpas_read_field(string filename, int id, double values[], int cur_time, int cur_level);
+void netcdf_mpas_read_full_field(string filename, int id, double values[]);
Added: branches/tools/mpas_draw/vec_utils.cpp
===================================================================
--- branches/tools/mpas_draw/vec_utils.cpp (rev 0)
+++ branches/tools/mpas_draw/vec_utils.cpp 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,274 @@
+#include <algorithm>
+
+using namespace std;
+
+double r8_huge ( ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // R8_HUGE returns a "huge" R8.
+ //
+ // Discussion:
+ //
+ // The value returned by this function is NOT required to be the
+ // maximum representable R8. This value varies from machine to machine,
+ // from compiler to compiler, and may cause problems when being printed.
+ // We simply want a "very large" but non-infinite number.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 06 October 2007
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Parameters:
+ //
+ // Output, double R8_HUGE, a "huge" R8 value.
+ //
+ double value;
+
+ value = 1.0E+30;
+
+ return value;
+}/*}}}*/
+//****************************************************************************80
+double r8_max ( double x, double y ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // R8_MAX returns the maximum of two R8's.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 18 August 2004
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Parameters:
+ //
+ // Input, double X, Y, the quantities to compare.
+ //
+ // Output, double R8_MAX, the maximum of X and Y.
+ //
+ double value;
+
+ if ( y < x )
+ {
+ value = x;
+ }
+ else
+ {
+ value = y;
+ }
+ return value;
+}/*}}}*/
+//****************************************************************************80
+double r8_min ( double x, double y ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // R8_MIN returns the minimum of two R8's.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 18 August 2004
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Parameters:
+ //
+ // Input, double X, Y, the quantities to compare.
+ //
+ // Output, double R8_MIN, the minimum of X and Y.
+ //
+ double value;
+
+ if ( y > x )
+ {
+ value = x;
+ }
+ else
+ {
+ value = y;
+ }
+ return value;
+}/*}}}*/
+//****************************************************************************80
+double r8vec_max ( int n, double r8vec[] ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // R8VEC_MAX returns the value of the maximum element in an R8VEC.
+ //
+ // Discussion:
+ //
+ // An R8VEC is a vector of R8's.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 22 August 2010
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Parameters:
+ //
+ // Input, int N, the number of entries in the array.
+ //
+ // Input, double R8VEC[N], a pointer to the first entry of the array.
+ //
+ // Output, double R8VEC_MAX, the value of the maximum element. This
+ // is set to 0.0 if N <= 0.
+ //
+ int i;
+ double value;
+
+ if ( n <= 0 )
+ {
+ value = 0.0;
+ }
+
+ value = r8vec[0];
+
+ for ( i = 1; i < n; i++ )
+ {
+ if(r8vec[i] > value)
+ value = r8vec[i];
+ }
+ return value;
+}/*}}}*/
+//****************************************************************************80
+double r8vec_min ( int n, double r8vec[] ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // R8VEC_MIN returns the value of the minimum element in an R8VEC.
+ //
+ // Discussion:
+ //
+ // An R8VEC is a vector of R8's.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 02 July 2005
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Parameters:
+ //
+ // Input, int N, the number of entries in the array.
+ //
+ // Input, double R8VEC[N], the array to be checked.
+ //
+ // Output, double R8VEC_MIN, the value of the minimum element.
+ //
+ int i;
+ double value;
+
+ value = r8_huge ( );
+
+ if ( n <= 0 )
+ {
+ return value;
+ }
+
+ for ( i = 0; i < n; i++ )
+ {
+ if(r8vec[i] < value)
+ value = r8vec[i];
+ }
+ return value;
+}/*}}}*/
+//****************************************************************************80
+void r8vec_min_max ( int n, double r8vec[], double &min, double&max ){/*{{{*/
+
+ //****************************************************************************80
+ //
+ // Purpose:
+ //
+ // R8VEC_MIN_MAX returns the value of the minimum and maximum element in an R8VEC.
+ //
+ // Discussion:
+ //
+ // An R8VEC is a vector of R8's.
+ //
+ // Licensing:
+ //
+ // This code is distributed under the GNU LGPL license.
+ //
+ // Modified:
+ //
+ // 02 July 2005
+ //
+ // Author:
+ //
+ // John Burkardt
+ //
+ // Parameters:
+ //
+ // Input, int N, the number of entries in the array.
+ //
+ // Input, double R8VEC[N], the array to be checked.
+ //
+ // Output, double R8VEC_MIN, the value of the minimum element.
+ //
+ // Output, double R8VEC_MAX, the value of the maximum element.
+ //
+ int i;
+
+ if(n <= 0){
+ min = r8_huge();
+ max = 0.0;
+ return;
+ }
+
+ min = r8vec[0];
+ max = r8vec[0];
+
+ for ( i = 1; i < n; i++ )
+ {
+ min = std::min(r8vec[i], min);
+ max = std::max(r8vec[i], max);
+ }
+
+ return;
+}/*}}}*/
Added: branches/tools/mpas_draw/vec_utils.h
===================================================================
--- branches/tools/mpas_draw/vec_utils.h (rev 0)
+++ branches/tools/mpas_draw/vec_utils.h 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,8 @@
+using namespace std;
+
+double r8_huge ( );
+double r8_max ( double x, double y );
+double r8_min ( double x, double y );
+double r8vec_max ( int n, double r8vec[] );
+double r8vec_min ( int n, double r8vec[] );
+void r8vec_min_max ( int n, double r8vec[], double &min, double&max );
Added: branches/tools/points-mpas/Makefile
===================================================================
--- branches/tools/points-mpas/Makefile (rev 0)
+++ branches/tools/points-mpas/Makefile 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,20 @@
+#Assumes ${NETCDF} is defined to the root of the netcdf library
+
+CC=g++
+CFLAGS= -O3
+EXECUTABLE= PointsToMpas.x
+
+NCINCDIR=/usr/include
+NCLIBDIR=/usr/lib
+
+all:
+ ${CC} points-mpas.cpp -I${NCINCDIR} -L${NCLIBDIR} -lnetcdf_c++ ${CFLAGS} -o ${EXECUTABLE}
+
+debug:
+ ${CC} points-mpas.cpp -I${NCINCDIR} -L${NCLIBDIR} -lnetcdf_c++ -g -o ${EXECUTABLE}
+
+clean:
+ rm -f grid.nc
+ rm -f graph.info
+ rm -f ${EXECUTABLE}
+
Added: branches/tools/points-mpas/Params
===================================================================
--- branches/tools/points-mpas/Params (rev 0)
+++ branches/tools/points-mpas/Params 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,14 @@
+Is the input Cartesian or Latitude-Longitude (0 - Cartesian, 1 - Lat-lon)
+0
+Are the triangles base zero or base one? (0 - base 0, 1 - base 1)
+0
+What is the radius of the sphere these points are defined on?
+1.0
+How many vertical levels do you want in the output grid?
+1
+How many tracers do you want in the output grid?
+1
+What was the convergence criteria used to make this grid?
+0.0
+Should this grid be vordraw compatible?
+0
Added: branches/tools/points-mpas/README
===================================================================
--- branches/tools/points-mpas/README (rev 0)
+++ branches/tools/points-mpas/README 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,83 @@
+Readme for points-mpas.cpp
+
+Author: Doug Jacobsen <dwj07@fsu.edu>
+
+Purpose:
+ points-mpas.cpp is a piece of software designed to read in a point set and a
+ triangulation defined over the entire sphere and output a netcdf grid file
+ that is capable of being used in the MPAS modeling framework.
+
+Requirements:
+ points-mpas requires the c++ netcdf libraries to be able to write a MPAS grid file.
+ It also requires the tr1 headers for c++, specifically unordered_set.
+ It has been tested using g++ version 4.4.5
+
+Source:
+ points-mpas has 2 source files. points-mpas.cpp and triangulation.h
+
+ points-mpas.cpp:
+ This file contains all of the main routines for reading the input files, building the grid file, and writing
+ the netcdf and graph files.
+
+ triangulation.h:
+ This file contains the classes used in points-mpas.cpp, as well as some useful triangulation functions.
+
+Input:
+ points-mpas has 3 input files, Params, SaveVertices and SaveTriangles. These will be described below.
+
+ Params:
+ This file contains the parameters used in points-mpas. If there isn't one created simply run the
+ compiled software once, and one will be generated for you with default values for the parameters.
+
+ SaveVertices:
+ This file contains the input vertices, or point set, that is to be converted. The file should only
+ contain the coordinates of each point in the set. Currently points-mpas requires this point set to be defined over the
+ entire sphere.
+
+ Each row represents one vertex. Each column in the file represents a coordinate of that point.
+ For cartesian points the order is x y z, for Lat-lon points the order is lat lon.
+
+ This point set can be given as either Cartesian or Latitude-Longitude.
+
+ SaveTriangles:
+ This file contains the triangulation of the points given in SaveVertices. Each row represents a triangle.
+ Each column represents an index to a point defined in SaveVertices. A triangle is made up of three
+ vertices, and the indices represent the vertices used to make up a triangle.
+
+ The triangle can be base zero (meaning that the lowest index is 0) or base one (meaning lowest index is 1).
+ To change between the two, see the Params file.
+
+Output:
+ points-mpas outputs two files. grid.nc and graph.info
+
+ grid.nc:
+ This is the mpas input grid file. I can be used for simulations in mpas.
+
+ graph.info:
+ This is the graph file used for domain decomposition using a program such as metis.
+ They are used to determine which processor will contain which cells.
+
+ Using this file, one can use metis to build the graph files for an arbitrary number of processors
+ which is needed if one is to use MPI with MPAS.
+
+Notes:
+ The point set and triangulation can be built in a variety of ways.
+ First, the point set needs to be defined, since the triangulation references the point set.
+ As mentioned previously, the point set needs to be defined over the entire sphere currently.
+ At a later point in time, points-mpas might support limited area domains, but currently there
+ are no plans of doing this.
+
+ After the point set has been created, the triangulation needs to be built.
+ To built a triangulation, one can use several available libraries.
+ The triangulation of a point set defined on the sphere is the same thing as
+ the convex hull of that same point set. From this, one can build either the Delaunay triangulation,
+ or the convex hull.
+
+ Some example libraries include: qhull, CGAL, and Stripack.
+
+ qhull can be used from inside matlab, or as a stand alone program.
+ qhull has a library that can be called from c or c++ but it's overly complicated
+ to make work correctly, and it is far easier to call it as an external program or from
+ matlab.
+ CGAL and Stripack are libraries for c++ and FORTRAN respectively.
+
Added: branches/tools/points-mpas/SaveTriangles
===================================================================
--- branches/tools/points-mpas/SaveTriangles (rev 0)
+++ branches/tools/points-mpas/SaveTriangles 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,80 @@
+22 41 31
+12 23 41
+9 23 12
+10 19 18
+7 24 22
+8 22 40
+0 21 33
+16 20 25
+6 30 19
+11 18 30
+5 16 25
+11 23 18
+14 38 27
+0 13 36
+8 12 41
+9 12 29
+18 23 37
+7 32 24
+7 31 15
+5 25 21
+6 27 15
+5 38 14
+17 35 29
+12 39 29
+2 40 24
+28 39 40
+6 14 27
+4 35 17
+4 33 25
+24 32 34
+4 20 35
+6 15 30
+0 36 21
+2 24 34
+15 27 32
+7 15 32
+1 32 27
+9 37 23
+8 41 22
+1 34 32
+13 34 36
+3 39 28
+2 13 28
+7 22 31
+1 27 38
+4 17 33
+3 17 29
+1 38 36
+21 25 33
+18 19 30
+3 28 26
+9 35 37
+2 28 40
+15 31 30
+3 26 17
+5 14 16
+21 36 38
+10 37 20
+10 18 37
+1 36 34
+10 20 16
+2 34 13
+20 37 35
+0 26 13
+0 33 26
+4 25 20
+11 41 23
+14 19 16
+5 21 38
+17 26 33
+13 26 28
+9 29 35
+11 30 31
+3 29 39
+6 19 14
+8 39 12
+11 31 41
+8 40 39
+10 16 19
+22 24 40
Added: branches/tools/points-mpas/SaveVertices
===================================================================
--- branches/tools/points-mpas/SaveVertices (rev 0)
+++ branches/tools/points-mpas/SaveVertices 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,42 @@
+0.1563072040384713 -0.0962758804584656 -0.9830050929713554
+-0.4072397669120782 0.6987500161062452 -0.5881361978631465
+-0.7846511913697883 -0.2819815773113404 -0.5520950080728684
+0.0189018261147366 -0.9325259815061092 -0.3606078407156394
+0.8929363271403706 -0.3538529410719103 -0.2783034526696073
+0.6295663386297677 0.6543310704843933 -0.4189237107911852
+-0.0189018261148352 0.9325259815061259 0.3606078407155913
+-0.8929363271403882 0.3538529410718601 0.2783034526696146
+-0.6295663386296628 -0.6543310704844766 0.4189237107912129
+0.4072397669121988 -0.6987500161062059 0.5881361978631097
+0.7846511913697121 0.2819815773114498 0.5520950080729208
+-0.1563072040385062 0.0962758804584336 0.9830050929713530
+-0.1306802800482702 -0.7953211078538781 0.5919349624820160
+-0.3693313291199989 -0.2223341552467363 -0.9023092001864855
+0.3589395945546560 0.9327311726563537 -0.0342772084050387
+-0.5359650189610209 0.7561145595512165 0.3755426357749311
+0.8312562076671715 0.5503507659093060 0.0782761245708598
+0.5359650189609594 -0.7561145595512457 -0.3755426357749599
+0.3693313291199290 0.2223341552467790 0.9023092001865036
+0.4500961838491371 0.7138696318707015 0.5364732742432123
+0.9860612028101775 -0.0422449276805475 0.1609306384682893
+0.4619248785472341 0.3280166106625165 -0.8240330756156564
+-0.8949046135156228 -0.1766166131051672 0.4098198441799928
+0.1474944597771114 -0.3541254118214183 0.9234936800202637
+-0.9860612028101833 0.0422449276805802 -0.1609306384682460
+0.8949046135156530 0.1766166131050846 -0.4098198441799624
+0.1029852839925218 -0.6047145619938074 -0.7897558672106321
+-0.2504797437696364 0.9588399738152266 -0.1337378128096654
+-0.4500961838492435 -0.7138696318706331 -0.5364732742432141
+0.2504797437696580 -0.9588399738152272 0.1337378128096201
+-0.1029852839926024 0.6047145619938111 0.7897558672106186
+-0.6167298736902294 0.2645790829272963 0.7413785617182255
+-0.7642243334672891 0.6187044947487504 -0.1821151183020540
+0.6167298736902037 -0.2645790829273468 -0.7413785617182289
+-0.7005759276188255 0.2449703419446037 -0.6702110870528524
+0.7642243334673138 -0.6187044947487278 0.1821151183020281
+-0.1474944597770587 0.3541254118214171 -0.9234936800202724
+0.7005759276188496 -0.2449703419445155 0.6702110870528596
+0.1306802800484028 0.7953211078538522 -0.5919349624820214
+-0.3589395945546398 -0.9327311726563604 0.0342772084050255
+-0.8312562076671763 -0.5503507659093050 -0.0782761245708147
+-0.4619248785471858 -0.3280166106625792 0.8240330756156585
Added: branches/tools/points-mpas/points-mpas.cpp
===================================================================
--- branches/tools/points-mpas/points-mpas.cpp (rev 0)
+++ branches/tools/points-mpas/points-mpas.cpp 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,2039 @@
+//**************************************************
+// points-mpas.cpp
+//
+// Purpose:
+//
+// points-mpas.cpp is supposed to take in a triangulation defined on a sphere and a point set defined on a sphere, and create
+// a mpas grid file out of the two.
+//
+// Licensing:
+//
+// This code is distributed under the GNU LGPL license.
+//
+// Modified:
+//
+// 03 December 2010
+//
+// Author:
+//
+// Doug Jacobsen
+//
+//**************************************************
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <inttypes.h>
+#include <iostream>
+#include <fstream>
+#include <tr1/unordered_set>
+#include <vector>
+#include <utility>
+#include <math.h>
+#include <assert.h>
+
+//#include <netcdfcpp.h>
+#include "netcdfcpp.h"
+#include "triangulation.h"
+
+#define SEED 3729
+
+using namespace std;
+using namespace tr1;
+
+struct int_hasher {/*{{{*/
+ size_t operator()(const int v) const { return v; }
+};/*}}}*/
+
+int pt_type;
+int tri_base;
+int vert_levs;
+int num_tracers;
+double radius;
+double eps;
+int vordraw;
+
+/*{{{*/ // Grid information, points, triangles, ccenters, edges
+//unordered_set<pnt,pnt::hasher> edges;
+unordered_set<pnt,pnt::edge_hasher> edges;
+vector<pnt> edge_vec;
+vector<pnt> points;
+vector<pnt> ccenters;
+vector<tri> triangles;
+/*}}}*/
+/*{{{*/ // Real connectivity arrays
+vector<vector<int> > cellsOnCell, cellsOnEdge, cellsOnVertex;
+vector<vector<int> > edgesOnCell, edgesOnEdge, edgesOnVertex;
+vector<vector<int> > verticesOnCell, verticesOnEdge;
+vector<vector<double> > weightsOnEdge;
+/*}}}*/
+/*{{{*/ // Unique connectivity holders
+vector<unordered_set<int, int_hasher> > cellsOnCell_u;
+vector<unordered_set<int, int_hasher> > cellsOnEdge_u;
+vector<unordered_set<int, int_hasher> > cellsOnVertex_u;
+vector<unordered_set<int, int_hasher> > edgesOnCell_u;
+vector<unordered_set<int, int_hasher> > edgesOnEdge_u;
+vector<unordered_set<int, int_hasher> > edgesOnVertex_u;
+vector<unordered_set<int, int_hasher> > verticesOnCell_u;
+vector<unordered_set<int, int_hasher> > verticesOnEdge_u;
+/*}}}*/
+/*{{{*/ // Grid Parameters
+vector<vector<double> > kiteAreasOnVertex;
+vector<double> areaCell;
+vector<double> areaTriangle;
+vector<double> angleEdge;
+vector<double> dcEdge;
+vector<double> dvEdge;
+vector<double> fCell;
+vector<double> fEdge;
+vector<double> fVertex;
+/*}}}*/
+/*{{{*/ // Iterators for STL containers
+vector<vector<int> >::iterator vec_int_itr;
+vector<int>::iterator int_itr;
+vector<vector<double> >::iterator vec_dbl_itr;
+vector<double>::iterator dbl_itr;
+unordered_set<int, int_hasher>::iterator u_int_itr;
+vector<unordered_set<int, int_hasher> >::iterator us_itr;
+unordered_set<pnt,pnt::hasher>::iterator edge_itr;
+vector<pnt>::iterator point_itr;
+vector<tri>::iterator tri_itr;
+/*}}}*/
+
+/*{{{*/ // Function Declarations
+void readParameters();
+void readPoints();
+void readTriangulation();
+void triangulatePoints();
+void buildConnectivityArrays();
+void orderConnectivityArrays();
+void makeWeightsOnEdge();
+int outputGridDimensions(const string outputFilename);
+int outputGridAttributes(const string outputFilename);
+int outputGridCoordinates(const string outputFilename);
+int outputCellConnectivity(const string outputFilename);
+int outputEdgeConnectivity(const string outputFilename);
+int outputVertexConnectivity(const string outputFilename);
+int outputCellParameters(const string outputFilename);
+int outputVertexParameters(const string outputFilename);
+int outputEdgeParameters(const string outputFilename);
+int outputInitialConditions(const string outputFilename);
+int outputVordrawArrays(const string outputFilename);
+int writeGraphFile();
+double coriolisParameter(const pnt &p);
+/*}}}*/
+
+int main(){
+ int error_code;
+ string name = "grid.nc";
+
+ cout << "Reading in paramters" << endl;
+ readParameters();
+
+ cout << " --- Points are defined on a sphere --- Radius = " << radius << endl;
+
+ cout << "Reading in points" << endl;
+ readPoints();
+ cout << "Reading in triangles" << endl;
+ readTriangulation();
+
+ cout << "Building connectivity arrays" << endl;
+ buildConnectivityArrays();
+ cout << "Ordering connectivity arrays" << endl;
+ orderConnectivityArrays();
+ cout << "Making weights on edge" << endl;
+ makeWeightsOnEdge();
+
+ cout << "Writing grid dimensions" << endl;
+ if(error_code = outputGridDimensions(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing grid attributes" << endl;
+ if(error_code = outputGridAttributes(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing grid coordinates" << endl;
+ if(error_code = outputGridCoordinates(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing cell connectivity" << endl;
+ if(error_code = outputCellConnectivity(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing edge connectivity" << endl;
+ if(error_code = outputEdgeConnectivity(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing vertex connectivity" << endl;
+ if(error_code = outputVertexConnectivity(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing cell parameters" << endl;
+ if(error_code = outputCellParameters(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing edge parameters" << endl;
+ if(error_code = outputEdgeParameters(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Writing vertex parameters" << endl;
+ if(error_code = outputVertexParameters(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ cout << "Making and writing initial conditions" << endl;
+ if(error_code = outputInitialConditions(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ if(vordraw){
+ cout << "Writing arrays for Vordraw" << endl;
+ if(error_code = outputVordrawArrays(name)){
+ cout << "Error - " << error_code << endl;
+ exit(error_code);
+ }
+ }
+
+ cout << "Writing graph.info file" << endl;
+ writeGraphFile();
+ cout << points.size() << " cells." << endl;
+ cout << edge_vec.size() << " edges." << endl;
+ cout << ccenters.size() << " vertices." << endl;
+
+ return 0;
+}
+
+void readParameters(){/*{{{*/
+ ifstream params("Params");
+ if(!params){
+ cout << "Params file not found. Writing default, and using default values." << endl;
+ radius = 1.0;
+ vert_levs = 1;
+ num_tracers = 1;
+ eps = 0.0;
+ vordraw = 0;
+ params.close();
+
+ ofstream pout("Params");
+ pout << "Is the input Cartesian or Latitude-Longitude (0 - Cartesian, 1 - Lat-lon)" << endl << "0" << endl;
+ pout << "Are the triangles base zero or base one? (0 - base 0, 1 - base 1)" << endl << "0" << endl;
+ pout << "What is the radius of the sphere these points are defined on?" << endl << "1.0" << endl;
+ pout << "How many vertical levels do you want in the output grid?" << endl << "1" << endl;
+ pout << "How many tracers do you want in the output grid?" << endl << "1" << endl;
+ pout << "What was the convergence criteria used to make this grid?" << endl << "0.0" << endl;
+ pout << "Should this grid be vordraw compatible?" << endl << "0" << endl;
+ pout.close();
+ } else {
+ params.ignore(10000,'</font>
<font color="blue">');
+ params >> pt_type;
+ params.ignore(10000,'</font>
<font color="blue">');
+ params.ignore(10000,'</font>
<font color="blue">');
+ params >> tri_base;
+ params.ignore(10000,'</font>
<font color="blue">');
+ params.ignore(10000,'</font>
<font color="blue">');
+ params >> radius;
+ params.ignore(10000,'</font>
<font color="blue">');
+ params.ignore(10000,'</font>
<font color="blue">');
+ params >> vert_levs;
+ params.ignore(10000,'</font>
<font color="blue">');
+ params.ignore(10000,'</font>
<font color="blue">');
+ params >> num_tracers;
+ params.ignore(10000,'</font>
<font color="blue">');
+ params.ignore(10000,'</font>
<font color="blue">');
+ params >> eps;
+ params.ignore(10000,'</font>
<font color="blue">');
+ params.ignore(10000,'</font>
<font color="blue">');
+ params >> vordraw;
+ params.close();
+ }
+
+}/*}}}*/
+void readPoints(){/*{{{*/
+ /******************************************************************
+ *
+ * This function reads in the point set from a file, and inserts it
+ * into a vector.
+ *
+ ******************************************************************/
+ pnt p;
+ double lat, lon;
+ int i;
+ ifstream pt_start("SaveVertices");
+
+#ifdef _DEBUG
+ cout << "Setting up points." << endl;
+#endif
+
+ i = 0;
+ while(!pt_start.eof()){
+ if(pt_type){
+ pt_start >> lat >> lon;
+ p = pntFromLatLon(lat,lon);
+ } else {
+ pt_start >> p;
+ }
+ p.idx = i;
+ pt_start.ignore(10000,'</font>
<font color="blue">');
+
+ if(pt_start.good()){
+ points.push_back(p);
+ }
+ i++;
+ }
+ pt_start.close();
+}/*}}}*/
+void readTriangulation(){/*{{{*/
+ /*****************************************************************
+ *
+ * This function reads in the triangulation from a file.
+ * It computes all of the circumcenters, and edges and adds them into
+ * corresponding vectors and hash tables.
+ *
+ * A hash table is used to store the edges, to ensure only one copy
+ * of each edge is kept.
+ *
+ *****************************************************************/
+
+ pnt a, b, c, ccent, edge;
+ pair<unordered_set<pnt,pnt::hasher>::iterator,bool> out_pair;
+ double jv1, jv2, jv3;
+ int vi1, vi2, vi3;
+ int ei1, ei2, ei3;
+ int i, j, junk;
+ int min_vi;
+
+ ifstream tris("SaveTriangles");
+
+ tri t;
+
+ i = 0;
+ j = 0;
+ while(!tris.eof()){
+ tris >> vi1 >> vi2 >> vi3;
+ ei1 = -1;
+ ei2 = -1;
+ ei3 = -1;
+
+ if(tri_base == 1){
+ vi1--;
+ vi2--;
+ vi3--;
+ }
+
+ tris.ignore(1000,'</font>
<font color="gray">');
+
+ if(tris.good()){
+ a = points.at(vi1);
+ b = points.at(vi2);
+ c = points.at(vi3);
+
+ vi1 = a.idx;
+ vi2 = b.idx;
+ vi3 = c.idx;
+
+ if(!isCcw(a,b,c)){
+ junk = vi2;
+ vi2 = vi3;
+ vi3 = junk;
+
+ b = points.at(vi2);
+ c = points.at(vi3);
+ }
+
+ circumcenter(a,b,c,ccent);
+
+ ccent.normalize();
+
+ ccent.idx = i;
+
+ ccenters.push_back(ccent);
+
+ edge = (a+b)/2.0;
+ edge.idx = j;
+ edge.isBdry = 0;
+ edge.normalize();
+ if(b.idx > a.idx){
+ edge.vert_idx1 = a.idx;
+ edge.vert_idx2 = b.idx;
+ } else {
+ edge.vert_idx2 = a.idx;
+ edge.vert_idx1 = b.idx;
+ }
+
+ out_pair = edges.insert(edge);
+ if(out_pair.second){
+ edge_vec.push_back(edge);
+ j++;
+ }
+ ei1 = (*out_pair.first).idx;
+
+ edge = (b+c)/2.0;
+ edge.idx = j;
+ edge.isBdry = 0;
+ edge.normalize();
+ if(c.idx > b.idx){
+ edge.vert_idx1 = b.idx;
+ edge.vert_idx2 = c.idx;
+ } else {
+ edge.vert_idx2 = b.idx;
+ edge.vert_idx1 = c.idx;
+ }
+
+ out_pair = edges.insert(edge);
+ if(out_pair.second){
+ edge_vec.push_back(edge);
+ j++;
+ }
+ ei2 = (*out_pair.first).idx;
+
+ edge = (c+a)/2.0;
+ edge.idx = j;
+ edge.isBdry = 0;
+ edge.normalize();
+ if(a.idx > c.idx){
+ edge.vert_idx1 = c.idx;
+ edge.vert_idx2 = a.idx;
+ } else {
+ edge.vert_idx2 = c.idx;
+ edge.vert_idx1 = a.idx;
+ }
+
+ out_pair = edges.insert(edge);
+ if(out_pair.second){
+ edge_vec.push_back(edge);
+ j++;
+ }
+ ei3 = (*out_pair.first).idx;
+
+ t = tri(vi1,vi2,vi3,i);
+
+ t.ei1 = ei1;
+ t.ei2 = ei2;
+ t.ei3 = ei3;
+ triangles.push_back(t);
+ i++;
+ }
+ }
+
+ edges.clear();
+ tris.close();
+}/*}}}*/
+void buildConnectivityArrays(){/*{{{*/
+ /*************************************************************************
+ *
+ * This function takes the triangulation and point set previously read in
+ * from files, and computes the unique connectivity arrays for use in the
+ * ordering function.
+ *
+ * This is done by adding every item into a hash table, which takes care
+ * of duplicates for us. Later, we will order this data, and transfer it
+ * into a vector
+ *
+ *************************************************************************/
+ pnt a, b, c;
+ pnt edge1, edge2, edge3;
+ double area_temp;
+ double angle;
+ int vi1, vi2, vi3, ei1, ei2, ei3;
+
+ cellsOnCell.resize(points.size());
+ edgesOnCell.resize(points.size());
+ verticesOnCell.resize(points.size());
+ cellsOnCell_u.resize(points.size());
+ edgesOnCell_u.resize(points.size());
+ verticesOnCell_u.resize(points.size());
+ areaCell.resize(points.size());
+ fCell.resize(points.size());
+
+ cellsOnEdge.resize(edge_vec.size());
+ edgesOnEdge.resize(edge_vec.size());
+ verticesOnEdge.resize(edge_vec.size());
+ cellsOnEdge_u.resize(edge_vec.size());
+ edgesOnEdge_u.resize(edge_vec.size());
+ verticesOnEdge_u.resize(edge_vec.size());
+ angleEdge.resize(edge_vec.size());
+ dcEdge.resize(edge_vec.size());
+ dvEdge.resize(edge_vec.size());
+ fEdge.resize(edge_vec.size());
+
+ cellsOnVertex_u.resize(ccenters.size());
+ edgesOnVertex_u.resize(ccenters.size());
+ cellsOnVertex.resize(ccenters.size());
+ edgesOnVertex.resize(ccenters.size());
+ areaTriangle.resize(ccenters.size());
+ kiteAreasOnVertex.resize(ccenters.size());
+ fVertex.resize(ccenters.size());
+
+ for(int i = 0; i < points.size(); i ++){
+ areaCell[i] = 0.0;
+ }
+
+ for(int i = 0; i < ccenters.size(); i ++){
+ areaTriangle[i] = 0.0;
+ kiteAreasOnVertex[i].resize(3);
+ }
+
+ for(tri_itr = triangles.begin(); tri_itr != triangles.end(); ++tri_itr){
+ vi1 = (*tri_itr).vi1;
+ vi2 = (*tri_itr).vi2;
+ vi3 = (*tri_itr).vi3;
+
+ a = points.at(vi1);
+ b = points.at(vi2);
+ c = points.at(vi3);
+
+ ei1 = (*tri_itr).ei1;
+ ei2 = (*tri_itr).ei2;
+ ei3 = (*tri_itr).ei3;
+
+ edge1 = edge_vec.at(ei1);
+ edge2 = edge_vec.at(ei2);
+ edge3 = edge_vec.at(ei3);
+
+ cellsOnCell_u[vi1].insert(vi2);
+ cellsOnCell_u[vi1].insert(vi3);
+ cellsOnCell_u[vi2].insert(vi3);
+ cellsOnCell_u[vi2].insert(vi1);
+ cellsOnCell_u[vi3].insert(vi1);
+ cellsOnCell_u[vi3].insert(vi2);
+
+ cellsOnEdge_u[ei1].insert(vi1);
+ cellsOnEdge_u[ei1].insert(vi2);
+ cellsOnEdge_u[ei2].insert(vi2);
+ cellsOnEdge_u[ei2].insert(vi3);
+ cellsOnEdge_u[ei3].insert(vi3);
+ cellsOnEdge_u[ei3].insert(vi1);
+
+ cellsOnVertex_u[(*tri_itr).idx].insert(vi1);
+ cellsOnVertex_u[(*tri_itr).idx].insert(vi2);
+ cellsOnVertex_u[(*tri_itr).idx].insert(vi3);
+
+ edgesOnCell_u[vi1].insert(ei1);
+ edgesOnCell_u[vi1].insert(ei3);
+ edgesOnCell_u[vi2].insert(ei1);
+ edgesOnCell_u[vi2].insert(ei2);
+ edgesOnCell_u[vi3].insert(ei2);
+ edgesOnCell_u[vi3].insert(ei3);
+
+ edgesOnEdge_u[ei1].insert(ei2);
+ edgesOnEdge_u[ei1].insert(ei3);
+ edgesOnEdge_u[ei2].insert(ei1);
+ edgesOnEdge_u[ei2].insert(ei3);
+ edgesOnEdge_u[ei3].insert(ei1);
+ edgesOnEdge_u[ei3].insert(ei2);
+
+ edgesOnVertex_u[(*tri_itr).idx].insert(ei1);
+ edgesOnVertex_u[(*tri_itr).idx].insert(ei2);
+ edgesOnVertex_u[(*tri_itr).idx].insert(ei3);
+
+ verticesOnCell_u[vi1].insert((*tri_itr).idx);
+ verticesOnCell_u[vi2].insert((*tri_itr).idx);
+ verticesOnCell_u[vi3].insert((*tri_itr).idx);
+
+ verticesOnEdge_u[ei1].insert((*tri_itr).idx);
+ verticesOnEdge_u[ei2].insert((*tri_itr).idx);
+ verticesOnEdge_u[ei3].insert((*tri_itr).idx);
+
+ //areaCell
+ area_temp = triArea(points.at(vi1),ccenters.at((*tri_itr).idx),edge3);
+ area_temp += triArea(points.at(vi1),edge1,ccenters.at((*tri_itr).idx));
+ areaCell[vi1] += area_temp;
+ area_temp = triArea(points.at(vi2),ccenters.at((*tri_itr).idx),edge1);
+ area_temp += triArea(points.at(vi2),edge2,ccenters.at((*tri_itr).idx));
+ areaCell[vi2] += area_temp;
+ area_temp = triArea(points.at(vi3),ccenters.at((*tri_itr).idx),edge2);
+ area_temp += triArea(points.at(vi3),edge3,ccenters.at((*tri_itr).idx));
+ areaCell[vi3] += area_temp;
+
+ //Coriolis Parameters for cells and vertices
+ fCell[vi1] = coriolisParameter(points.at(vi1));
+ fCell[vi2] = coriolisParameter(points.at(vi2));
+ fCell[vi3] = coriolisParameter(points.at(vi3));
+ fVertex[(*tri_itr).idx] = coriolisParameter(ccenters.at((*tri_itr).idx));
+ }
+}/*}}}*/
+void orderConnectivityArrays(){/*{{{*/
+ /******************************************************************
+ *
+ * This function takes all of the hash tables that uniquely define the connectivity
+ * arrays, and orders them to be CCW and adjacent.
+ *
+ ******************************************************************/
+ int i, j, k;
+
+ // Order cellsOnEdge and verticesOnEdge, and make angleEdge
+ j = 0;
+ for(i = 0; i < edge_vec.size(); i++){/*{{{*/
+ int cell1, cell2, vert1, vert2;
+ pnt u, v, cross;
+ pnt np;
+ pnt edge;
+ double sign;
+
+ angleEdge[i] = 0;
+ // Ensure that u (cellsOnEdge2 - cellsOnEdge1) crossed with
+ // v (verticesOnEdge2 - verticesOnEdge1) is a right handed
+ // cross product
+ //
+ // The vectors u and v represent the perpendicular and parallel velocity
+ // directions along an edge
+ assert(cellsOnEdge_u[i].size() == 2);
+ u_int_itr = cellsOnEdge_u[i].begin();
+ cell1 = (*u_int_itr);
+ u_int_itr++;
+ cell2 = (*u_int_itr);
+
+ assert(verticesOnEdge_u[i].size() == 2);
+ u_int_itr = verticesOnEdge_u[i].begin();
+ vert1 = (*u_int_itr);
+ u_int_itr++;
+ vert2 = (*u_int_itr);
+
+ edge = edge_vec.at(i);
+ edge = gcIntersect(points.at(cell1),points.at(cell2),ccenters.at(vert1),ccenters.at(vert2));
+ edge.idx = i;
+ edge.isBdry = 0;
+ edge_vec.at(i) = edge;
+
+ u = points.at(cell2)-points.at(cell1);
+ v = ccenters.at(vert2)-ccenters.at(vert1);
+
+ dcEdge[i] = points.at(cell1).sphereDistance(points.at(cell2));
+ dvEdge[i] = ccenters.at(vert1).sphereDistance(ccenters.at(vert2));
+
+ fEdge[i] = coriolisParameter(edge);
+
+ cross = u.cross(v);
+ sign = cross.dot(edge);
+ cellsOnEdge[i].push_back(cell1);
+ cellsOnEdge[i].push_back(cell2);
+
+
+ if(sign < 0){
+ verticesOnEdge[i].push_back(vert2);
+ verticesOnEdge[i].push_back(vert1);
+
+ v = ccenters.at(vert1)-ccenters.at(vert2);
+ } else{
+ verticesOnEdge[i].push_back(vert1);
+ verticesOnEdge[i].push_back(vert2);
+ }
+
+ // angleEdge is either:
+ // 1. The angle the positive tangential direction (v)
+ // makes with the local northward direction.
+ // or
+ // 2. The angles the positive normal direction (u)
+ // makes with the local eastward direction.
+ np = pntFromLatLon(edge.getLat()+0.05, edge.getLon());
+ angleEdge[i] = acos((ccenters.at(vert2).getLat() - ccenters.at(vert1).getLat())/dvEdge[i]);
+
+ sign = planeAngle(edge, np, ccenters.at(vert2), edge);
+ sign = sign/fabs(sign);
+ angleEdge[i] = fabs(angleEdge[i]) * sign;
+ if (angleEdge[i] > M_PI) angleEdge[i] = angleEdge[i] - 2.0*M_PI;
+ if (angleEdge[i] < -M_PI) angleEdge[i] = angleEdge[i] + 2.0*M_PI;
+ }/*}}}*/
+
+ //Order cellsOnVertex and edgesOnVertex, areaTriangle
+ for(i = 0; i < ccenters.size(); i++){/*{{{*/
+
+ /*
+ * Since all of the *OnVertex arrays should only have 3 items, it's easy to order CCW
+ *
+ * Then using the triangles, build kitesAreasOnVertex, and areaTrianlge
+ */
+ pnt a, b, c;
+ pnt edge1, edge2, edge3;;
+ int c1, c2, c3;
+ int e1, e2, e3;
+ int swp_int;
+
+ u_int_itr = cellsOnVertex_u[i].begin();
+ c1 = (*u_int_itr);
+ u_int_itr++;
+ c2 = (*u_int_itr);
+ u_int_itr++;
+ c3 = (*u_int_itr);
+
+ a = points.at(c1);
+ b = points.at(c2);
+ c = points.at(c3);
+
+ if(!isCcw(a,b,c)){
+ swp_int = c2;
+ c2 = c3;
+ c3 = swp_int;
+ }
+
+ cellsOnVertex[i].clear();
+ cellsOnVertex[i].push_back(c1);
+ cellsOnVertex[i].push_back(c2);
+ cellsOnVertex[i].push_back(c3);
+
+ u_int_itr = edgesOnVertex_u[i].begin();
+ e1 = (*u_int_itr);
+ u_int_itr++;
+ e2 = (*u_int_itr);
+ u_int_itr++;
+ e3 = (*u_int_itr);
+
+ edge1 = edge_vec.at(e1);
+ edge2 = edge_vec.at(e2);
+ edge3 = edge_vec.at(e3);
+
+ if(!isCcw(edge1, edge2, edge3)){
+ swp_int = e2;
+ e2 = e3;
+ e3 = swp_int;
+ }
+
+ edgesOnVertex[i].clear();
+ edgesOnVertex[i].push_back(e1);
+ edgesOnVertex[i].push_back(e2);
+ edgesOnVertex[i].push_back(e3);
+
+ areaTriangle[i] = triArea(points.at(c1),points.at(c2),points.at(c3));
+ }/*}}}*/
+
+ //Order cellsOnCell, edgesOnCell and verticesOnCell
+ for(i = 0; i < points.size(); i++){/*{{{*/
+ // *
+ // * Since we have the *OnEdge arrays ordered correctly, and the *OnVertex arrays ordered correctly,
+ // * it should be easy to order the *OnCell arrays
+ // *
+ pnt cell_center;
+ int cur_edge;
+ int vert1, vert2;
+ int found;
+ size_t erased;
+
+ // Choose a starting edge on cell.
+ u_int_itr = edgesOnCell_u[i].begin();
+ cur_edge = (*u_int_itr);
+ cell_center = points.at(i);
+
+ while(!edgesOnCell_u[i].empty()){
+ //push this edge into edgesOnCell
+ edgesOnCell[i].push_back(cur_edge);
+ erased = edgesOnCell_u[i].erase(cur_edge);
+ if(erased != 1){
+ cout << " Edge " << cur_edge << " not valid" << endl;
+ cout << " On cell " << cell_center << endl;
+ cout << " Available edges on cell..." << endl;
+ for(u_int_itr = edgesOnCell_u[i].begin(); u_int_itr != edgesOnCell_u[i].end(); ++u_int_itr){
+ cout << (*u_int_itr) << " ";
+ }
+ cout << endl;
+ assert((int)erased == 1);
+ }
+
+ //Add the cell across the edge to cellsOnCell
+ if(cellsOnEdge[cur_edge].at(0) == i){
+ cellsOnCell[i].push_back(cellsOnEdge[cur_edge].at(1));
+ } else {
+ cellsOnCell[i].push_back(cellsOnEdge[cur_edge].at(0));
+ }
+
+ // Get the correct vertices on current edge, vert1 = starting vertex, vert2 = ending vertex
+ // at least in the ccw order here
+ vert1 = verticesOnEdge[cur_edge].at(0);
+ vert2 = verticesOnEdge[cur_edge].at(1);
+
+ if(!isCcw(cell_center, ccenters.at(vert1), ccenters.at(vert2))){
+ vert1 = verticesOnEdge[cur_edge].at(1);
+ vert2 = verticesOnEdge[cur_edge].at(0);
+ }
+
+ //Push the end vertex back into verticesOnCell (in the correct order).
+ verticesOnCell[i].push_back(vert2);
+
+ // Find the next edge by cycling over the edges connceted to the ending vertex.
+ found = 0;
+ for(int_itr = edgesOnVertex[vert2].begin(); int_itr != edgesOnVertex[vert2].end(); ++int_itr){
+ if((cellsOnEdge[(*int_itr)].at(0) == i || cellsOnEdge[(*int_itr)].at(1) == i)
+ && ((*int_itr) != cur_edge)){
+ cur_edge = (*int_itr);
+ found = 1;
+ break;
+ }
+ }
+
+ if(found != 1){
+ break;
+ }
+ }
+ }/*}}}*/
+
+ //Order edgesOnEdge
+ for(i = 0; i < edge_vec.size(); i++){/*{{{*/
+ /*
+ * Edges on edge should be easily built now that all of the other connectivity arrays are built
+ */
+ int cell1, cell2;
+ int found;
+
+ // Get cells connected to current edge
+ cell1 = cellsOnEdge[i].at(0);
+ cell2 = cellsOnEdge[i].at(1);
+
+ // Cell 1 loops
+ //Find current edge on cell 1, and add all edges from there to the end of edgesOnCell to edgesOnEdge
+ //Then, add all edges from the beginning of edgesOnCell to edgesOnEdge
+ //Since edgesOnCell should be CCW by now, these will all be CCW as well, and will iterate from cur_edge, around cell 1,
+ //and then back ccw around cell 2
+ found = 0;
+ for(int_itr = edgesOnCell.at(cell1).begin(); int_itr != edgesOnCell.at(cell1).end(); ++int_itr){
+ if((*int_itr) == i){
+ found = 1;
+ } else if (found && (*int_itr) != i){
+ edgesOnEdge[i].push_back((*int_itr));
+ }
+ }
+ assert(found == 1);
+ for(int_itr = edgesOnCell.at(cell1).begin(); int_itr != edgesOnCell.at(cell1).end(); ++int_itr){
+ if((*int_itr) == i){
+ found = 0;
+ } else if(found && (*int_itr) != i){
+ edgesOnEdge[i].push_back((*int_itr));
+ }
+ }
+ assert(found == 0);
+ //Cell 2 loops
+ //Do the same thing done in the cell1 loops, just in cell 2 now.
+ for(int_itr = edgesOnCell.at(cell2).begin(); int_itr != edgesOnCell.at(cell2).end(); ++int_itr){
+ if((*int_itr) == i){
+ found = 1;
+ } else if (found && (*int_itr) != i){
+ edgesOnEdge[i].push_back((*int_itr));
+ }
+ }
+ assert(found == 1);
+ for(int_itr = edgesOnCell.at(cell2).begin(); int_itr != edgesOnCell.at(cell2).end(); ++int_itr){
+ if((*int_itr) == i){
+ found = 0;
+ } else if(found && (*int_itr) != i){
+ edgesOnEdge[i].push_back((*int_itr));
+ }
+ }
+ assert(found == 0);
+ }/*}}}*/
+
+ cellsOnCell_u.clear();
+ cellsOnEdge_u.clear();
+ cellsOnVertex_u.clear();
+ edgesOnCell_u.clear();
+ edgesOnEdge_u.clear();
+ edgesOnVertex_u.clear();
+ verticesOnEdge_u.clear();
+ verticesOnCell_u.clear();
+}/*}}}*/
+void makeWeightsOnEdge(){/*{{{*/
+ /************************************************************
+ *
+ * This function computes weightsOnEdge based on kiteAreasOnVertex
+ * and areaCell.
+ *
+ * The weights correspond to edgesOnEdge and allow MPAS to reconstruct the
+ * edge perpendicular (previously defined as v) velocity using the edge
+ * neighbors
+ *
+ * Weight formulation is defined in J. Thurburn, et al. JCP 2009
+ * Numerical representation of geostrophic modes on arbitrarily
+ * structured C-grids
+ *
+ * I'm not entirely sure it's correct yet
+ ************************************************************/
+ int i, j, k;
+
+ weightsOnEdge.resize(edge_vec.size());
+
+ for(i = 0; i < edge_vec.size(); i++){
+ int jj;
+ int cur_edge;
+ int prev_edge;
+ int nei;
+ int cell1, cell2;
+ int edge1, edge2;
+ int vert1, vert2;
+ int neoc1, neoc2;
+ double de;
+ double area;
+ double sum_r;
+
+ cell1 = cellsOnEdge[i].at(0);
+ cell2 = cellsOnEdge[i].at(1);
+
+ neoc1 = edgesOnCell[cell1].size();
+ neoc2 = edgesOnCell[cell2].size();
+
+ weightsOnEdge[i].resize(edgesOnEdge[i].size());
+
+ sum_r = 0.0;
+ prev_edge = i;
+ de = dcEdge.at(i);
+ jj = 0;
+ for(j = 0; j < neoc1-1; j++){
+ cur_edge = edgesOnEdge[i].at(jj);
+
+ //Find the vertex that is shared between prev_edge and cur_edge
+ if(verticesOnEdge[prev_edge].at(0) == verticesOnEdge[cur_edge].at(0) ||
+ verticesOnEdge[prev_edge].at(0) == verticesOnEdge[cur_edge].at(1)){
+ vert1 = verticesOnEdge[prev_edge].at(0);
+ } else if (verticesOnEdge[prev_edge].at(1) == verticesOnEdge[cur_edge].at(0) ||
+ verticesOnEdge[prev_edge].at(1) == verticesOnEdge[cur_edge].at(1)){
+ vert1 = verticesOnEdge[prev_edge].at(1);
+ } else {
+ cout << "Edge " << prev_edge << " doesn't share a vertex with edge " << cur_edge << endl;
+ exit(1);
+ }
+
+ //Using the vertex, cell center, prev_edge and cur_edge compute the kite area using
+ //the two sub triangles, CC-PE-V and CC-V-CE
+ //
+ //This order is CCW
+ area = triArea(points.at(cell1),edge_vec.at(prev_edge), ccenters.at(vert1));
+ area += triArea(points.at(cell1), ccenters.at(vert1), edge_vec.at(cur_edge));
+
+ for(k = 0; k < 3; k++){
+ if(cellsOnVertex[vert1].at(k) == cell1){
+ kiteAreasOnVertex[vert1][k] = area;
+ }
+ }
+
+ //Compute running sum of area ratios for edges (using kites)
+ sum_r = sum_r + area/areaCell.at(cell1);
+
+ //Compute indicator function. -1 means inward edge normal, 1 mean outward edge normal.
+ //Inward means cell center is end 0 of current edge, where as
+ //Outward mean cell center is end 1 of current edge
+ if(cell1 == cellsOnEdge[cur_edge].at(0)){
+ nei = 1;
+ } else {
+ nei = -1;
+ }
+
+ //weightsOnEdge as defined in Thuburn paper referenced above.
+ //nei is indicator function, 0.5 is alpha (in equation 26)
+ //sum_r is running sum of area ratios
+ //dvEdge and de are the le and de terms used to scale the weights.
+ weightsOnEdge[i][jj] = (0.5 - sum_r)*nei*dvEdge[cur_edge]/de;
+
+ prev_edge = cur_edge;
+ jj++;
+ }
+
+ sum_r = 0.0;
+ prev_edge = i;
+ for(j = 0; j < neoc2-1; j++){
+ cur_edge = edgesOnEdge[i].at(jj);
+ if(verticesOnEdge[prev_edge].at(0) == verticesOnEdge[cur_edge].at(0) ||
+ verticesOnEdge[prev_edge].at(0) == verticesOnEdge[cur_edge].at(1)){
+ vert1 = verticesOnEdge[prev_edge].at(0);
+ } else if (verticesOnEdge[prev_edge].at(1) == verticesOnEdge[cur_edge].at(0) ||
+ verticesOnEdge[prev_edge].at(1) == verticesOnEdge[cur_edge].at(1)){
+ vert1 = verticesOnEdge[prev_edge].at(1);
+ } else {
+ cout << "Edge " << prev_edge << " doesn't share a vertex with edge " << cur_edge << endl;
+ cout << "Edge " << prev_edge << " has vertices " << verticesOnEdge[prev_edge].at(0) << " " << verticesOnEdge[prev_edge].at(1) << endl;
+ cout << "Edge " << cur_edge << " has vertices " << verticesOnEdge[cur_edge].at(0) << " " << verticesOnEdge[cur_edge].at(1) << endl;
+ exit(1);
+ }
+
+
+ area = triArea(points.at(cell2),edge_vec.at(prev_edge), ccenters.at(vert1));
+ area += triArea(points.at(cell2), ccenters.at(vert1), edge_vec.at(cur_edge));
+
+ for(k = 0; k < 3; k++){
+ if(cellsOnVertex[vert1].at(k) == cell2){
+ kiteAreasOnVertex[vert1][k] = area;
+ }
+ }
+
+ sum_r = sum_r + area/areaCell.at(cell2);
+
+ if(cell2 == cellsOnEdge[cur_edge].at(0)){
+ nei = -1;
+ } else {
+ nei = 1;
+ }
+
+ weightsOnEdge[i][jj] = (0.5 - sum_r)*nei*dvEdge[cur_edge]/de;
+ prev_edge = cur_edge;
+ jj++;
+ }
+ }
+}/*}}}*/
+
+int outputGridDimensions( const string outputFilename ){/*{{{*/
+ /************************************************************************
+ *
+ * This function writes the grid dimensions to the netcdf file named
+ * outputFilename
+ *
+ * **********************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Replace, NULL, 0, NcFile::Offset64Bits);
+
+ int nCells, maxEdges;
+ int junk;
+
+ nCells = points.size();
+
+ maxEdges = 0;
+
+ for(vec_int_itr = edgesOnCell.begin(); vec_int_itr != edgesOnCell.end(); ++vec_int_itr){
+ maxEdges = std::max(maxEdges, (int)(*vec_int_itr).size());
+ }
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // define dimensions
+ NcDim *nCellsDim;
+ NcDim *nEdgesDim;
+ NcDim *nVerticesDim;
+ NcDim *maxEdgesDim;
+ NcDim *maxEdges2Dim;
+ NcDim *TWODim;
+ NcDim *THREEDim;
+ NcDim *vertexDegreeDim;
+ NcDim *nVertLevelsDim;
+ NcDim *nTracersDim;
+ NcDim *timeDim;
+
+ // write dimensions
+ if (!(nCellsDim = grid.add_dim( "nCells", nCells) )) return NC_ERR;
+ if (!(nEdgesDim = grid.add_dim( "nEdges", (3*nCells)-6) )) return NC_ERR;
+ if (!(nVerticesDim = grid.add_dim( "nVertices", (2*nCells)-4) )) return NC_ERR;
+ if (!(maxEdgesDim = grid.add_dim( "maxEdges", maxEdges) )) return NC_ERR;
+ if (!(maxEdges2Dim = grid.add_dim( "maxEdges2", maxEdges*2) )) return NC_ERR;
+ if (!(TWODim = grid.add_dim( "TWO", 2) )) return NC_ERR;
+ if (!(vertexDegreeDim = grid.add_dim( "vertexDegree", 3) )) return NC_ERR;
+ if (!(nVertLevelsDim = grid.add_dim( "nVertLevels", vert_levs) )) return NC_ERR;
+ if (!(nTracersDim = grid.add_dim( "nTracers", num_tracers) )) return NC_ERR;
+ if (!(timeDim = grid.add_dim( "Time") )) return NC_ERR;
+
+ cout << " nCells --- " << nCells << endl;
+ cout << " nEdges --- " << (3*nCells)-6 << " " << edge_vec.size() << endl;
+ cout << " nVertices --- " << (2*nCells)-4 << " " << triangles.size() << endl;
+ cout << " maxEdges --- " << maxEdges << endl;
+ cout << " maxEdges2 --- " << maxEdges*2 << endl;
+ cout << " nVertLevels --- " << vert_levs << endl;
+ cout << " nTracers --- " << num_tracers << endl;
+
+ // file closed when file obj goes out of scope
+ return 0;
+}/*}}}*/
+int outputGridAttributes( const string outputFilename ){/*{{{*/
+ /************************************************************************
+ *
+ * This function writes the grid dimensions to the netcdf file named
+ * outputFilename
+ *
+ * **********************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+ NcBool sphereAtt, radiusAtt;
+ NcBool convAtt;
+
+ // write attributes
+ if (!(sphereAtt = grid.add_att( "on_a_sphere", "YES \0"))) return NC_ERR;
+ if (!(radiusAtt = grid.add_att( "sphere_radius", radius))) return NC_ERR;
+ if (!(convAtt = grid.add_att( "eps", eps ))) return NC_ERR;
+
+ // file closed when file obj goes out of scope
+ return 0;
+}/*}}}*/
+int outputGridCoordinates( const string outputFilename) {/*{{{*/
+ /************************************************************************
+ *
+ * This function writes the grid coordinates to the netcdf file named
+ * outputFilename
+ * This includes all cell centers, vertices, and edges.
+ * Both cartesian and lat,lon, as well as all of their indices
+ *
+ * **********************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nCellsDim = grid.get_dim( "nCells" );
+ NcDim *nEdgesDim = grid.get_dim( "nEdges" );
+ NcDim *nVerticesDim = grid.get_dim( "nVertices" );
+
+ int nCells = nCellsDim->size();
+ int nEdges = nEdgesDim->size();
+ int nVertices = nVerticesDim->size();
+
+ //Define nc variables
+ NcVar *xCellVar, *yCellVar, *zCellVar, *xEdgeVar, *yEdgeVar, *zEdgeVar, *xVertexVar, *yVertexVar, *zVertexVar;
+ NcVar *lonCellVar, *latCellVar, *lonEdgeVar, *latEdgeVar, *lonVertexVar, *latVertexVar;
+ NcVar *idx2cellVar, *idx2edgeVar, *idx2vertexVar;
+
+ int i;
+
+ double *x, *y, *z, *lat, *lon;
+ int *idxTo;
+
+ // Build and write cell coordinate arrays
+ x = new double[nCells];
+ y = new double[nCells];
+ z = new double[nCells];
+ lat = new double[nCells];
+ lon = new double[nCells];
+ idxTo = new int[nCells];
+ i = 0;
+ for(point_itr = points.begin(); point_itr != points.end(); ++point_itr){
+ x[i] = (*point_itr).x;
+ y[i] = (*point_itr).y;
+ z[i] = (*point_itr).z;
+ lat[i] = (*point_itr).getLat();
+ lon[i] = (*point_itr).getLon();
+ idxTo[i] = (*point_itr).idx+1;
+
+ i++;
+ }
+ if (!(latCellVar = grid.add_var("latCell", ncDouble, nCellsDim))) return NC_ERR;
+ if (!latCellVar->put(lat,nCells)) return NC_ERR;
+ if (!(lonCellVar = grid.add_var("lonCell", ncDouble, nCellsDim))) return NC_ERR;
+ if (!lonCellVar->put(lon,nCells)) return NC_ERR;
+ if (!(xCellVar = grid.add_var("xCell", ncDouble, nCellsDim))) return NC_ERR;
+ if (!xCellVar->put(x,nCells)) return NC_ERR;
+ if (!(yCellVar = grid.add_var("yCell", ncDouble, nCellsDim))) return NC_ERR;
+ if (!yCellVar->put(y,nCells)) return NC_ERR;
+ if (!(zCellVar = grid.add_var("zCell", ncDouble, nCellsDim))) return NC_ERR;
+ if (!zCellVar->put(z,nCells)) return NC_ERR;
+ if (!(idx2cellVar = grid.add_var("indexToCellID", ncInt, nCellsDim))) return NC_ERR;
+ if (!idx2cellVar->put(idxTo,nCells)) return NC_ERR;
+ free(x);
+ free(y);
+ free(z);
+ free(lat);
+ free(lon);
+ free(idxTo);
+
+ //Build and write edge coordinate arrays
+ x = new double[nEdges];
+ y = new double[nEdges];
+ z = new double[nEdges];
+ lat = new double[nEdges];
+ lon = new double[nEdges];
+ idxTo = new int[nEdges];
+
+ i = 0;
+ for(point_itr = edge_vec.begin(); point_itr != edge_vec.end(); ++point_itr){
+ x[i] = (*point_itr).x;
+ y[i] = (*point_itr).y;
+ z[i] = (*point_itr).z;
+ lat[i] = (*point_itr).getLat();
+ lon[i] = (*point_itr).getLon();
+ idxTo[i] = (*point_itr).idx+1;
+
+ i++;
+ }
+ if (!(latEdgeVar = grid.add_var("latEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!latEdgeVar->put(lat,nEdges)) return NC_ERR;
+ if (!(lonEdgeVar = grid.add_var("lonEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!lonEdgeVar->put(lon,nEdges)) return NC_ERR;
+ if (!(xEdgeVar = grid.add_var("xEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!xEdgeVar->put(x,nEdges)) return NC_ERR;
+ if (!(yEdgeVar = grid.add_var("yEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!yEdgeVar->put(y,nEdges)) return NC_ERR;
+ if (!(zEdgeVar = grid.add_var("zEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!zEdgeVar->put(z,nEdges)) return NC_ERR;
+ if (!(idx2edgeVar = grid.add_var("indexToEdgeID", ncInt, nEdgesDim))) return NC_ERR;
+ if (!idx2edgeVar->put(idxTo, nEdges)) return NC_ERR;
+ free(x);
+ free(y);
+ free(z);
+ free(lat);
+ free(lon);
+ free(idxTo);
+
+ //Build and write vertex coordinate arrays
+ x = new double[nVertices];
+ y = new double[nVertices];
+ z = new double[nVertices];
+ lat = new double[nVertices];
+ lon = new double[nVertices];
+ idxTo = new int[nVertices];
+
+ i = 0;
+ for(point_itr = ccenters.begin(); point_itr != ccenters.end(); ++point_itr){
+ x[i] = (*point_itr).x;
+ y[i] = (*point_itr).y;
+ z[i] = (*point_itr).z;
+ lat[i] = (*point_itr).getLat();
+ lon[i] = (*point_itr).getLon();
+ idxTo[i] = (*point_itr).idx+1;
+
+ i++;
+ }
+ if (!(latVertexVar = grid.add_var("latVertex", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!latVertexVar->put(lat,nVertices)) return NC_ERR;
+ if (!(lonVertexVar = grid.add_var("lonVertex", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!lonVertexVar->put(lon,nVertices)) return NC_ERR;
+ if (!(xVertexVar = grid.add_var("xVertex", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!xVertexVar->put(x,nVertices)) return NC_ERR;
+ if (!(yVertexVar = grid.add_var("yVertex", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!yVertexVar->put(y,nVertices)) return NC_ERR;
+ if (!(zVertexVar = grid.add_var("zVertex", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!zVertexVar->put(z,nVertices)) return NC_ERR;
+ if (!(idx2vertexVar = grid.add_var("indexToVertexID", ncInt, nVerticesDim))) return NC_ERR;
+ if (!idx2vertexVar->put(idxTo, nVertices)) return NC_ERR;
+ free(x);
+ free(y);
+ free(z);
+ free(lat);
+ free(lon);
+ free(idxTo);
+
+ return 0;
+}/*}}}*/
+int outputCellConnectivity( const string outputFilename) {/*{{{*/
+ /*****************************************************************
+ *
+ * This function writes all of the *OnCell arrays. Including
+ * cellsOnCell
+ * edgesOnCell
+ * verticesOnCell
+ * nEdgesonCell
+ *
+ * ***************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nCellsDim = grid.get_dim( "nCells" );
+ NcDim *maxEdgesDim = grid.get_dim( "maxEdges" );
+
+ int nCells = nCellsDim->size();
+ int maxEdges = maxEdgesDim->size();
+ int i, j;
+
+ // define nc variables
+ NcVar *cocVar, *nEocVar, *eocVar, *vocVar;
+
+ int *tmp_arr;
+
+ // Build and write COC array
+ tmp_arr = new int[nCells*maxEdges];
+
+ for(i = 0; i < nCells; i++){
+ for(j = 0; j < maxEdges; j++){
+ tmp_arr[i*maxEdges + j] = 0;
+ }
+ }
+
+ i = 0;
+ for(vec_int_itr = cellsOnCell.begin(); vec_int_itr != cellsOnCell.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*maxEdges + j] = (*int_itr) + 1;
+ j++;
+ }
+ i++;
+ }
+ if (!(cocVar = grid.add_var("cellsOnCell", ncInt, nCellsDim, maxEdgesDim))) return NC_ERR;
+ if (!cocVar->put(tmp_arr,nCells,maxEdges)) return NC_ERR;
+
+ // Build and write EOC array
+ for(i = 0; i < nCells; i++){
+ for(j = 0; j < maxEdges; j++){
+ tmp_arr[i*maxEdges + j] = 0;
+ }
+ }
+
+ i = 0;
+ for(vec_int_itr = edgesOnCell.begin(); vec_int_itr != edgesOnCell.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*maxEdges + j] = (*int_itr) + 1;
+ j++;
+ }
+
+ i++;
+ }
+
+ if (!(eocVar = grid.add_var("edgesOnCell", ncInt, nCellsDim, maxEdgesDim))) return NC_ERR;
+ if (!eocVar->put(tmp_arr,nCells,maxEdges)) return NC_ERR;
+
+ // Build and write VOC array
+ for(i = 0; i < nCells; i++){
+ for(j = 0; j < maxEdges; j++){
+ tmp_arr[i*maxEdges + j] = 0;
+ }
+ }
+
+ i = 0;
+ for(vec_int_itr = verticesOnCell.begin(); vec_int_itr != verticesOnCell.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*maxEdges + j] = (*int_itr) + 1;
+ j++;
+ }
+ i++;
+ }
+
+ if (!(vocVar = grid.add_var("verticesOnCell", ncInt, nCellsDim, maxEdgesDim))) return NC_ERR;
+ if (!vocVar->put(tmp_arr,nCells,maxEdges)) return NC_ERR;
+ free(tmp_arr);
+
+ //Build and write nEOC array
+ tmp_arr = new int[nCells];
+
+ i = 0;
+ for(vec_int_itr = edgesOnCell.begin(); vec_int_itr != edgesOnCell.end(); ++vec_int_itr){
+ tmp_arr[i] = (*vec_int_itr).size();
+ i++;
+ }
+
+ if (!(nEocVar = grid.add_var("nEdgesOnCell", ncInt, nCellsDim))) return NC_ERR;
+ if (!nEocVar->put(tmp_arr,nCells)) return NC_ERR;
+ verticesOnCell.clear();
+ edgesOnCell.clear();
+
+ free(tmp_arr);
+
+ return 0;
+}/*}}}*/
+int outputEdgeConnectivity( const string outputFilename) {/*{{{*/
+ /*****************************************************************
+ *
+ * This function writes all of the *OnEdge arrays. Including
+ * cellsOnEdge
+ * edgesOnEdge
+ * verticesOnEdge
+ * nEdgesOnEdge
+ *
+ * ***************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nEdgesDim = grid.get_dim( "nEdges" );
+ NcDim *maxEdges2Dim = grid.get_dim( "maxEdges2" );
+ NcDim *vertexDegreeDim = grid.get_dim( "vertexDegree" );
+ NcDim *twoDim = grid.get_dim( "TWO" );
+
+ // define nc variables
+ NcVar *coeVar, *nEoeVar, *eoeVar, *voeVar, *bdryEdgeVar;
+
+ int nEdges = nEdgesDim->size();
+ int maxEdges2 = maxEdges2Dim->size();
+ int vertexDegree = vertexDegreeDim->size();
+ int two = twoDim->size();
+ int i, j;
+
+ int *tmp_arr;
+
+ // Build and write EOE array
+ tmp_arr = new int[nEdges*maxEdges2];
+
+ for(i = 0; i < nEdges; i++){
+ for(j = 0; j < maxEdges2; j++){
+ tmp_arr[i*maxEdges2 + j] = 0;
+ }
+ }
+
+ i = 0;
+ for(vec_int_itr = edgesOnEdge.begin(); vec_int_itr != edgesOnEdge.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*maxEdges2 + j] = (*int_itr) + 1;
+ j++;
+ }
+
+ i++;
+ }
+
+ if (!(eoeVar = grid.add_var("edgesOnEdge", ncInt, nEdgesDim, maxEdges2Dim))) return NC_ERR;
+ if (!eoeVar->put(tmp_arr,nEdges,maxEdges2)) return NC_ERR;
+ free(tmp_arr);
+
+ // Build and write COE array
+ tmp_arr = new int[nEdges*two];
+ for(i = 0; i < nEdges; i++){
+ for(j = 0; j < two; j++){
+ tmp_arr[i*two + j] = 0;
+ }
+ }
+ i = 0;
+ for(vec_int_itr = cellsOnEdge.begin(); vec_int_itr != cellsOnEdge.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*two + j] = (*int_itr) + 1;
+ j++;
+ }
+
+ i++;
+ }
+
+ if (!(coeVar = grid.add_var("cellsOnEdge", ncInt, nEdgesDim, twoDim))) return NC_ERR;
+ if (!coeVar->put(tmp_arr,nEdges,two)) return NC_ERR;
+
+ // Build VOE array
+ i = 0;
+ for(vec_int_itr = verticesOnEdge.begin(); vec_int_itr != verticesOnEdge.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*two + j] = (*int_itr) + 1;
+ j++;
+ }
+
+ i++;
+ }
+
+ if (!(voeVar = grid.add_var("verticesOnEdge", ncInt, nEdgesDim, twoDim))) return NC_ERR;
+ if (!voeVar->put(tmp_arr,nEdges,two)) return NC_ERR;
+ free(tmp_arr);
+
+ // Build and write nEoe array
+ tmp_arr = new int[nEdges];
+ i = 0;
+ for(vec_int_itr = edgesOnEdge.begin(); vec_int_itr != edgesOnEdge.end(); ++vec_int_itr){
+ tmp_arr[i] = (*vec_int_itr).size();
+ i++;
+ }
+
+ if (!(nEoeVar = grid.add_var("nEdgesOnEdge", ncInt, nEdgesDim))) return NC_ERR;
+ if (!nEoeVar->put(tmp_arr,nEdges)) return NC_ERR;
+
+ // Build and write bdryEdge array
+ i = 0;
+ for(vec_int_itr = cellsOnEdge.begin(); vec_int_itr != cellsOnEdge.end(); ++vec_int_itr){
+ if((*vec_int_itr).size() != 2){
+ tmp_arr[i] = 1;
+ } else {
+ tmp_arr[i] = 0;
+ }
+ }
+
+ if (!(bdryEdgeVar = grid.add_var("boundaryEdge",ncInt, nEdgesDim))) return NC_ERR;
+ if (!bdryEdgeVar->put(tmp_arr,nEdges)) return NC_ERR;
+ free(tmp_arr);
+
+ cellsOnEdge.clear();
+// verticesOnEdge.clear(); // Needed for Initial conditions.
+ edgesOnEdge.clear();
+
+ return 0;
+}/*}}}*/
+int outputVertexConnectivity( const string outputFilename) {/*{{{*/
+ /*****************************************************************
+ *
+ * This function writes all of the *OnVertex arrays. Including
+ * cellsOnVertex
+ * edgesOnVertex
+ *
+ * ***************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nVerticesDim = grid.get_dim( "nVertices" );
+ NcDim *vertexDegreeDim = grid.get_dim( "vertexDegree" );
+
+ // define nc variables
+ NcVar *covVar, *eovVar, *bdryVertVar;
+
+ int nVertices = nVerticesDim->size();
+ int vertexDegree = vertexDegreeDim->size();
+ int i, j;
+
+ int *tmp_arr;
+
+ // Build and write COV array
+ tmp_arr = new int[nVertices*vertexDegree];
+
+ for(i = 0; i < nVertices; i++){
+ for(j = 0; j < vertexDegree; j++){
+ tmp_arr[i*vertexDegree + j] = 0;
+ }
+ }
+
+ i = 0;
+ for(vec_int_itr = cellsOnVertex.begin(); vec_int_itr != cellsOnVertex.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*vertexDegree + j] = (*int_itr) + 1;
+ j++;
+ }
+ i++;
+ }
+
+ if (!(covVar = grid.add_var("cellsOnVertex", ncInt, nVerticesDim, vertexDegreeDim))) return NC_ERR;
+ if (!covVar->put(tmp_arr,nVertices,vertexDegree)) return NC_ERR;
+
+ // Build and write EOV array
+ for(i = 0; i < nVertices; i++){
+ for(j = 0; j < vertexDegree; j++){
+ tmp_arr[i*vertexDegree + j] = 0;
+ }
+ }
+ i = 0;
+ for(vec_int_itr = edgesOnVertex.begin(); vec_int_itr != edgesOnVertex.end(); ++vec_int_itr){
+ j = 0;
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ tmp_arr[i*vertexDegree + j] = (*int_itr) + 1;
+ j++;
+ }
+
+ i++;
+ }
+ if (!(eovVar = grid.add_var("edgesOnVertex", ncInt, nVerticesDim, vertexDegreeDim))) return NC_ERR;
+ if (!eovVar->put(tmp_arr,nVertices,vertexDegree)) return NC_ERR;
+ free(tmp_arr);
+
+ // Build and write bdryVert array
+ tmp_arr = new int[nVertices];
+
+ i = 0;
+ for(vec_int_itr = cellsOnVertex.begin(); vec_int_itr != cellsOnVertex.end(); ++vec_int_itr){
+ if((*vec_int_itr).size() == vertexDegree){
+ tmp_arr[i] = 0;
+ } else {
+ tmp_arr[i] = 1;
+ }
+ i++;
+ }
+
+ if (!(bdryVertVar = grid.add_var("boundaryVertex", ncInt, nVerticesDim))) return NC_ERR;
+ if (!bdryVertVar->put(tmp_arr, nVertices)) return NC_ERR;
+
+ free(tmp_arr);
+
+ cellsOnVertex.clear();
+ edgesOnVertex.clear();
+
+ return 0;
+}/*}}}*/
+int outputCellParameters( const string outputFilename) {/*{{{*/
+ /*********************************************************
+ *
+ * This function writes all cell parameters, including
+ * areaCell
+ * fCell
+ *
+ * *******************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nCellsDim = grid.get_dim( "nCells" );
+
+ // define nc variables
+ NcVar *fCellVar;
+ NcVar *areacVar;
+
+ int nCells = nCellsDim->size();
+ int i, j;
+
+ if (!(fCellVar = grid.add_var("fCell", ncDouble, nCellsDim))) return NC_ERR;
+ if (!fCellVar->put(&fCell[0],nCells)) return NC_ERR;
+
+ if (!(areacVar = grid.add_var("areaCell", ncDouble, nCellsDim))) return NC_ERR;
+ if (!areacVar->put(&areaCell[0],nCells)) return NC_ERR;
+
+ fCell.clear();
+ areaCell.clear();
+
+ return 0;
+}/*}}}*/
+int outputVertexParameters( const string outputFilename) {/*{{{*/
+ /*********************************************************
+ *
+ * This function writes all vertex parameters, including
+ * areaTriangle
+ * kiteAreasOnVertex
+ * fVertex
+ *
+ * *******************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nVerticesDim = grid.get_dim( "nVertices" );
+ NcDim *vertexDegreeDim = grid.get_dim( "vertexDegree" );
+
+ // define nc variables
+ NcVar *fVertexVar;
+ NcVar *areatVar;
+ NcVar *kareaVar;
+
+ int nVertices = nVerticesDim->size();
+ int vertexDegree = vertexDegreeDim->size();
+ int i, j;
+
+ double *tmp_arr;
+
+ // Build and write fVertex
+ if (!(fVertexVar = grid.add_var("fVertex", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!fVertexVar->put(&fVertex[0],nVertices)) return NC_ERR;
+
+ // Build and write areaTriangle
+ if (!(areatVar = grid.add_var("areaTriangle", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!areatVar->put(&areaTriangle[0],nVertices)) return NC_ERR;
+
+ // Build and write kiteAreasOnVertex
+ tmp_arr = new double[nVertices*vertexDegree];
+ i = 0;
+ for(vec_dbl_itr = kiteAreasOnVertex.begin(); vec_dbl_itr != kiteAreasOnVertex.end(); ++vec_dbl_itr){
+ j = 0;
+ for(dbl_itr = (*vec_dbl_itr).begin(); dbl_itr != (*vec_dbl_itr).end(); ++dbl_itr){
+ tmp_arr[i*vertexDegree + j] = (*dbl_itr);
+ j++;
+ }
+ i++;
+ }
+
+ if (!(kareaVar = grid.add_var("kiteAreasOnVertex", ncDouble, nVerticesDim, vertexDegreeDim))) return NC_ERR;
+ if (!kareaVar->put(tmp_arr,nVertices,vertexDegree)) return NC_ERR;
+
+ free(tmp_arr);
+
+ fVertex.clear();
+ areaTriangle.clear();
+ kiteAreasOnVertex.clear();
+
+ return 0;
+}/*}}}*/
+int outputEdgeParameters( const string outputFilename) {/*{{{*/
+ /*********************************************************
+ *
+ * This function writes all grid parameters, including
+ * fEdge
+ * angleEdge
+ * dcEdge
+ * dvEdge
+ * weightsOnEdge
+ *
+ * *******************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nEdgesDim = grid.get_dim( "nEdges" );
+ NcDim *maxEdges2Dim = grid.get_dim( "maxEdges2" );
+
+ // define nc variables
+ NcVar *fEdgeVar;
+ NcVar *angleVar;
+ NcVar *kareaVar, *dcEdgeVar, *dvEdgeVar;
+ NcVar *woeVar;
+
+ int nEdges = nEdgesDim->size();
+ int maxEdges2 = maxEdges2Dim->size();
+ int i, j;
+
+ double *tmp_arr;
+
+ // Build and write fEdges
+ if (!(fEdgeVar = grid.add_var("fEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!fEdgeVar->put(&fEdge[0],nEdges)) return NC_ERR;
+
+ //Build and write angleEdge
+ if (!(angleVar = grid.add_var("angleEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!angleVar->put(&angleEdge[0],nEdges)) return NC_ERR;
+
+ //Build and write dcEdge
+ if (!(dcEdgeVar = grid.add_var("dcEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!dcEdgeVar->put(&dcEdge[0],nEdges)) return NC_ERR;
+
+ //Build and write dvEdge
+ if (!(dvEdgeVar = grid.add_var("dvEdge", ncDouble, nEdgesDim))) return NC_ERR;
+ if (!dvEdgeVar->put(&dvEdge[0],nEdges)) return NC_ERR;
+
+ //Build and write weightsOnEdge
+ tmp_arr = new double[nEdges*maxEdges2];
+
+ i = 0;
+ for(vec_dbl_itr = weightsOnEdge.begin(); vec_dbl_itr != weightsOnEdge.end(); ++vec_dbl_itr){
+ for(j = 0; j < maxEdges2; j++){
+ tmp_arr[i*maxEdges2 + j] = 0.0;
+ }
+
+ j = 0;
+ for(dbl_itr = (*vec_dbl_itr).begin(); dbl_itr != (*vec_dbl_itr).end(); ++dbl_itr){
+ tmp_arr[i*maxEdges2 + j] = (*dbl_itr);
+ j++;
+ }
+ i++;
+ }
+
+ if (!(woeVar = grid.add_var("weightsOnEdge", ncDouble, nEdgesDim, maxEdges2Dim))) return NC_ERR;
+ if (!woeVar->put(tmp_arr,nEdges,maxEdges2)) return NC_ERR;
+
+ free(tmp_arr);
+
+ fEdge.clear();
+ angleEdge.clear();
+ dcEdge.clear();
+ weightsOnEdge.clear();
+
+ return 0;
+}/*}}}*/
+int outputInitialConditions( const string outputFilename) {/*{{{*/
+ /***************************************************************************
+ *
+ * This function writes all initial conditions to the grid file, including
+ * (required)
+ * h_s
+ * u
+ * h
+ * tracers
+ *
+ * Any extra initial conditions can be written as needed, or added to tracers.
+ *
+ * *************************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ cout << "******************************************* " << endl;
+ cout << "* Using default initial conditions." << endl << "* For more fine grained control, edit the function outputInitialConditions." << endl;
+ cout << "******************************************* " << endl;
+
+ // fetch dimensions
+ NcDim *nCellsDim = grid.get_dim( "nCells" );
+ NcDim *nEdgesDim = grid.get_dim( "nEdges" );
+ NcDim *nVerticesDim = grid.get_dim( "nVertices" );
+ NcDim *nVertLevelsDim = grid.get_dim( "nVertLevels" );
+ NcDim *nTracersDim = grid.get_dim( "nTracers" );
+ NcDim *timeDim = grid.get_dim( "Time" );
+
+ // define nc variables
+ NcVar *uVar, *usrcVar, *hVar, *hsVar, *tracerVar;
+
+ int nCells = nCellsDim->size();
+ int nEdges = nEdgesDim->size();
+ int nVertices = nVerticesDim->size();
+ int nVertLevels = nVertLevelsDim->size();
+ int nTracers = nTracersDim->size();
+ int vert1, vert2;
+ int i, j, k;
+
+ double *tmp_arr;
+ double *tmp_arr2;
+ double h_s_val;
+
+ //Build and write h and h_s
+ tmp_arr = new double[nCells]; // h_s
+ tmp_arr2 = new double[nCells*nVertLevels]; // h
+
+ for(i = 0; i < nCells; i++){
+ h_s_val = 0.0;
+ for(j = 0; j < nVertLevels; j++){
+ tmp_arr2[i*nVertLevels + j] = 1000.0;
+ h_s_val -= tmp_arr2[i*nVertLevels + j];
+ }
+ tmp_arr[i] = h_s_val;
+ }
+
+ if (!(hsVar = grid.add_var("h_s", ncDouble, nCellsDim))) return NC_ERR;
+ if (!hsVar->put(tmp_arr,nCells)) return NC_ERR;
+ if (!(hVar = grid.add_var("h", ncDouble, timeDim, nCellsDim, nVertLevelsDim))) return NC_ERR;
+ if (!hVar->put(tmp_arr2,1,nCells,nVertLevels)) return NC_ERR;
+ free(tmp_arr);
+ free(tmp_arr2);
+
+ //Build and write u and u_src
+ tmp_arr = new double[nEdges*nVertLevels]; // u
+ tmp_arr2 = new double[nEdges*nVertLevels]; // u_src
+
+ for(i = 0; i < nEdges; i++){
+ for(j = 0; j < nVertLevels; j++){
+ vert1 = verticesOnEdge.at(i).at(0);
+ vert2 = verticesOnEdge.at(i).at(1);
+ if(j == 0){
+ // Setup solid body rotation
+ tmp_arr[i*nVertLevels + j] = (ccenters.at(vert2).z - ccenters.at(vert1).z) / dvEdge.at(i);
+ tmp_arr2[i*nVertLevels + j] = (ccenters.at(vert2).z - ccenters.at(vert1).z) / dvEdge.at(i);
+ } else {
+ tmp_arr2[i*nVertLevels + j] = 0.0;
+ }
+ }
+ }
+
+ if (!(uVar = grid.add_var("u", ncDouble, timeDim, nEdgesDim, nVertLevelsDim))) return NC_ERR;
+ if (!uVar->put(tmp_arr,1,nEdges,nVertLevels)) return NC_ERR;
+ if (!(usrcVar = grid.add_var("u_src", ncDouble, nEdgesDim, nVertLevelsDim))) return NC_ERR;
+ if (!usrcVar->put(tmp_arr2,nEdges,nVertLevels)) return NC_ERR;
+ free(tmp_arr);
+ free(tmp_arr2);
+
+ //Build and write tracers
+ tmp_arr = new double[nCells*nVertLevels*nTracers];
+ for(i = 0; i < nCells; i++){
+ for(j = 0; j < nVertLevels; j++){
+ for(k = 0; k < nTracers; k++){
+ tmp_arr[i*nVertLevels + j*nTracers + k] = 0.0;
+ }
+ }
+ }
+
+ if (!(tracerVar = grid.add_var("tracers", ncDouble, timeDim, nCellsDim, nVertLevelsDim, nTracersDim))) return NC_ERR;
+ if (!tracerVar->put(tmp_arr,1,nCells,nVertLevels,nTracers)) return NC_ERR;
+ free(tmp_arr);
+
+ return 0;
+}/*}}}*/
+int outputVordrawArrays( const string outputFilename) {/*{{{*/
+ /***************************************************************************
+ *
+ * This function writes all of the arrays for Vordraw compatibilty to the grid file, including
+ * cellProxy
+ * cellDensity
+ * vertexProxy
+ * vertexDensity
+ * nBorder
+ * xBorder
+ * yBorder
+ * zBorder
+ *
+ * *************************************************************************/
+ // Return this code to the OS in case of failure.
+ static const int NC_ERR = 2;
+
+ // set error behaviour (matches fortran behaviour)
+ NcError err(NcError::verbose_nonfatal);
+
+ // open the scvtmesh file
+ NcFile grid(outputFilename.c_str(), NcFile::Write);
+
+ // check to see if the file was opened
+ if(!grid.is_valid()) return NC_ERR;
+
+ // fetch dimensions
+ NcDim *nCellsDim = grid.get_dim( "nCells" );
+ NcDim *nVerticesDim = grid.get_dim( "nVertices" );
+
+ NcDim *nBorderDim;
+ if (!(nBorderDim = grid.add_dim( "nBorder", 1) )) return NC_ERR;
+
+ // define nc variables
+ NcVar *cDensVar, *cProxyVar, *vDensVar, *vProxyVar;
+ NcVar *xBordVar, *yBordVar, *zBordVar;
+
+ int nCells = nCellsDim->size();
+ int nVertices = nVerticesDim->size();
+ int nBorder = nBorderDim->size();
+ int i, j, k;
+ int junk_int;
+ double junk_dbl;
+
+ vector<double> dbl_tmp_arr;
+ vector<int> int_tmp_arr;
+ vector<double> xBorder, yBorder, zBorder;
+
+
+ //Build and write cellProxy and cellDens
+ int_tmp_arr.resize(nCells);
+ dbl_tmp_arr.resize(nCells);
+ ifstream cellproxy_in("SaveCellProxy");
+ ifstream celldens_in("SaveCellDensity");
+ if(!cellproxy_in){
+ for(i = 0 ; i < nCells; i++){
+ int_tmp_arr.at(i) = 1;
+ }
+ } else {
+ for(i = 0; i < nCells; i++){
+ cellproxy_in >> int_tmp_arr.at(i);
+
+ }
+ }
+
+ if(!celldens_in){
+ for(i = 0 ; i < nCells; i++){
+ dbl_tmp_arr.at(i) = 1;
+ }
+ } else {
+ for(i = 0; i < nCells; i++){
+ cellproxy_in >> dbl_tmp_arr.at(i);
+
+ }
+ }
+
+ cellproxy_in.close();
+ celldens_in.close();
+
+ if (!(cDensVar = grid.add_var("cellDensity", ncDouble, nCellsDim))) return NC_ERR;
+ if (!cDensVar->put(&dbl_tmp_arr.at(0),nCells)) return NC_ERR;
+ if (!(cProxyVar = grid.add_var("cellProxy", ncInt, nCellsDim))) return NC_ERR;
+ if (!cProxyVar->put(&int_tmp_arr.at(0),nCells)) return NC_ERR;
+
+ //Build and write vertProxy and vertDens
+ int_tmp_arr.clear();
+ dbl_tmp_arr.clear();
+ int_tmp_arr.resize(nVertices);
+ dbl_tmp_arr.resize(nVertices);
+ ifstream vertproxy_in("SaveVertexProxy");
+ ifstream vertdens_in("SaveVertexDensity");
+
+ if(!vertproxy_in){
+ for(i = 0 ; i < nVertices; i++){
+ int_tmp_arr.at(i) = 1;
+ }
+ } else {
+ for(i = 0; i < nVertices; i++){
+ vertproxy_in >> int_tmp_arr.at(i);
+
+ }
+ }
+
+ if(!vertdens_in){
+ for(i = 0 ; i < nVertices; i++){
+ dbl_tmp_arr.at(i) = 1;
+ }
+ } else {
+ for(i = 0; i < nVertices; i++){
+ vertproxy_in >> dbl_tmp_arr.at(i);
+
+ }
+ }
+
+ vertproxy_in.close();
+ vertdens_in.close();
+ if (!(vDensVar = grid.add_var("vertexDensity", ncDouble, nVerticesDim))) return NC_ERR;
+ if (!vDensVar->put(&dbl_tmp_arr.at(0),nVertices)) return NC_ERR;
+ if (!(vProxyVar = grid.add_var("vertexProxy", ncInt, nVerticesDim))) return NC_ERR;
+ if (!vProxyVar->put(&int_tmp_arr.at(0),nVertices)) return NC_ERR;
+
+ int_tmp_arr.clear();
+ dbl_tmp_arr.clear();
+ xBorder.push_back(0.0);
+ yBorder.push_back(0.0);
+ zBorder.push_back(0.0);
+
+ if (!(xBordVar = grid.add_var("xBorder", ncDouble, nBorderDim))) return NC_ERR;
+ if (!xBordVar->put(&xBorder.at(0),nBorder)) return NC_ERR;
+
+ if (!(yBordVar = grid.add_var("yBorder", ncDouble, nBorderDim))) return NC_ERR;
+ if (!yBordVar->put(&yBorder.at(0),nBorder)) return NC_ERR;
+
+ if (!(zBordVar = grid.add_var("zBorder", ncDouble, nBorderDim))) return NC_ERR;
+ if (!zBordVar->put(&zBorder.at(0),nBorder)) return NC_ERR;
+
+ return 0;
+}/*}}}*/
+
+int writeGraphFile() {/*{{{*/
+ //This function writes out the graph.info file, for use with metis domain decomposition software.
+ ofstream graph("graph.info");
+
+ graph << points.size() << " " << edge_vec.size() << endl;
+
+ for(vec_int_itr = cellsOnCell.begin(); vec_int_itr != cellsOnCell.end(); ++vec_int_itr){
+ for(int_itr = (*vec_int_itr).begin(); int_itr != (*vec_int_itr).end(); ++int_itr){
+ graph << (*int_itr)+1 << " ";
+ }
+ graph << endl;
+ }
+
+ graph.close();
+ cellsOnCell.clear();
+ return 0;
+}/*}}}*/
+
+double coriolisParameter(const pnt &p){/*{{{*/
+ /******************************************************
+ * This function returns the Coriolis parameter at a point.
+ ******************************************************/
+ return 2.0*7.292E-5*sin(p.getLat());
+}/*}}}*/
+
Added: branches/tools/points-mpas/triangulation.h
===================================================================
--- branches/tools/points-mpas/triangulation.h (rev 0)
+++ branches/tools/points-mpas/triangulation.h 2011-09-01 17:45:41 UTC (rev 972)
@@ -0,0 +1,533 @@
+#include <vector>
+#include <iostream>
+#include <fstream>
+
+class pnt {/*{{{*/
+ public:
+ double x, y, z;
+ int isBdry;
+ int idx;
+ int vert_idx1, vert_idx2;
+
+ pnt(double x_, double y_, double z_, int isBdry_, int idx_)
+ : x(x_), y(y_), z(z_), isBdry(isBdry_), idx(idx_) { }
+
+ pnt(double x_, double y_, double z_)
+ : x(x_), y(y_), z(z_), isBdry(0), idx(0) { }
+
+ pnt(double x_, double y_, double z_, int isBdry_)
+ : x(x_), y(y_), z(z_), isBdry(isBdry_), idx(0) { }
+
+ pnt()
+ : x(0.0), y(0.0), z(0.0), isBdry(0), idx(0) { }
+
+ friend pnt operator*(const double d, const pnt &p);
+ friend std::ostream & operator<<(std::ostream &os, const pnt &p);
+ friend std::istream & operator>>(std::istream &is, pnt &p);
+
+ pnt& operator=(const pnt &p){/*{{{*/
+ x = p.x;
+ y = p.y;
+ z = p.z;
+ isBdry = p.isBdry;
+ idx = p.idx;
+ return *this;
+ }/*}}}*/
+ bool operator==(const pnt &p) const {/*{{{*/
+ return (x == p.x) & (y == p.y) & (z == p.z);
+ }/*}}}*/
+ pnt operator-(const pnt &p) const {/*{{{*/
+ double x_, y_, z_;
+
+ x_ = x-p.x;
+ y_ = y-p.y;
+ z_ = z-p.z;
+
+ return pnt(x_,y_,z_,0,0);
+ }/*}}}*/
+ pnt operator+(const pnt &p) const {/*{{{*/
+ double x_, y_, z_;
+
+ x_ = x+p.x;
+ y_ = y+p.y;
+ z_ = z+p.z;
+
+ return pnt(x_,y_,z_,0,0);
+ }/*}}}*/
+ pnt operator*(double d) const {/*{{{*/
+ double x_, y_, z_;
+ x_ = x*d;
+ y_ = y*d;
+ z_ = z*d;
+ return pnt(x_,y_,z_,0,0);
+ }/*}}}*/
+ pnt operator/(double d) const {/*{{{*/
+ double x_, y_, z_;
+
+ if(d == 0.0){
+ std::cout << "pnt: operator/" << std::endl;
+ std::cout << (*this) << std::endl;
+ }
+
+ assert(d != 0.0);
+ x_ = x/d;
+ y_ = y/d;
+ z_ = z/d;
+ return pnt(x_,y_,z_,0,0);
+ }/*}}}*/
+ pnt& operator/=(double d){/*{{{*/
+ if(d == 0.0){
+ std::cout << "pnt: operator /=" << std::endl << (*this) << std::endl;
+ }
+ assert(d != 0.0);
+ x = x/d;
+ y = y/d;
+ z = z/d;
+ return *this;
+ }/*}}}*/
+ pnt& operator+=(const pnt &p){/*{{{*/
+ x += p.x;
+ y += p.y;
+ z += p.z;
+ return *this;
+ }/*}}}*/
+ double operator[](int i) const {/*{{{*/
+ if(i == 0){
+ return x;
+ } else if(i == 1){
+ return y;
+ } else {
+ return z;
+ }
+ }/*}}}*/
+ void normalize(){/*{{{*/
+ double norm;
+
+ norm = x*x + y*y + z*z;
+ if(norm == 0){
+ std::cout << "pnt: normalize" << std::endl;
+ std::cout << x << " " << y << " " << z << " " << isBdry << " " << idx << std::endl;
+
+ assert(norm != 0);
+ }
+ norm = sqrt(norm);
+
+ x = x/norm;
+ y = y/norm;
+ z = z/norm;
+ }/*}}}*/
+ double dot(const pnt &p) const {/*{{{*/
+ double junk;
+ junk = x*p.x+y*p.y+z*p.z;
+
+ return junk;
+ }/*}}}*/
+ double dotForAngle(const pnt &p) const {/*{{{*/
+ double junk;
+ junk = x*p.x+y*p.y+z*p.z;
+ if(junk > 1.0){
+ junk = 1.0;
+ }
+
+ if(junk < -1.0){
+ junk = -1.0;
+ }
+ return acos(junk);
+ }/*}}}*/
+ pnt cross(const pnt &p) const {/*{{{*/
+ double x_, y_, z_;
+
+ x_ = y*p.z - p.y*z;
+ y_ = z*p.x - p.z*x;
+ z_ = x*p.y - p.x*y;
+
+ return pnt(x_,y_,z_,0,0);
+ }/*}}}*/
+ double magnitude() const{/*{{{*/
+ return sqrt(x*x + y*y + z*z);
+ }/*}}}*/
+ double magnitude2() const {/*{{{*/
+ return x*x + y*y + z*z;
+ }/*}}}*/
+ void rotate(const pnt &vec, const double angle){/*{{{*/
+ double x_, y_, z_;
+ double cos_t, sin_t;
+
+ cos_t = cos(angle);
+ sin_t = sin(angle);
+
+ x_ = (cos_t + vec.x*vec.x *(1.0-cos_t)) * x
+ +(vec.x*vec.y*(1.0-cos_t) - vec.z * sin_t) * y
+ +(vec.x*vec.z*(1.0-cos_t) + vec.y * sin_t) * z;
+
+ y_ = (vec.y*vec.x*(1.0-cos_t) + vec.z*sin_t) * x
+ +(cos_t + vec.y*vec.x*(1.0 - cos_t)) * y
+ +(vec.y*vec.z*(1.0-cos_t) - vec.x*sin_t) * z;
+
+ z_ = (vec.z*vec.x*(1.0-cos_t)-vec.y*sin_t)*x
+ +(vec.z*vec.y*(1.0-cos_t)-vec.x*sin_t)*y
+ +(cos_t + vec.x*vec.x*(1.0-cos_t))*z;
+
+ x = x_;
+ y = y_;
+ z = z_;
+ }/*}}}*/
+ double getLat() const {/*{{{*/
+ return asin(z);
+ }/*}}}*/
+ double getLon() const {/*{{{*/
+ double lon;
+
+ lon = atan2(y,x);
+
+ if(lon < 0){
+ return 2.0 * M_PI + lon;
+ } else {
+ return lon;
+ }
+ }/*}}}*/
+ double sphereDistance(const pnt &p) const {/*{{{*/
+ pnt temp, cross;
+ double dot;
+ temp.x = x;
+ temp.y = y;
+ temp.z = z;
+
+ cross = temp.cross(p);
+ dot = temp.dot(p);
+
+ return atan2(cross.magnitude(),dot);
+ }/*}}}*/
+ struct hasher {/*{{{*/
+ size_t operator()(const pnt &p) const {
+ uint32_t hash;
+ size_t i, key[3] = { (size_t)p.x, (size_t)p.y, (size_t)p.z };
+ for(hash = i = 0; i < sizeof(key); ++i) {
+ hash += ((uint8_t *)key)[i];
+ hash += (hash << 10);
+ hash ^= (hash >> 6);
+ }
+ hash += (hash << 3);
+ hash ^= (hash >> 11);
+ hash += (hash << 15);
+ return hash;
+ }
+ };/*}}}*/
+ struct idx_hasher {/*{{{*/
+ size_t operator()(const pnt &p) const {
+ return (size_t)p.idx;
+ }
+ };/*}}}*/
+ struct edge_hasher {/*{{{*/
+ size_t operator()(const pnt &p) const {
+ uint32_t hash;
+ size_t i, key[2] = { (size_t)p.vert_idx1, (size_t)p.vert_idx2 };
+ for(hash = i = 0; i < sizeof(key); ++i) {
+ hash += ((uint8_t *)key)[i];
+ hash += (hash << 10);
+ hash ^= (hash >> 6);
+ }
+ hash += (hash << 3);
+ hash ^= (hash >> 11);
+ hash += (hash << 15);
+ return hash;
+ }
+ };/*}}}*/
+};/*}}}*/
+class tri {/*{{{*/
+ public:
+ int vi1, vi2, vi3;
+ int ei1, ei2, ei3;
+ int idx;
+
+ tri() : vi1(0), vi2(0), vi3(0), idx(0) { }
+
+ tri(int vi1_, int vi2_, int vi3_)
+ : vi1(vi1_), vi2(vi2_), vi3(vi3_), idx(0) { }
+
+ tri(int vi1_, int vi2_, int vi3_, int idx_)
+ : vi1(vi1_), vi2(vi2_), vi3(vi3_), idx(idx_) { }
+
+
+ friend std::ostream & operator<<(std::ostream &os, const tri &t);
+ friend std::istream & operator>>(std::istream &is, tri &t);
+
+ tri sortedTri(){/*{{{*/
+ int v1, v2, v3, swp_v;
+ v1 = vi1;
+ v2 = vi2;
+ v3 = vi3;
+
+ if(v1 > v2){
+ swp_v = v1;
+ v1 = v2;
+ v2 = swp_v;
+ }
+
+ if(v2 > v3){
+ swp_v = v2;
+ v2 = v3;
+ v3 = swp_v;
+ }
+
+ if(v1 > v2){
+ swp_v = v1;
+ v1 = v2;
+ v2 = swp_v;
+ }
+
+ return tri(v1,v2,v3);
+ }/*}}}*/
+ bool operator==(const tri &t) const {/*{{{*/
+ return (vi1 == t.vi1) & (vi2 == t.vi2) & (vi3 == t.vi3);
+ }/*}}}*/
+ struct hasher {/*{{{*/
+ size_t operator()(const tri &t) const {
+ uint32_t hash;
+ size_t i, key[3] = { t.vi1, t.vi2, t.vi3 };
+ for(hash = i = 0; i < sizeof(key); ++i) {
+ hash += ((uint8_t *)key)[i];
+ hash += (hash << 10);
+ hash ^= (hash >> 6);
+ }
+ hash += (hash << 3);
+ hash ^= (hash >> 11);
+ hash += (hash << 15);
+ return hash;
+ }
+ };/*}}}*/
+};/*}}}*/
+
+inline pnt operator*(const double d, const pnt &p){/*{{{*/
+ return pnt(d*p.x, d*p.y, d*p.z, 0, 0);
+}/*}}}*/
+
+inline std::ostream & operator<<(std::ostream &os, const pnt &p){/*{{{*/
+ os << '(' << p.x << ", " << p.y << ", " << p.z << ") bdry=" << p.isBdry << " idx=" << p.idx << ' ';
+ //os << p.x << " " << p.y << " " << p.z;
+ return os;
+}/*}}}*/
+inline std::istream & operator>>(std::istream &is, pnt &p){/*{{{*/
+ //is >> p.x >> p.y >> p.z >> p.isBdry >> p.idx;
+ is >> p.x >> p.y >> p.z;
+ return is;
+}/*}}}*/
+
+inline std::ostream & operator<<(std::ostream &os, const tri &t){/*{{{*/
+ //return os << t.vi1 << " " << t.vi2 << " " << t.vi3;
+ return os << t.vi1 << " " << t.vi2 << " " << t.vi3;
+}/*}}}*/
+inline std::istream & operator>>(std::istream &is, tri &t){/*{{{*/
+ //return is >> t.vi1 >> t.vi2 >> t.vi3;
+ return is >> t.vi1 >> t.vi2 >> t.vi3;
+}/*}}}*/
+
+void circumcenter(const pnt &A,const pnt &B,const pnt &C, pnt &cent){/*{{{*/
+ double a, b, c;
+ double pbc, apc, abp;
+
+ a = (B-C).magnitude2();
+ b = (C-A).magnitude2();
+ c = (A-B).magnitude2();
+
+ pbc = a*(-a + b + c);
+ apc = b*( a - b + c);
+ abp = c*( a + b - c);
+
+ cent = (pbc*A + apc*B + abp*C)/(pbc + apc + abp);
+}/*}}}*/
+double circumradius(const pnt &A, const pnt &B, const pnt &C){/*{{{*/
+
+ pnt ccenter;
+ pnt ac;
+
+ circumcenter(A,B,C,ccenter);
+ ac = A - ccenter;
+
+ return ac.magnitude();
+
+/*
+ pnt a(0.0,0.0,0.0,0);
+ pnt b(0.0,0.0,0.0,0);
+ pnt sub(0.0,0.0,0.0,0);
+ pnt cross(0.0,0.0,0.0,0);
+ double top, bot;
+
+ a = A-C;
+ b = B-C;
+
+ sub = a-b;
+ cross = a.cross(b);
+
+ top = a.magnitude() * b.magnitude() * sub.magnitude();
+ bot = 2.0 * cross.magnitude();
+
+ return top/bot; // */
+
+/* dx = a.x - b.x;
+ dy = a.y - b.y;
+ dz = a.z - b.z;
+
+ sa = sqrt(dx*dx + dy*dy + dz*dz);
+
+ dx = b.x - c.x;
+ dy = b.y - c.y;
+ dz = b.z - c.z;
+
+ sb = sqrt(dx*dx + dy*dy + dz*dz);
+
+ dx = c.x - a.x;
+ dy = c.y - a.y;
+ dz = c.z - a.z;
+
+ sc = sqrt(dx*dx + dy*dy + dz*dz);
+
+ dx = (sa+sb+sc)/2.0; // Semiperimeter
+ dy = 4.0*sqrt(dx*(sa+sb-dx)*(sa+sc-dx)*(sb+sc-dx)); // Bottom of circumradius computation
+
+ return (sa*sb*sc)/dy;*/
+}/*}}}*/
+double triAreaBAK(const pnt &A, const pnt &B, const pnt &C){/*{{{*/
+ /**************************************************************************
+ * - This function calculates the area of the triangle A,B,C on the
+ * surface of a sphere.
+ *
+ * Input: A, B, C
+ * A: vertex 1 of triangle
+ * B: vertex 2 of triangle
+ * C: vertex 3 of triangle
+ * Output: (returned value area)
+ * area: surface area of triangle on sphere.
+ **************************************************************************/
+ pnt u12, u23, u31;
+ double a, b, c, s, tanqe, area;
+ double sign;
+
+ area = 0.0;
+
+ //Compute Surface normal for triangle to get "sign" of area
+ u12 = B - A;
+ u23 = C - A;
+
+ u31 = u12.cross(u23);
+
+ u31.normalize();
+ sign = u31.magnitude();
+ assert(sign != 0.0);
+
+ //dot the surface norm with one of the points to get sign.
+ sign = u31.dot(A);
+ sign = sign/fabs(sign);
+
+ a = A.dotForAngle(B);
+ b = B.dotForAngle(C);
+ c = C.dotForAngle(A);
+
+ s = 0.5*(a+b+c);
+
+ tanqe = sqrt(tan(0.5*s)*tan(0.5*(s-a))*tan(0.5*(s-b))*tan(0.5*(s-c)));
+
+// area = sign*4.0*atan(tanqe);
+ area = 4.0*atan(tanqe);
+
+ if(isnan(area))
+ area = 0.0;
+
+ return area;
+}/*}}}*/
+double triArea(const pnt &A, const pnt &B, const pnt &C){/*{{{*/
+ double tanqe, s, a, b, c, e;
+
+ a = A.sphereDistance(B);
+ b = B.sphereDistance(C);
+ c = C.sphereDistance(A);
+ s = 0.5*(a+b+c);
+
+ tanqe = sqrt(tan(0.5*s)*tan(0.5*(s-a))*tan(0.5*(s-b))*tan(0.5*(s-c)));
+ e = 4.*atan(tanqe);
+
+ return e;
+}/*}}}*/
+int isCcw(const pnt &p1, const pnt &p2, const pnt &p3){/*{{{*/
+ // */
+ double sign;
+ pnt a, b, cross;
+
+ a = p2-p1;
+ b = p3-p1;
+ cross = a.cross(b);
+
+ sign = cross.dot(p1);
+
+ if(sign > 0){
+ return 1;
+ } else {
+ return 0;
+ }// */
+
+ /*
+ double side_check;
+ pnt a, cross;
+
+ a = p3 - p1;
+ a.normalize();
+ cross = p1.cross(a);
+ side_check = cross.dot(p2);
+
+ if(side_check < 0){
+ return 1;
+ } else {
+ return 0;
+ } // */
+}/*}}}*/
+pnt gcIntersect(const pnt &c1, const pnt &c2, const pnt &v1, const pnt &v2){/*{{{*/
+ pnt n, m,c ;
+ double dot;
+
+
+ //n = c1 cross c2
+ n = c1.cross(c2);
+
+ //m = v1 cross v2
+ m = v1.cross(v2);
+
+ //c = n cross m
+ c = n.cross(m);
+
+ dot = c1.dot(c);
+
+ dot = dot/fabs(dot);
+
+ c = c*dot;
+
+ c.normalize();
+
+ return c;
+}/*}}}*/
+
+double planeAngle(const pnt &A, const pnt &B, const pnt &C, const pnt &n){/*{{{*/
+ pnt ab;
+ pnt ac;
+ pnt cross;
+ double cos_angle;
+
+ ab = B - A;
+ ac = C - A;
+ cross = ab.cross(ac);
+ cos_angle = ab.dot(ac)/(ab.magnitude() * ac.magnitude());
+
+ cos_angle = std::min(std::max(cos_angle,-1.0),1.0);
+
+ if(cross.x*n.x + cross.y*n.y + cross.z*n.z >= 0){
+ return acos(cos_angle);
+ } else {
+ return -acos(cos_angle);
+ }
+}/*}}}*/
+pnt pntFromLatLon(const double &lat, const double &lon){/*{{{*/
+ pnt temp;
+ temp.x = cos(lon) * cos(lat);
+ temp.y = sin(lon) * cos(lat);
+ temp.z = sin(lat);
+ return temp;
+}/*}}}*/
</font>
</pre>