# include
# include
# include
# include
# include "square_grid.h"
/******************************************************************************/
double *square_grid ( int n, int ns[], double a[], double b[], int c[] )
/******************************************************************************/
/*
Purpose:
SQUARE_GRID: grid points over the interior of a square in 2D.
Discussion:
In 2D, a logically rectangular grid is to be created.
In the I-th dimension, the grid will use S(I) points.
The total number of grid points is
N = product ( 1 <= I <= 2 ) S(I)
Over the interval [A(i),B(i)], we have 5 choices for grid centering:
1: 0, 1/3, 2/3, 1
2: 1/5, 2/5, 3/5, 4/5
3: 0, 1/4, 2/4, 3/4
4: 1/4, 2/4, 3/4, 1
5: 1/8, 3/8, 5/8, 7/8
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
31 August 2014
Author:
John Burkardt
Parameters:
Input, int N, the number of points.
N = product ( 1 <= I <= 2 ) NS(I).
Input, int NS[2], the number of points along
each dimension.
Input, double A[2], B[2], the endpoints for each dimension.
Input, int C[2], the grid centering for each dimension.
1 <= C(*) <= 5.
Output, double SQUARE_GRID[2*N] = X(2*S(0)*S(1)), the points.
*/
{
int i;
int j;
static int m = 2;
int s;
double *x;
double *xs;
x = ( double * ) malloc ( m * n * sizeof ( double ) );
/*
Create the 1D grids in each dimension.
*/
for ( i = 0; i < m; i++ )
{
s = ns[i];
xs = ( double * ) malloc ( s * sizeof ( double ) );
for ( j = 0; j < s; j++ )
{
if ( c[i] == 1 )
{
if ( s == 1 )
{
xs[j] = 0.5 * ( a[i] + b[i] );
}
else
{
xs[j] = ( ( double ) ( s - j - 1 ) * a[i]
+ ( double ) ( j ) * b[i] )
/ ( double ) ( s - 1 );
}
}
else if ( c[i] == 2 )
{
xs[j] = ( ( double ) ( s - j ) * a[i]
+ ( double ) ( j + 1 ) * b[i] )
/ ( double ) ( s + 1 );
}
else if ( c[i] == 3 )
{
xs[j] = ( ( double ) ( s - j ) * a[i]
+ ( double ) ( j - 2 ) * b[i] )
/ ( double ) ( s );
}
else if ( c[i] == 4 )
{
xs[j] = ( ( double ) ( s - j - 1 ) * a[i]
+ ( double ) ( j + 1 ) * b[i] )
/ ( double ) ( s );
}
else if ( c[i] == 5 )
{
xs[j] = ( ( double ) ( 2 * s - 2 * j - 1 ) * a[i]
+ ( double ) ( 2 * j + 1 ) * b[i] )
/ ( double ) ( 2 * s );
}
}
r8vec_direct_product ( i, s, xs, m, n, x );
free ( xs );
}
return x;
}
/******************************************************************************/
void r8mat_transpose_print ( int m, int n, double a[], char *title )
/******************************************************************************/
/*
Purpose:
R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed.
Discussion:
An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
in column-major order.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
28 May 2008
Author:
John Burkardt
Parameters:
Input, int M, N, the number of rows and columns.
Input, double A[M*N], an M by N matrix to be printed.
Input, char *TITLE, a title.
*/
{
r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title );
return;
}
/******************************************************************************/
void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo,
int ihi, int jhi, char *title )
/******************************************************************************/
/*
Purpose:
R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
Discussion:
An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
in column-major order.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
10 September 2013
Author:
John Burkardt
Parameters:
Input, int M, N, the number of rows and columns.
Input, double A[M*N], an M by N matrix to be printed.
Input, int ILO, JLO, the first row and column to print.
Input, int IHI, JHI, the last row and column to print.
Input, char *TITLE, a title.
*/
{
# define INCX 5
int i;
int i2;
int i2hi;
int i2lo;
int i2lo_hi;
int i2lo_lo;
int inc;
int j;
int j2hi;
int j2lo;
fprintf ( stdout, "\n" );
fprintf ( stdout, "%s\n", title );
if ( m <= 0 || n <= 0 )
{
fprintf ( stdout, "\n" );
fprintf ( stdout, " (None)\n" );
return;
}
if ( ilo < 1 )
{
i2lo_lo = 1;
}
else
{
i2lo_lo = ilo;
}
if ( ihi < m )
{
i2lo_hi = m;
}
else
{
i2lo_hi = ihi;
}
for ( i2lo = i2lo_lo; i2lo <= i2lo_hi; i2lo = i2lo + INCX )
{
i2hi = i2lo + INCX - 1;
if ( m < i2hi )
{
i2hi = m;
}
if ( ihi < i2hi )
{
i2hi = ihi;
}
inc = i2hi + 1 - i2lo;
fprintf ( stdout, "\n" );
fprintf ( stdout, " Row:" );
for ( i = i2lo; i <= i2hi; i++ )
{
fprintf ( stdout, " %7d ", i - 1 );
}
fprintf ( stdout, "\n" );
fprintf ( stdout, " Col\n" );
fprintf ( stdout, "\n" );
if ( jlo < 1 )
{
j2lo = 1;
}
else
{
j2lo = jlo;
}
if ( n < jhi )
{
j2hi = n;
}
else
{
j2hi = jhi;
}
for ( j = j2lo; j <= j2hi; j++ )
{
fprintf ( stdout, "%5d:", j - 1 );
for ( i2 = 1; i2 <= inc; i2++ )
{
i = i2lo - 1 + i2;
fprintf ( stdout, " %14g", a[(i-1)+(j-1)*m] );
}
fprintf ( stdout, "\n" );
}
}
return;
# undef INCX
}
/******************************************************************************/
void r8vec_direct_product ( int factor_index, int factor_order,
double factor_value[], int factor_num, int point_num, double x[] )
/******************************************************************************/
/*
Purpose:
R8VEC_DIRECT_PRODUCT creates a direct product of R8VEC's.
Discussion:
An R8VEC is a vector of R8's.
To explain what is going on here, suppose we had to construct
a multidimensional quadrature rule as the product of K rules
for 1D quadrature.
The product rule will be represented as a list of points and weights.
The J-th item in the product rule will be associated with
item J1 of 1D rule 1,
item J2 of 1D rule 2,
...,
item JK of 1D rule K.
In particular,
X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK))
and
W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK)
So we can construct the quadrature rule if we can properly
distribute the information in the 1D quadrature rules.
This routine carries out that task.
Another way to do this would be to compute, one by one, the
set of all possible indices (J1,J2,...,JK), and then index
the appropriate information. An advantage of the method shown
here is that you can process the K-th set of information and
then discard it.
Example:
Rule 1:
Order = 4
X(1:4) = ( 1, 2, 3, 4 )
Rule 2:
Order = 3
X(1:3) = ( 10, 20, 30 )
Rule 3:
Order = 2
X(1:2) = ( 100, 200 )
Product Rule:
Order = 24
X(1:24) =
( 1, 10, 100 )
( 2, 10, 100 )
( 3, 10, 100 )
( 4, 10, 100 )
( 1, 20, 100 )
( 2, 20, 100 )
( 3, 20, 100 )
( 4, 20, 100 )
( 1, 30, 100 )
( 2, 30, 100 )
( 3, 30, 100 )
( 4, 30, 100 )
( 1, 10, 200 )
( 2, 10, 200 )
( 3, 10, 200 )
( 4, 10, 200 )
( 1, 20, 200 )
( 2, 20, 200 )
( 3, 20, 200 )
( 4, 20, 200 )
( 1, 30, 200 )
( 2, 30, 200 )
( 3, 30, 200 )
( 4, 30, 200 )
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
31 May 2009
Author:
John Burkardt
Parameters:
Input, int FACTOR_INDEX, the index of the factor being processed.
The first factor processed must be factor 0.
Input, int FACTOR_ORDER, the order of the factor.
Input, double FACTOR_VALUE[FACTOR_ORDER], the factor values
for factor FACTOR_INDEX.
Input, int FACTOR_NUM, the number of factors.
Input, int POINT_NUM, the number of elements in the direct product.
Input/output, double X[FACTOR_NUM*POINT_NUM], the elements of the
direct product, which are built up gradually.
Local Parameters:
Local, int START, the first location of a block of values to set.
Local, int CONTIG, the number of consecutive values to set.
Local, int SKIP, the distance from the current value of START
to the next location of a block of values to set.
Local, int REP, the number of blocks of values to set.
*/
{
static int contig = 0;
int i;
int j;
int k;
static int rep = 0;
static int skip = 0;
int start;
if ( factor_index == 0 )
{
contig = 1;
skip = 1;
rep = point_num;
for ( j = 0; j < point_num; j++ )
{
for ( i = 0; i < factor_num; i++ )
{
x[i+j*factor_num] = 0.0;
}
}
}
rep = rep / factor_order;
skip = skip * factor_order;
for ( i = 0; i < factor_order; i++ )
{
start = 0 + i * contig;
for ( k = 1; k <= rep; k++ )
{
for ( j = start; j < start + contig; j++ )
{
x[factor_index+j*factor_num] = factor_value[i];
}
start = start + skip;
}
}
contig = contig * factor_order;
return;
}
/******************************************************************************/
void timestamp ( )
/******************************************************************************/
/*
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:
24 September 2003
Author:
John Burkardt
Parameters:
None
*/
{
# define TIME_SIZE 40
static char time_buffer[TIME_SIZE];
const struct tm *tm;
time_t now;
now = time ( NULL );
tm = localtime ( &now );
strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
fprintf ( stdout, "%s\n", time_buffer );
return;
# undef TIME_SIZE
}