# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <string.h>
# include "cube_arbq_rule.h"
/******************************************************************************/
void cube_arbq ( int degree, int n, double x[], double w[] )
/******************************************************************************/
/*
Purpose:
CUBE_ARBQ returns a quadrature rule for the symmetric cube.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
08 July 2014
Author:
Original FORTRAN77 version by Hong Xiao, Zydrunas Gimbutas.
This C version by John Burkardt.
Reference:
Hong Xiao, Zydrunas Gimbutas,
A numerical algorithm for the construction of efficient quadrature
rules in two and higher dimensions,
Computers and Mathematics with Applications,
Volume 59, 2010, pages 663-676.
Parameters:
Input, int DEGREE, the degree of the quadrature rule.
1 <= DEGREE <= 15.
Input, int N, the number of nodes.
This can be determined by a call to CUBE_ARBQ_SIZE(DEGREE).
Output, double X[3*N], the coordinates of the nodes.
Output, double W[N], the weights.
*/
{
int i;
double w_sum;
if ( degree == 1 )
{
rule01 ( n, x, w );
}
else if ( degree == 2 )
{
rule02 ( n, x, w );
}
else if ( degree == 3 )
{
rule03 ( n, x, w );
}
else if ( degree == 4 )
{
rule04 ( n, x, w );
}
else if ( degree == 5 )
{
rule05 ( n, x, w );
}
else if ( degree == 6 )
{
rule06 ( n, x, w );
}
else if ( degree == 7 )
{
rule07 ( n, x, w );
}
else if ( degree == 8 )
{
rule08 ( n, x, w );
}
else if ( degree == 9 )
{
rule09 ( n, x, w );
}
else if ( degree == 10 )
{
rule10 ( n, x, w );
}
else if ( degree == 11 )
{
rule11 ( n, x, w );
}
else if ( degree == 12 )
{
rule12 ( n, x, w );
}
else if ( degree == 13 )
{
rule13 ( n, x, w );
}
else if ( degree == 14 )
{
rule14 ( n, x, w );
}
else if ( degree == 15 )
{
rule15 ( n, x, w );
}
else
{
fprintf ( stderr, "\n" );
fprintf ( stderr, "CUBE_ARBQ - Fatal error\n" );
fprintf ( stderr, " Illegal value of DEGREE.\n" );
exit ( 1 );
}
w_sum = r8vec_sum ( n, w );
for ( i = 0; i < n; i++ )
{
w[i] = 8.0 * w[i] / w_sum;
}
return;
}
/******************************************************************************/
void cube_arbq_gnuplot ( int n, double x[], char *header )
/******************************************************************************/
/*
Purpose:
CUBE_ARBQ_GNUPLOT: plot of a quadrature rule for the symmetric cube.
Licensing:
This code is distributed under the GNU GPL license.
Modified:
11 July 2014
Author:
John Burkardt
Reference:
Hong Xiao, Zydrunas Gimbutas,
A numerical algorithm for the construction of efficient quadrature
rules in two and higher dimensions,
Computers and Mathematics with Applications,
Volume 59, 2010, pages 663-676.
Parameters:
Input, double N, the number of nodes.
Input, double X[3*N], the coordinates of the nodes.
Input, char *HEADER, a string to be used to identify
the files created.
*/
{
char command_filename[255];
FILE *command_unit;
int j;
char node_filename[255];
FILE *node_unit;
char plot_filename[255];
char vertex_filename[255];
FILE *vertex_unit;
/*
Create the vertex file.
*/
strcpy ( vertex_filename, header );
strcat ( vertex_filename, "_vertices.txt" );
vertex_unit = fopen ( vertex_filename, "wt" );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, -1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, -1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, +1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, +1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, -1.0, -1.0 );
fprintf ( vertex_unit, "\n" );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, -1.0, +1.0 );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, -1.0, +1.0 );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, +1.0, +1.0 );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, +1.0, +1.0 );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, -1.0, +1.0 );
fprintf ( vertex_unit, "\n" );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, -1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, -1.0, +1.0 );
fprintf ( vertex_unit, "\n" );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, -1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, -1.0, +1.0 );
fprintf ( vertex_unit, "\n" );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, +1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", +1.0, +1.0, +1.0 );
fprintf ( vertex_unit, "\n" );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, +1.0, -1.0 );
fprintf ( vertex_unit, "%g %g %g\n", -1.0, +1.0, +1.0 );
fclose ( vertex_unit );
printf ( "\n" );
printf ( " Created vertex file '%s'\n", vertex_filename );
/*
Create node file.
*/
strcpy ( node_filename, header );
strcat ( node_filename, "_nodes.txt" );
node_unit = fopen ( node_filename, "wt" );
for ( j = 0; j < n; j++ )
{
fprintf ( node_unit, "%g %g %g\n", x[0+j*3], x[1+j*3], x[2+j*3] );
}
fclose ( node_unit );
printf ( " Created node file '%s'\n", node_filename );
/*
Create graphics command file.
*/
strcpy ( command_filename, header );
strcat ( command_filename, "_commands.txt" );
command_unit = fopen ( command_filename, "wt" );
fprintf ( command_unit, "# %s\n", command_filename );
fprintf ( command_unit, "#\n" );
fprintf ( command_unit, "# Usage:\n" );
fprintf ( command_unit, "# gnuplot < %s\n", command_filename );
fprintf ( command_unit, "#\n" );
fprintf ( command_unit, "set term png\n" );
strcpy ( plot_filename, header );
strcat ( plot_filename, ".png" );
fprintf ( command_unit, "set output '%s'\n", plot_filename );
fprintf ( command_unit, "set xlabel '<--- X --->'\n" );
fprintf ( command_unit, "set ylabel '<--- Y --->'\n" );
fprintf ( command_unit, "set zlabel '<--- Z --->'\n" );
fprintf ( command_unit, "set title '%s'\n", header );
fprintf ( command_unit, "set grid\n" );
fprintf ( command_unit, "set key off\n" );
fprintf ( command_unit, "set view equal xyz\n" );
fprintf ( command_unit, "set style data lines\n" );
fprintf ( command_unit, "set timestamp\n" );
fprintf ( command_unit, "splot '%s' with lines lw 3, \\\n", vertex_filename );
fprintf ( command_unit, " '%s' with points pt 7 lt 0\n", node_filename );
fclose ( command_unit );
printf ( " Created command file '%s'\n", command_filename );
return;
}
/******************************************************************************/
int cube_arbq_size ( int degree )
/******************************************************************************/
/*
Purpose:
CUBE_ARBQ_SIZE returns the size of quadrature rule for a cube.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
08 July 2014
Author:
John Burkardt.
Reference:
Hong Xiao, Zydrunas Gimbutas,
A numerical algorithm for the construction of efficient quadrature
rules in two and higher dimensions,
Computers and Mathematics with Applications,
Volume 59, 2010, pages 663-676.
Parameters:
Input, int DEGREE, the desired degree of exactness.
1 <= DEGREE <= 15.
Output, int CUBE_ARBQ_SIZE, the number of points in the
corresponding rule.
*/
{
int n;
const int n_save[15] = {
1, 4, 6, 10, 13,
22, 26, 42, 50, 73,
84, 116, 130, 172, 190 };
if ( degree < 1 || 15 < degree )
{
fprintf ( stderr, "\n" );
fprintf ( stderr, "CUBE_ARBQ_SIZE - Fatal error!\n" );
fprintf ( stderr, " Illegal value of DEGREE.\n" );
exit ( 1 );
}
n = n_save[degree-1];
return n;
}
/******************************************************************************/
double *lege3eva ( int degree, double z[] )