# include # include # include # include # include # include "fem2d_bvp_serene.h" /******************************************************************************/ double *basis_serene ( double xq, double yq, double xw, double ys, double xe, double yn, double xx[], double yy[] ) /******************************************************************************/ /* Purpose: BASIS_SERENE evaluates the serendipity basis functions. Discussion: This procedure assumes that a serendipity element has been defined, whose sides are parallel to coordinate axes. The local element numbering is YN 3--2--1 | | | | 4 8 | | | YS 5--6--7 | +--XW---XE-- We note that each basis function can be written as the product of three linear terms, which never result in an x^2y^2 term. Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, double XQ, YQ, the evaluation point. Input, double XW, YS, the coordinates of the lower left corner. Input, double XE, YN, the coordinates of the upper right corner. Input, double XX[8], YY[8], the coordinates of the 8 nodes. Output, double BASIS_SERENE[8], the value of the basis functions at (XQ,YQ). */ { double *v; v = ( double * ) malloc ( 8 * sizeof ( double ) ); v[0] = not1 ( xq, xw, xx[0] ) * not1 ( yq, ys, yy[0] ) * not2 ( xq, yq, xx[7], yy[7], xx[1], yy[1], xx[0], yy[0] ); v[1] = not1 ( xq, xw, xx[1] ) * not1 ( xq, xe, xx[1] ) * not1 ( yq, ys, yy[1] ); v[2] = not1 ( xq, xe, xx[2] ) * not1 ( yq, ys, yy[2] ) * not2 ( xq, yq, xx[1], yy[1], xx[3], yy[3], xx[2], yy[2] ); v[3] = not1 ( xq, xe, xx[3] ) * not1 ( yq, yn, yy[3] ) * not1 ( yq, ys, yy[3] ); v[4] = not1 ( xq, xe, xx[4] ) * not1 ( yq, yn, yy[4] ) * not2 ( xq, yq, xx[3], yy[3], xx[5], yy[5], xx[4], yy[4] ); v[5] = not1 ( xq, xe, xx[5] ) * not1 ( xq, xw, xx[5] ) * not1 ( yq, yn, yy[5] ); v[6] = not1 ( xq, xw, xx[6] ) * not1 ( yq, yn, yy[6] ) * not2 ( xq, yq, xx[5], yy[5], xx[7], yy[7], xx[6], yy[6] ); v[7] = not1 ( yq, ys, yy[7] ) * not1 ( yq, yn, yy[7] ) * not1 ( xq, xw, xx[7] ); return v; } /******************************************************************************/ double *basis_dx_serene ( double xq, double yq, double xw, double ys, double xe, double yn, double xx[], double yy[] ) /******************************************************************************/ /* Purpose: BASIS_DX_SERENE differentiates the serendipity basis functions for X. Discussion: This procedure assumes that a serendipity element has been defined, whose sides are parallel to coordinate axes. The local element numbering is YN 3--2--1 | | | | 4 8 | | | YS 5--6--7 | +--XW---XE-- We note that each basis function can be written as the product of three linear terms, which never result in an x^2y^2 term. Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, double XQ, YQ, the evaluation point. Input, double XW, YS, the coordinates of the lower left corner. Input, double XE, YN, the coordinates of the upper right corner. Input, double XX[8], YY[8], the coordinates of the 8 nodes. Output, double BASIS_DX_SERENE[8], the derivatives of the basis functions at (XQ,YQ) with respect to X. */ { double *vx; vx = ( double * ) malloc ( 8 * sizeof ( double ) ); vx[0] = not1d ( xw, xx[0] ) * not1 ( yq, ys, yy[0] ) * not2 ( xq, yq, xx[7], yy[7], xx[1], yy[1], xx[0], yy[0] ) + not1 ( xq, xw, xx[0] ) * not1 ( yq, ys, yy[0] ) * not2dx ( xx[7], yy[7], xx[1], yy[1], xx[0], yy[0] ); vx[1] = not1d ( xw, xx[1] ) * not1 ( xq, xe, xx[1] ) * not1 ( yq, ys, yy[1] ) + not1 ( xq, xw, xx[1] ) * not1d ( xe, xx[1] ) * not1 ( yq, ys, yy[1] ); vx[2] = not1d ( xe, xx[2] ) * not1 ( yq, ys, yy[2] ) * not2 ( xq, yq, xx[1], yy[1], xx[3], yy[3], xx[2], yy[2] ) + not1 ( xq, xe, xx[2] ) * not1 ( yq, ys, yy[2] ) * not2dx ( xx[1], yy[1], xx[3], yy[3], xx[2], yy[2] ); vx[3] = not1d ( xe, xx[3] ) * not1 ( yq, yn, yy[3] ) * not1 ( yq, ys, yy[3] ); vx[4] = not1d ( xe, xx[4] ) * not1 ( yq, yn, yy[4] ) * not2 ( xq, yq, xx[3], yy[3], xx[5], yy[5], xx[4], yy[4] ) + not1 ( xq, xe, xx[4] ) * not1 ( yq, yn, yy[4] ) * not2dx ( xx[3], yy[3], xx[5], yy[5], xx[4], yy[4] ); vx[5] = not1d ( xe, xx[5] ) * not1 ( xq, xw, xx[5] ) * not1 ( yq, yn, yy[5] ) + not1 ( xq, xe, xx[5] ) * not1d ( xw, xx[5] ) * not1 ( yq, yn, yy[5] ); vx[6] = not1d ( xw, xx[6] ) * not1 ( yq, yn, yy[6] ) * not2 ( xq, yq, xx[5], yy[5], xx[7], yy[7], xx[6], yy[6] ) + not1 ( xq, xw, xx[6] ) * not1 ( yq, yn, yy[6] ) * not2dx ( xx[5], yy[5], xx[7], yy[7], xx[6], yy[6] ); vx[7] = not1 ( yq, ys, yy[7] ) * not1 ( yq, yn, yy[7] ) * not1d ( xw, xx[7] ); return vx; } /******************************************************************************/ double *basis_dy_serene ( double xq, double yq, double xw, double ys, double xe, double yn, double xx[], double yy[] ) /******************************************************************************/ /* Purpose: BASIS_DY_SERENE differentiates the serendipity basis functions for Y. Discussion: This procedure assumes that a serendipity element has been defined, whose sides are parallel to coordinate axes. The local element numbering is YN 3--2--1 | | | | 4 8 | | | YS 5--6--7 | +--XW---XE-- We note that each basis function can be written as the product of three linear terms, which never result in an x^2y^2 term. Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, double XQ, YQ, the evaluation point. Input, double XW, YS, the coordinates of the lower left corner. Input, double XE, YN, the coordinates of the upper right corner. Input, double XX[8], YY[8], the coordinates of the 8 nodes. Output, double BASIS_DY_SERENE[8], the derivatives of the basis functions at (XQ,YQ) with respect to Y. */ { double *vy; vy = ( double * ) malloc ( 8 * sizeof ( double ) ); vy[0] = not1 ( xq, xw, xx[0] ) * not1d ( ys, yy[0] ) * not2 ( xq, yq, xx[7], yy[7], xx[1], yy[1], xx[0], yy[0] ) + not1 ( xq, xw, xx[0] ) * not1 ( yq, ys, yy[0] ) * not2dy ( xx[7], yy[7], xx[1], yy[1], xx[0], yy[0] ); vy[1] = not1 ( xq, xw, xx[1] ) * not1 ( xq, xe, xx[1] ) * not1d ( ys, yy[1] ); vy[2] = not1 ( xq, xe, xx[2] ) * not1d ( ys, yy[2] ) * not2 ( xq, yq, xx[1], yy[1], xx[3], yy[3], xx[2], yy[2] ) + not1 ( xq, xe, xx[2] ) * not1 ( yq, ys, yy[2] ) * not2dy ( xx[1], yy[1], xx[3], yy[3], xx[2], yy[2] ); vy[3] = not1 ( xq, xe, xx[3] ) * not1d ( yn, yy[3] ) * not1 ( yq, ys, yy[3] ) + not1 ( xq, xe, xx[3] ) * not1 ( yq, yn, yy[3] ) * not1d ( ys, yy[3] ); vy[4] = not1 ( xq, xe, xx[4] ) * not1d ( yn, yy[4] ) * not2 ( xq, yq, xx[3], yy[3], xx[5], yy[5], xx[4], yy[4] ) + not1 ( xq, xe, xx[4] ) * not1 ( yq, yn, yy[4] ) * not2dy ( xx[3], yy[3], xx[5], yy[5], xx[4], yy[4] ); vy[5] = not1 ( xq, xe, xx[5] ) * not1 ( xq, xw, xx[5] ) * not1d ( yn, yy[5] ); vy[6] = not1 ( xq, xw, xx[6] ) * not1d ( yn, yy[6] ) * not2 ( xq, yq, xx[5], yy[5], xx[7], yy[7], xx[6], yy[6] ) + not1 ( xq, xw, xx[6] ) * not1 ( yq, yn, yy[6] ) * not2dy ( xx[5], yy[5], xx[7], yy[7], xx[6], yy[6] ); vy[7] = not1d ( ys, yy[7] ) * not1 ( yq, yn, yy[7] ) * not1 ( xq, xw, xx[7] ) + not1 ( yq, ys, yy[7] ) * not1d ( yn, yy[7] ) * not1 ( xq, xw, xx[7] ); return vy; } /******************************************************************************/ double *fem2d_bvp_serene ( int nx, int ny, double a ( double x, double y ), double c ( double x, double y ), double f ( double x, double y ), double x[], double y[], int show11 ) /******************************************************************************/ /* Purpose: FEM2D_BVP_SERENE solves boundary value problem on a rectangle. Discussion: The program uses the finite element method, with piecewise serendipity basis functions to solve a 2D boundary value problem over a rectangle. The following differential equation is imposed inside the region: - d/dx a(x,y) du/dx - d/dy a(x,y) du/dy + c(x,y) * u(x,y) = f(x,y) where a(x,y), c(x,y), and f(x,y) are given functions. On the boundary, the solution is constrained to have the value 0. The finite element method will use a regular grid of NX nodes in X, and NY nodes in Y. Both NX and NY must be odd. The local element numbering is 3--2--1 | | 4 8 | | 5--6--7 The serendipity element mass matrix is a multiple of: 6.0, -6.0, 2.0, -8.0, 3.0, -8.0, 2.0, -6.0 -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0 2.0, -6.0, 6.0, -6.0, 2.0, -8.0, 3.0, -8.0 -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0 3.0, -8.0, 2.0, -6.0, 6.0, -6.0, 2.0, -8.0 -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0 2.0, -8.0, 3.0, -8.0, 2.0, -6.0, 6.0, -6.0 -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, int NX, NY, the number of X and Y grid values. NX and NY must be odd and at least 3. Input, function A(X,Y), evaluates a(x,y) Input, function C(X,Y), evaluates c(x,y) Input, function F(X,Y), evaluates f(x,y) Input, double X[NX], Y[NY], the mesh points. Input, int SHOW11, is 1 to print out the element matrix for the element in row 1, column 1. Output, double FEM2D_BVP_SERENE[MN], the finite element coefficients, which are also the value of the computed solution at the mesh points. */ { # define QUAD_NUM 3 double abscissa[QUAD_NUM] = { -0.774596669241483377035853079956, 0.000000000000000000000000000000, 0.774596669241483377035853079956 }; double *ae; double *amat; double aq; double *b; double *be; int cc; double cq; int e; int ex; int ex_num; int ey; int ey_num; double fq; int i; int ierror; int ii; int inc; int j; int jj; int k; int mm; int mn; int n; int node[8]; int quad_num = QUAD_NUM; int qx; int qy; int s; double scale; double *u; double *v; double *vx; double *vy; int w; double weight[QUAD_NUM] = { 0.555555555555555555555555555556, 0.888888888888888888888888888889, 0.555555555555555555555555555556 }; double wq; double xe; double xq; double xw; double xx[8]; double yn; double yq; double ys; double yy[8]; /* Make room for the matrix A and right hand side b. */ mn = fem2d_bvp_serene_node_num ( nx, ny ); amat = r8mat_zero_new ( mn, mn ); b = r8vec_zero_new ( mn ); /* Compute the matrix entries by integrating over each element. */ ex_num = ( nx - 1 ) / 2; ey_num = ( ny - 1 ) / 2; for ( ey = 0; ey < ey_num; ey++ ) { s = 2 * ey; mm = 2 * ey + 1; n = 2 * ey + 2; ys = y[s]; yn = y[n]; yy[0] = y[n]; yy[1] = y[n]; yy[2] = y[n]; yy[3] = y[mm]; yy[4] = y[s]; yy[5] = y[s]; yy[6] = y[s]; yy[7] = y[mm]; for ( ex = 0; ex < ex_num; ex++ ) { w = 2 * ex; cc = 2 * ex + 1; e = 2 * ex + 2; xe = x[e]; xw = x[w]; xx[0] = x[e]; xx[1] = x[cc]; xx[2] = x[w]; xx[3] = x[w]; xx[4] = x[w]; xx[5] = x[cc]; xx[6] = x[e]; xx[7] = x[e]; /* Node indices 3 2 1 wn cn en 4 8 wm em 5 6 7 ws cs es */ node[0] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 4; node[1] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 3; node[2] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 2; node[3] = ( 3 * ey + 2 ) * ey_num + 2 * ey + ex + 1; node[4] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex; node[5] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex + 1; node[6] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex + 2; node[7] = ( 3 * ey + 2 ) * ey_num + 2 * ey + ex + 2; if ( show11 ) { if ( ey == 0 && ex == 0 ) { ae = r8mat_zero_new ( 8, 8 ); be = r8vec_zero_new ( 8 ); } } for ( qx = 0; qx < quad_num; qx++ ) { xq = ( ( 1.0 - abscissa[qx] ) * x[e] + ( 1.0 + abscissa[qx] ) * x[w] ) / 2.0; for ( qy = 0; qy < quad_num; qy++ ) { yq = ( ( 1.0 - abscissa[qy] ) * y[n] + ( 1.0 + abscissa[qy] ) * y[s] ) / 2.0; wq = weight[qx] * ( x[e] - x[w] ) / 2.0 * weight[qy] * ( y[n] - y[s] ) / 2.0; v = basis_serene ( xq, yq, xw, ys, xe, yn, xx, yy ); vx = basis_dx_serene ( xq, yq, xw, ys, xe, yn, xx, yy ); vy = basis_dy_serene ( xq, yq, xw, ys, xe, yn, xx, yy ); aq = a ( xq, yq ); cq = c ( xq, yq ); fq = f ( xq, yq ); /* Build the element matrix. */ if ( show11 ) { if ( ey == 0 && ex == 0 ) { for ( i = 0; i < 8; i++ ) { for ( j = 0; j < 8; j++ ) { ae[i+j*8] = ae[i+j*8] + wq * ( vx[i] * aq * vx[j] + vy[i] * aq * vy[j] + v[i] * cq * v[j] ); } be[i] = be[i] + wq * ( v[i] * fq ); } } } for ( i = 0; i < 8; i++ ) { ii = node[i]; for ( j = 0; j < 8; j++ ) { jj = node[j]; amat[ii+jj*mn] = amat[ii+jj*mn] + wq * ( vx[i] * aq * vx[j] + vy[i] * aq * vy[j] + v[i] * cq * v[j] ); } b[ii] = b[ii] + wq * ( v[i] * fq ); } free ( v ); free ( vx ); free ( vy ); } } /* Print a sample element matrix. */ if ( show11 ) { if ( ey == 0 && ex == 0 ) { scale = 0.5 * ae[0+2*8]; for ( j = 0; j < 8; j++ ) { for ( i = 0; i < 8; i++ ) { ae[i+j*8] = ae[i+j*8] / scale; } } r8mat_print ( 8, 8, ae, " Wathen elementary mass matrix:" ); free ( ae ); free ( be ); } } } } /* Where a node is on the boundary, replace the finite element equation by a boundary condition. */ k = 0; for ( j = 0; j < ny; j++ ) { if ( ( j % 2 ) == 0 ) { inc = 1; } else { inc = 2; } for ( i = 0; i < nx; i = i + inc ) { if ( i == 0 || i == nx - 1 || j == 0 || j == ny - 1 ) { for ( jj = 0; jj < mn; jj++ ) { amat[k+jj*mn] = 0.0; } for ( ii = 0; ii < mn; ii++ ) { amat[ii+k*mn] = 0.0; } amat[k+k*mn] = 1.0; b[k] = 0.0; } k = k + 1; } } /* Solve the linear system. */ u = r8mat_solve2 ( mn, amat, b, &ierror ); free ( amat ); free ( b ); return u; # undef QUAD_NUM } /******************************************************************************/ int fem2d_bvp_serene_node_num ( int nx, int ny ) /******************************************************************************/ /* Purpose: FEM2D_BVP_SERENE_NODE_NUM counts the number of nodes. Discussion: The program uses the finite element method, with piecewise serendipity basis functions to solve a 2D boundary value problem over a rectangle. The grid uses NX nodes in the X direction and NY nodes in the Y direction. Both NX and NY must be odd. Because of the peculiar shape of the serendipity elements, counting the number of nodes and variables is a little tricky. Here is a grid for the case when NX = 7 and NY = 5, for which there are 29 nodes and variables. 23 24 25 26 27 28 29 19 20 21 22 12 13 14 15 16 17 18 8 9 10 11 1 2 3 4 5 6 7 Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, int NX, NY, the number of X and Y grid values. NX and NY must be odd and at least 3. Output, int FEM2D_BVP_SERENE_NODE_NUM, the number of nodes and variables. */ { int value; value = nx * ( ny + 1 ) / 2 + ( nx + 1 ) / 2 * ( ny - 1 ) / 2; return value; } /******************************************************************************/ double fem2d_h1s_error_serene ( int nx, int ny, double x[], double y[], double u[], double exact_ux ( double x, double y ), double exact_uy ( double x, double y ) ) /******************************************************************************/ /* Purpose: FEM2D_H1S_ERROR_SERENE: seminorm error of a finite element solution. Discussion: We assume the finite element method has been used, over a product region involving a grid of NX*NY nodes, with serendipity elements used for the basis. The finite element solution U(x,y) has been computed, and formulas for the exact derivatives Vx(x,y) and Vy(x,y) are known. This function estimates the H1 seminorm of the error: H1S = sqrt ( integral ( x, y ) ( Ux(x,y) - Vx(x,y) )^2 + ( Uy(x,y) - Vy(x,y) )^2 dx dy ) Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, int NX, NY, the number of nodes. Input, double X[NX], Y[NY], the grid coordinates. Input, double U[*], the finite element coefficients. Input, function EQX = EXACT_UX(X,Y), function EQY = EXACT_UY(X,Y) returns the exact derivatives with respect to X and Y. Output, double FEM2D_BVP_H1S, the estimated seminorm of the error. */ { # define QUAD_NUM 3 double abscissa[QUAD_NUM] = { -0.774596669241483377035853079956, 0.000000000000000000000000000000, 0.774596669241483377035853079956 }; int cc; int e; int ex; int ex_num; int ey; int ey_num; double exq; double eyq; double h1s; int k; int mm; int n; int node[8]; int quad_num = QUAD_NUM; int qx; int qy; int s; double uxq; double uyq; double *vx; double *vy; int w; double weight[QUAD_NUM] = { 0.555555555555555555555555555556, 0.888888888888888888888888888889, 0.555555555555555555555555555556 }; double wq; double xe; double xq; double xw; double xx[8]; double yn; double yq; double ys; double yy[8]; h1s = 0.0; /* Quadrature definitions. */ ex_num = ( nx - 1 ) / 2; ey_num = ( ny - 1 ) / 2; for ( ey = 0; ey < ey_num; ey++ ) { s = 2 * ey; mm = 2 * ey + 1; n = 2 * ey + 2; ys = y[s]; yn = y[n]; yy[0] = y[n]; yy[1] = y[n]; yy[2] = y[n]; yy[3] = y[mm]; yy[4] = y[s]; yy[5] = y[s]; yy[6] = y[s]; yy[7] = y[mm]; for ( ex = 0; ex < ex_num; ex++ ) { w = 2 * ex; cc = 2 * ex + 1; e = 2 * ex + 2; xe = x[e]; xw = x[w]; xx[0] = x[e]; xx[1] = x[cc]; xx[2] = x[w]; xx[3] = x[w]; xx[4] = x[w]; xx[5] = x[cc]; xx[6] = x[e]; xx[7] = x[e]; /* Node indices 3 2 1 wn cn en 4 8 wm em 5 6 7 ws cs es */ node[0] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 4; node[1] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 3; node[2] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 2; node[3] = ( 3 * ey + 2 ) * ey_num + 2 * ey + ex + 1; node[4] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex; node[5] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex + 1; node[6] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex + 2; node[7] = ( 3 * ey + 2 ) * ey_num + 2 * ey + ex + 2; for ( qx = 0; qx < quad_num; qx++ ) { xq = ( ( 1.0 - abscissa[qx] ) * x[e] + ( 1.0 + abscissa[qx] ) * x[w] ) / 2.0; for ( qy = 0; qy < quad_num; qy++ ) { yq = ( ( 1.0 - abscissa[qy] ) * y[n] + ( 1.0 + abscissa[qy] ) * y[s] ) / 2.0; wq = weight[qx] * ( x[e] - x[w] ) / 2.0 * weight[qy] * ( y[n] - y[s] ) / 2.0; vx = basis_dx_serene ( xq, yq, xw, ys, xe, yn, xx, yy ); vy = basis_dy_serene ( xq, yq, xw, ys, xe, yn, xx, yy ); uxq = 0.0; uyq = 0.0; for ( k = 0; k < 8; k++ ) { uxq = uxq + u[node[k]] * vx[k]; uyq = uyq + u[node[k]] * vy[k]; } exq = exact_ux ( xq, yq ); eyq = exact_uy ( xq, yq ); h1s = h1s + wq * ( pow ( uxq - exq, 2 ) + pow ( uyq - eyq, 2 ) ); free ( vx ); free ( vy ); } } } } h1s = sqrt ( h1s ); return h1s; # undef QUAD_NUM } /******************************************************************************/ double fem2d_l1_error_serene ( int nx, int ny, double x[], double y[], double u[], double exact ( double x, double y ) ) /******************************************************************************/ /* Purpose: FEM2D_L1_ERROR_SERENE: l1 error norm of a finite element solution. Discussion: We assume the finite element method has been used, over a product region with NX*NY nodes and the serendipity element. The coefficients U(*) have been computed, and a formula for the exact solution is known. Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, int NX, NY, the number of X and Y grid values. Input, double X[NX], Y[NY], the grid coordinates. Input, double U[*], the finite element coefficients. Input, function EQ = EXACT(X,Y), returns the value of the exact solution at the point (X,Y). Output, double FEM2D_L1_ERROR_SERENE, the little l1 norm of the error. */ { int i; int inc; int j; int k; double e1; e1 = 0.0; k = 0; for ( j = 0; j < ny; j++ ) { if ( ( j % 2 ) == 0 ) { inc = 1; } else { inc = 2; } for ( i = 0; i < nx; i = i + inc ) { e1 = e1 + fabs ( u[k] - exact ( x[i], y[j] ) ); k = k + 1; } } e1 = e1 / ( double ) ( k ); return e1; } /******************************************************************************/ double fem2d_l2_error_serene ( int nx, int ny, double x[], double y[], double u[], double exact ( double x, double y ) ) /******************************************************************************/ /* Purpose: FEM2D_L2_ERROR_SERENE: L2 error norm of a finite element solution. Discussion: The finite element method has been used, over a rectangle, involving a grid of NXxNY nodes, with serendipity elements used for the basis. The finite element coefficients have been computed, and a formula for the exact solution is known. This function estimates E2, the L2 norm of the error: E2 = Integral ( X, Y ) ( U(X,Y) - EXACT(X,Y) )^2 dX dY Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, int NX, NY, the number of nodes in the X and Y directions. Input, double X[NX], Y[NY], the grid coordinates. Input, double U[*], the finite element coefficients. Input, function EQ = EXACT(X,Y), returns the value of the exact solution at the point (X,Y). Output, double E2, the estimated L2 norm of the error. */ { # define QUAD_NUM 3 double abscissa[QUAD_NUM] = { -0.774596669241483377035853079956, 0.000000000000000000000000000000, 0.774596669241483377035853079956 }; int cc; int e; double e2; double eq; int ex; int ex_num; int ey; int ey_num; int k; int mm; int n; int node[8]; int quad_num = QUAD_NUM; int qx; int qy; int s; double uq; double *v; int w; double weight[QUAD_NUM] = { 0.555555555555555555555555555556, 0.888888888888888888888888888889, 0.555555555555555555555555555556 }; double wq; double xe; double xq; double xw; double xx[8]; double yn; double yq; double ys; double yy[8]; e2 = 0.0; /* Compute the matrix entries by integrating over each element. */ ex_num = ( nx - 1 ) / 2; ey_num = ( ny - 1 ) / 2; for ( ey = 0; ey < ey_num; ey++ ) { s = 2 * ey; mm = 2 * ey + 1; n = 2 * ey + 2; ys = y[s]; yn = y[n]; yy[0] = y[n]; yy[1] = y[n]; yy[2] = y[n]; yy[3] = y[mm]; yy[4] = y[s]; yy[5] = y[s]; yy[6] = y[s]; yy[7] = y[mm]; for ( ex = 0; ex < ex_num; ex++ ) { w = 2 * ex; cc = 2 * ex + 1; e = 2 * ex + 2; xe = x[e]; xw = x[w]; xx[0] = x[e]; xx[1] = x[cc]; xx[2] = x[w]; xx[3] = x[w]; xx[4] = x[w]; xx[5] = x[cc]; xx[6] = x[e]; xx[7] = x[e]; /* Node indices 3 2 1 wn cn en 4 8 wm em 5 6 7 ws cs es */ node[0] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 4; node[1] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 3; node[2] = ( 3 * ey + 3 ) * ey_num + 2 * ey + 2 * ex + 2; node[3] = ( 3 * ey + 2 ) * ey_num + 2 * ey + ex + 1; node[4] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex; node[5] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex + 1; node[6] = ( 3 * ey ) * ey_num + 2 * ey + 2 * ex + 2; node[7] = ( 3 * ey + 2 ) * ey_num + 2 * ey + ex + 2; for ( qx = 0; qx < quad_num; qx++ ) { xq = ( ( 1.0 - abscissa[qx] ) * x[e] + ( 1.0 + abscissa[qx] ) * x[w] ) / 2.0; for ( qy = 0; qy < quad_num; qy++ ) { yq = ( ( 1.0 - abscissa[qy] ) * y[n] + ( 1.0 + abscissa[qy] ) * y[s] ) / 2.0; wq = weight[qx] * ( x[e] - x[w] ) / 2.0 * weight[qy] * ( y[n] - y[s] ) / 2.0; v = basis_serene ( xq, yq, xw, ys, xe, yn, xx, yy ); uq = 0.0; for ( k = 0; k < 8; k++ ) { uq = uq + u[node[k]] * v[k]; } eq = exact ( xq, yq ); e2 = e2 + wq * pow ( uq - eq, 2 ); free ( v ); } } } } e2 = sqrt ( e2 ); return e2; # undef QUAD_NUM } /******************************************************************************/ int *i4vec_zero_new ( int n ) /******************************************************************************/ /* Purpose: I4VEC_ZERO_NEW creates and zeroes an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 05 September 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, int I4VEC_ZERO_NEW[N], a vector of zeroes. */ { int *a; int i; a = ( int * ) malloc ( n * sizeof ( int ) ); for ( i = 0; i < n; i++ ) { a[i] = 0; } return a; } /******************************************************************************/ double not1 ( double x1, double x2, double x3 ) /******************************************************************************/ /* Purpose: NOT1 evaluates a factor for serendipity basis functions. Discussion: not1(x1,x2,x3) evaluates at the point x1, the basis factor that is 0 at x2 and 1 at x3: not1(x1,x2,x3) = ( x1 - x2 ) / ( x3 - x2 ) Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, double X1, the evaluation point. Input, double X2, X3, values that define the factor. Output, double NOT1, the value of the basis function factor. */ { double value; value = ( x1 - x2 ) / ( x3 - x2 ); return value; } /******************************************************************************/ double not1d ( double x2, double x3 ) /******************************************************************************/ /* Purpose: NOT1D differentiates a factor for serendipity basis functions. Discussion: not1(x1,x2,x3) evaluates at the point x1, the basis factor that is 0 at x2 and 1 at x3: not1(x1,x2,x3) = ( x1 - x2 ) / ( x3 - x2 ) This function returns the derivative of the factor with respect to x1: not1d(x1,x2,x3) = 1 / ( x3 - x2 ) Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, double X2, X3, values that define the factor. Output, double NOT1D, the derivative of the basis function factor. */ { double value; value = 1.0 / ( x3 - x2 ); return value; } /******************************************************************************/ double not2 ( double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4 ) /******************************************************************************/ /* Purpose: NOT2 evaluates a factor for serendipity basis functions. Discussion: not2(x1,y1,x2,y2,x3,y3,x4,y4) evaluates at the point (x1,y1), the basis factor that is 0 at (x2,y2) and (x3,y3) and 1 at (x4,y4): ( ( x1 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y1 - y2 ) ) / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) ) Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, double X1, Y2, the evaluation point. Input, double X2, Y2, X3, Y3, values that define the factor. Output, double NOT2, the value of the basis function factor. */ { double value; value = ( ( x1 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y1 - y2 ) ) / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) ); return value; } /******************************************************************************/ double not2dx ( double x2, double y2, double x3, double y3, double x4, double y4 ) /******************************************************************************/ /* Purpose: NOT2DX evaluates a factor for serendipity basis functions. Discussion: not2(x1,y1,x2,y2,x3,y3,x4,y4) evaluates at the point (x1,y1), the basis factor that is 0 at (x2,y2) and (x3,y3) and 1 at (x4,y4): ( ( x1 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y1 - y2 ) ) / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) ) not2dx returns the derivative of this function with respect to X1. Licensing: This code is distributed under the MIT license. Modified: 07 July 2014 Author: John Burkardt Parameters: Input, double X2, Y2, X3, Y3, values that define the factor. Output, double NOT2DX, the derivative of the basis function factor with respect to X1. */ { double value; value = ( 1.0 * ( y3 - y2 ) + 0.0 ) / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) ); return value; } /******************************************************************************/ double not2dy ( double x2, double y2, double x3, double y3, double x4, double y4 ) /******************************************************************************/ /* Purpose: NOT2DY evaluates a factor for serendipity basis functions. Discussion: not2(x1,y1,x2,y2,x3,y3,x4,y4) evaluates at the point (x1,y1), the basis factor that is 0 at (x2,y2) and (x3,y3) and 1 at (x4,y4): ( ( x1 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y1 - y2 ) ) / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) ) not2dy returns the derivatives of this function with respect to Y1. Licensing: This code is distributed under the MIT license. Modified: 01 July 2014 Author: John Burkardt Parameters: Input, double X2, Y2, X3, Y3, values that define the factor. Output, double NOT2DY, the derivative of the basis function factor with respect to Y1. */ { double value; value = ( 0.0 - ( x3 - x2 ) * 1.0 ) / ( ( x4 - x2 ) * ( y3 - y2 ) - ( x3 - x2 ) * ( y4 - y2 ) ); return value; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1]. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { const int i4_huge = 2147483647; int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT_SOME prints some of an R8MAT. 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 MIT license. Modified: 26 June 2013 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, double A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2hi; int i2lo; 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; } /* Print the columns of the matrix, in strips of 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ double *r8mat_solve2 ( int n, double a[], double b[], int *ierror ) /******************************************************************************/ /* Purpose: R8MAT_SOLVE2 computes the solution of an N by N linear system. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. The linear system may be represented as A*X = B If the linear system is singular, but consistent, then the routine will still produce a solution. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int N, the number of equations. Input/output, double A[N*N]. On input, A is the coefficient matrix to be inverted. On output, A has been overwritten. Input/output, double B[N]. On input, B is the right hand side of the system. On output, B has been overwritten. Output, double R8MAT_SOLVE2[N], the solution of the linear system. Output, int *IERROR. 0, no error detected. 1, consistent singularity. 2, inconsistent singularity. */ { double amax; int i; int imax; int j; int k; int *piv; double *x; *ierror = 0; piv = i4vec_zero_new ( n ); x = r8vec_zero_new ( n ); /* Process the matrix. */ for ( k = 1; k <= n; k++ ) { /* In column K: Seek the row IMAX with the properties that: IMAX has not already been used as a pivot; A(IMAX,K) is larger in magnitude than any other candidate. */ amax = 0.0; imax = 0; for ( i = 1; i <= n; i++ ) { if ( piv[i-1] == 0 ) { if ( amax < fabs ( a[i-1+(k-1)*n] ) ) { imax = i; amax = fabs ( a[i-1+(k-1)*n] ); } } } /* If you found a pivot row IMAX, then, eliminate the K-th entry in all rows that have not been used for pivoting. */ if ( imax != 0 ) { piv[imax-1] = k; for ( j = k+1; j <= n; j++ ) { a[imax-1+(j-1)*n] = a[imax-1+(j-1)*n] / a[imax-1+(k-1)*n]; } b[imax-1] = b[imax-1] / a[imax-1+(k-1)*n]; a[imax-1+(k-1)*n] = 1.0; for ( i = 1; i <= n; i++ ) { if ( piv[i-1] == 0 ) { for ( j = k+1; j <= n; j++ ) { a[i-1+(j-1)*n] = a[i-1+(j-1)*n] - a[i-1+(k-1)*n] * a[imax-1+(j-1)*n]; } b[i-1] = b[i-1] - a[i-1+(k-1)*n] * b[imax-1]; a[i-1+(k-1)*n] = 0.0; } } } } /* Now, every row with nonzero PIV begins with a 1, and all other rows are all zero. Begin solution. */ for ( j = n; 1 <= j; j-- ) { imax = 0; for ( k = 1; k <= n; k++ ) { if ( piv[k-1] == j ) { imax = k; } } if ( imax == 0 ) { x[j-1] = 0.0; if ( b[j-1] == 0.0 ) { *ierror = 1; printf ( "\n" ); printf ( "R8MAT_SOLVE2 - Warning:\n" ); printf ( " Consistent singularity, equation = %d\n", j ); } else { *ierror = 2; printf ( "\n" ); printf ( "R8MAT_SOLVE2 - Warning:\n" ); printf ( " Inconsistent singularity, equation = %d\n", j ); } } else { x[j-1] = b[imax-1]; for ( i = 1; i <= n; i++ ) { if ( i != imax ) { b[i-1] = b[i-1] - a[i-1+(j-1)*n] * x[j-1]; } } } } free ( piv ); return x; } /******************************************************************************/ double *r8mat_zero_new ( int m, int n ) /******************************************************************************/ /* Purpose: R8MAT_ZERO_NEW returns a new zeroed R8MAT. Licensing: This code is distributed under the MIT license. Modified: 26 September 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Output, double R8MAT_ZERO[M*N], the new zeroed matrix. */ { double *a; int i; int j; a = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a[i+j*m] = 0.0; } } return a; } /******************************************************************************/ double *r8vec_linspace_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: R8VEC_LINSPACE_NEW creates a vector of linearly spaced values. Discussion: An R8VEC is a vector of R8's. 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. In other words, the interval is divided into N-1 even subintervals, and the endpoints of intervals are used as the points. Licensing: This code is distributed under the MIT license. Modified: 29 March 2011 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A, B, the first and last entries. Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. */ { int i; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { x[0] = ( a + b ) / 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = ( ( double ) ( n - 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( n - 1 ); } } return x; } /******************************************************************************/ double r8vec_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SUM returns the sum of an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 26 August 2008 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ double *r8vec_zero_new ( int n ) /******************************************************************************/ /* Purpose: R8VEC_ZERO_NEW creates and zeroes an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 25 March 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, double R8VEC_ZERO_NEW[N], a vector of zeroes. */ { double *a; int i; a = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { a[i] = 0.0; } return a; } /******************************************************************************/ 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 MIT 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 } /******************************************************************************/ double *wathen ( int nx, int ny, int n ) /******************************************************************************/ /* Purpose: WATHEN returns the WATHEN matrix. Discussion: The Wathen matrix is a finite element matrix which is sparse. The entries of the matrix depend in part on a physical quantity related to density. That density is here assigned random values between 0 and 100. The matrix order N is determined by the input quantities NX and NY, which would usually be the number of elements in the X and Y directions. The value of N is N = 3*NX*NY + 2*NX + 2*NY + 1, and sufficient storage in A must have been set aside to hold the matrix. A is the consistent mass matrix for a regular NX by NY grid of 8 node serendipity elements. The local element numbering is 3--2--1 | | 4 8 | | 5--6--7 Here is an illustration for NX = 3, NY = 2: 23-24-25-26-27-28-29 | | | | 19 20 21 22 | | | | 12-13-14-15-16-17-18 | | | | 8 9 10 11 | | | | 1--2--3--4--5--6--7 For this example, the total number of nodes is, as expected, N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29 Properties: A is symmetric positive definite for any positive values of the density RHO(NX,NY), which is here given the value 1. The problem could be reprogrammed so that RHO is nonconstant, but positive. Licensing: This code is distributed under the MIT license. Modified: 02 July 2014 Author: John Burkardt Reference: Nicholas Higham, Algorithm 694: A Collection of Test Matrices in MATLAB, ACM Transactions on Mathematical Software, Volume 17, Number 3, September 1991, pages 289-305. Andrew Wathen, Realistic eigenvalue bounds for the Galerkin mass matrix, IMA Journal of Numerical Analysis, Volume 7, Number 4, October 1987, pages 449-457. Parameters: Input, int NX, NY, values which determine the size of A. Input, int N, the order of the matrix. Output, double WATHEN[N*N], the matrix. */ { double *a; static double em[8*8] = { 6.0, -6.0, 2.0, -8.0, 3.0, -8.0, 2.0, -6.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, 2.0, -6.0, 6.0, -6.0, 2.0, -8.0, 3.0, -8.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, -8.0, 16.0, 3.0, -8.0, 2.0, -6.0, 6.0, -6.0, 2.0, -8.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0, -6.0, 20.0, 2.0, -8.0, 3.0, -8.0, 2.0, -6.0, 6.0, -6.0, -6.0, 20.0, -8.0, 16.0, -8.0, 20.0, -6.0, 32.0 }; int i; int j; int kcol; int krow; int node[8]; double rho; int seed; a = ( double * ) malloc ( n * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < n; i++ ) { a[i+j*n] = 0.0; } } seed = 123456789; for ( j = 1; j <= ny; j++ ) { for ( i = 1; i <= nx; i++ ) { /* For the element (I,J), determine the indices of the 8 nodes. */ node[0] = 3 * j * nx + 2 * j + 2 * i; node[1] = node[0] - 1; node[2] = node[0] - 2; node[3] = ( 3 * j - 1 ) * nx + 2 * j + i - 2; node[4] = ( 3 * j - 3 ) * nx + 2 * j + 2 * i - 4; node[5] = node[4] + 1; node[6] = node[4] + 2; node[7] = node[3] + 1; rho = 100.0 * r8_uniform_01 ( &seed ); for ( krow = 0; krow < 8; krow++ ) { for ( kcol = 0; kcol < 8; kcol++ ) { a[node[krow]+node[kcol]*n] = a[node[krow]+node[kcol]*n] + rho * em[krow+kcol*8]; } } } } return a; }