weight.txt.txt

边界单元法应用于电磁场及弹性力学中权函数的代码

function [w, dwdx, dwdxx] =3D Weight(type, para, di, dmi)
% EVALUATE WEIGHT FUNCTION
%
% SYNTAX: [w, dwdr, dwdrr] =3D GaussWeight(type, para, di, dmi)
%
% INPUT PARAMETERS
%    type - Type of weight function
%    para - Weight function parameter
%    di   - Distance
%    dmi  - Support size
% OUTPUT PARAMETERS
%    w    - Value of weight function at r
%    dwdx - Value of first order derivative of weight function with =
respect to x at r
%    dwdxx- Value of Second order derivative of weight function with =
respect to x at r
%
r  =3D abs(di) / dmi;

if (di >=3D 0.0)
   drdx =3D 1.0/dmi;
else
   drdx =3D -1.0/dmi;
end

% EVALUATE WEIGHT FUNCTION AND ITS FIRST AND SECOND ORDER OF DERIVATIVES =
WITH RESPECT r AT r
if (type =3D=3D 'GAUSS')
   [w,dwdr,dwdrr] =3D Gauss(para,r);
elseif (type =3D=3D 'CUBIC')
   [w,dwdr,dwdrr] =3D Cubic(r);
elseif (type =3D=3D 'SPLIN')
   [w,dwdr,dwdrr] =3D Spline(r);
elseif (type =3D=3D 'CSRBF')
   [w,dwdr,dwdrr] =3D CSRBF2(r);
else
   error('Invalid type of weight function.');
end

dwdx  =3D dwdr * drdx;
dwdxx =3D dwdrr * drdx * drdx;


function [w,dwdr,dwdrr] =3D Gauss(beta,r)
if (r>1.0)
   w     =3D 0.0;
   dwdr  =3D 0.0;
   dwdrr =3D 0.0;
else
   b2 =3D beta*beta;
   r2 =3D r*r;
   eb2 =3D exp(-b2);

   w     =3D (exp(-b2*r2) - eb2) / (1.0 - eb2);
   dwdr  =3D -2*b2*r*exp(-b2*r2) / (1.0 - eb2);
   dwdrr =3D -2*b2*exp(-b2*r2)*(1-2*b2*r2) / (1.0 - eb2);
end

function [w,dwdr,dwdrr] =3D Cubic(r)
if (r>1.0)
   w     =3D 0.0;
   dwdr  =3D 0.0;
   dwdrr =3D 0.0;
else
   w     =3D 1-6*r^2+8*r^3-3*r^4;
   dwdr  =3D -12*r+24*r^2-12*r^3;
   dwdrr =3D -12+48*r-36*r^2;
end

function [w,dwdr,dwdrr] =3D Spline(r)
if (r>1.0)
   w     =3D 0.0;
   dwdr  =3D 0.0;
   dwdrr =3D 0.0;
elseif (r<=3D0.5)
   w     =3D 2/3 - 4*r^2 + 4*r^3;
   dwdr  =3D -8*r + 12*r^2;
   dwdrr =3D -8 + 24*r;
else
   w     =3D 4/3 - 4*r + 4*r^2 - 4*r^3/3;
   dwdr  =3D -4 + 8*r -4*r^2;
   dwdrr =3D 8 - 8*r;
end

function [w,dwdr,dwdrr] =3D CSRBF2(r)
if (r>1.0)
   w     =3D 0.0;
   dwdr  =3D 0.0;
   dwdrr =3D 0.0;
else
	w     =3D (1-r)^6*(6+36*r+82*r^2+72*r^3+30*r^4+5*r^5);
	dwdr  =3D 11*r*(r+2)*(5*r^3+15*r^2+18*r+4)*(r-1)^5;
   dwdrr =3D 22*(25*r^5+100*r^4+142*r^3+68*r^2-16*r-4)*(r-1)^4;
end