# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include "fem2d_bvp_linear.h"
/******************************************************************************/
double *fem2d_bvp_linear ( 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[] )
/******************************************************************************/
/*
Purpose:
FEM2D_BVP_LINEAR solves a boundary value problem on a rectangle.
Discussion:
The procedure uses the finite element method, with piecewise linear 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.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
20 June 2014
Author:
John Burkardt
Parameters:
Input, int NX, NY, the number of X and Y grid values.
Input, double A ( double X, double Y ), evaluates a(x,y);
Input, double C ( double X, double Y ), evaluates c(x,y);
Input, double F ( double X, double Y ), evaluates f(x,y);
Input, double X[NX], Y[NY], the grid coordinates.
Output, double FEM1D_BVP_LINEAR[NX*NY], 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 *amat;
double aq;
double *b;
double cq;
int e;
int ex;
int ey;
double fq;
int i;
int ierror;
int ii;
int j;
int jj;
int k;
int mn;
int n;
int ne;
int nw;
int quad_num = QUAD_NUM;
int qx;
int qy;
int s;
int se;
int sw;
int w;
double weight[QUAD_NUM] = {
0.555555555555555555555555555556,
0.888888888888888888888888888889,
0.555555555555555555555555555556 };
double wq;
double *u;
double vne;
double vnex;
double vney;
double vnw;
double vnwx;
double vnwy;
double vse;
double vsex;
double vsey;
double vsw;
double vswx;
double vswy;
double xe;
double xq;
double xw;
double yn;
double yq;
double ys;
mn = nx * ny;
amat = ( double * ) malloc ( mn * mn * sizeof ( double ) );
b = ( double * ) malloc ( mn * sizeof ( double ) );
for ( jj = 0; jj < mn; jj++ )
{
for ( ii = 0; ii < mn; ii++ )
{
amat[ii+jj*mn] = 0.0;
}
}
for ( ii = 0; ii < mn; ii++ )
{
b[ii] = 0.0;
}
for ( ex = 0; ex < nx - 1; ex++ )
{
w = ex;
e = ex + 1;
xw = x[w];
xe = x[e];
for ( ey = 0; ey < ny - 1; ey++ )
{
s = ey;
n = ey + 1;
ys = y[s];
yn = y[n];
sw = ey * nx + ex;
se = ey * nx + ex + 1;
nw = ( ey + 1 ) * nx + ex;
ne = ( ey + 1 ) * nx + ex + 1;
for ( qx = 0; qx < quad_num; qx++ )
{
xq = ( ( 1.0 - abscissa[qx] ) * xw
+ ( 1.0 + abscissa[qx] ) * xe )
/ 2.0;
for ( qy = 0; qy < quad_num; qy++ )
{
yq = ( ( 1.0 - abscissa[qy] ) * ys
+ ( 1.0 + abscissa[qy] ) * yn )
/ 2.0;
wq = weight[qx] * ( xe - xw ) / 2.0
* weight[qy] * ( yn - ys ) / 2.0;
aq = a ( xq, yq );
cq = c ( xq, yq );
fq = f ( xq, yq );
vsw = ( xe - xq ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
vswx = (-1.0 ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
vswy = ( xe - xq ) / ( xe - xw ) * (-1.0 ) / ( yn - ys );
vse = ( xq - xw ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
vsex = ( 1.0 ) / ( xe - xw ) * ( yn - yq ) / ( yn - ys );
vsey = ( xq - xw ) / ( xe - xw ) * (-1.0 ) / ( yn - ys );
vnw = ( xe - xq ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
vnwx = (-1.0 ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
vnwy = ( xe - xq ) / ( xe - xw ) * ( 1.0 ) / ( yn - ys );
vne = ( xq - xw ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
vnex = ( 1.0 ) / ( xe - xw ) * ( yq - ys ) / ( yn - ys );
vney = ( xq - xw ) / ( xe - xw ) * ( 1.0 ) / ( yn - ys );
amat[sw+sw*mn] = amat[sw+sw*mn] + wq * ( vswx * aq * vswx
+ vswy * aq * vswy
+ vsw * cq * vsw );
amat[sw+se*mn] = amat[sw+se*mn] + wq * ( vswx * aq * vsex
+ vswy * aq * vsey
+ vsw * cq * vse );
amat[sw+nw*mn] = amat[sw+nw*mn] + wq * ( vswx * aq * vnwx
+ vswy * aq * vnwy
+ vsw * cq * vnw );
amat[sw+ne*mn] = amat[sw+ne*mn] + wq * ( vswx * aq * vnex
+ vswy * aq * vney
+ vsw * cq * vne );
b[sw] = b[sw] + wq * ( vsw * fq );
amat[se+sw*mn] = amat[se+sw*mn] + wq * ( vsex * aq * vswx
+ vsey * aq * vswy
+ vse * cq * vsw );
amat[se+se*mn] = amat[se+se*mn] + wq * ( vsex * aq * vsex
+ vsey * aq * vsey
+ vse * cq * vse );
amat[se+nw*mn] = amat[se+nw*mn] + wq * ( vsex * aq * vnwx
+ vsey * aq * vnwy
+ vse * cq * vnw );
amat[se+ne*mn] = amat[se+ne*mn] + wq * ( vsex * aq * vnex
+ vsey * aq * vney
+ vse * cq * vne );
b[se] = b[se] + wq * ( vse * fq );
amat[nw+sw*mn] = amat[nw+sw*mn] + wq * ( vnwx * aq * vswx
+ vnwy * aq * vswy
+ vnw * cq * vsw );
amat[nw+se*mn] = amat[nw+se*mn] + wq * ( vnwx * aq * vsex
+ vnwy * aq * vsey
+ vnw * cq * vse );
amat[nw+nw*mn] = amat[nw+nw*mn] + wq * ( vnwx * aq * vnwx
+ vnwy * aq * vnwy
+ vnw * cq * vnw );
amat[nw+ne*mn] = amat[nw+ne*mn] + wq * ( vnwx * aq * vnex
+ vnwy * aq * vney
+ vnw * cq * vne );
b[nw] = b[nw] + wq * ( vnw * fq );
amat[ne+sw*mn] = amat[ne+sw*mn] + wq * ( vnex * aq * vswx
+ vney * aq * vswy
+ vne * cq * vsw );
amat[ne+se*mn] = amat[ne+se*mn] + wq * ( vnex * aq * vsex
+ vney * aq * vsey
+ vne * cq * vse );
amat[ne+nw*mn] = amat[ne+nw*mn] + wq * ( vnex * aq * vnwx
+ vney * aq * vnwy
+ vne * cq * vnw );
amat[ne+ne*mn] = amat[ne+ne*mn] + wq * ( vnex * aq * vnex