C 代码 应用有限元法（FEM），与偶然性元素， 矩形中的二维边界值问题 （BVP）.rar

# include <stdlib.h> # include <stdio.h> # include <math.h> # include <time.h> # include <string.h> # 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 (