我正在尝试将 Fortran 中的代码翻译成 C++。显然，我在某个地方误解了 Fortran 语法，这就是我在第 94 行遇到问题的原因。我查看了第 94 行组件的值，意识到我很可能错误地计算了二维数组 A .但是我多次检查了它的计算，并没有发现任何错误。请告诉我需要做什么才能使代码正常工作？ （我使用了调试器，我倾向于相信我因为语言语法的差异而搞砸了一些东西）
С++：
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
void initl(double& sxxs, double& sxxn, double& syys, double& syyn, double& sxys, double& sxyn, double& uxs, double& uxn, double& uys, double& uyn)
{
sxxs = 0.0;
sxxn = 0.0;
syys = 0.0;
syyn = 0.0;
sxys = 0.0;
sxyn = 0.0;
uxs = 0.0;
uxn = 0.0;
uys = 0.0;
uyn = 0.0;
}
void coeff(double x, double y, double cx, double cy, double a, double cosb, double sinb, int msym, double pi, double con, double cond, double pr1, double pr2, double pr3, double& uxs, double& uxn, double& uys, double& uyn, double& sxxs, double& sxxn, double& syys, double& syyn, double& sxys, double& sxyn) {
double cos2b, sin2b, xb, yb, r1s, r2s, fl1, fl2, fb1, fb2, fb3, fb4, fb5, uxps, uxpn, uyps, uypn;
double sxxps, sxxpn, syyps, syypn, sxyps, sxypn;
cos2b = cosb * cosb - sinb * sinb;
sin2b = 2.0 * sinb * cosb;
xb = (x - cx) * cosb + (y - cy) * sinb;
yb = -(x - cx) * sinb + (y - cy) * cosb;
r1s = (xb - a) * (xb - a) + yb * yb;
r2s = (xb + a) * (xb + a) + yb * yb;
fl1 = 0.5 * log(r1s);//log - натуральный логарифм
fl2 = 0.5 * log(r2s);
fb2 = con * (fl1 - fl2);
if (yb != 0.0) {
fb3 = -con * (atan((xb + a) / yb) - atan((xb - a) / yb));
}
else
{
fb3 = 0.0;
if (abs(xb) < a) {
fb3 = con * pi;
}
}
fb1 = yb * fb3 + con * ((xb - a) * fl1 - (xb + a) * fl2);
fb4 = con * (yb / r1s - yb / r2s);
fb5 = con * ((xb - a) / r1s - (xb + a) / r2s);
uxps = cond * (pr3 * cosb * fb1 + yb * (sinb * fb2 + cosb * fb3));
uxpn = cond * (-pr3 * sinb * fb1 - yb * (cosb * fb2 - sinb * fb3));
uyps = cond * (pr3 * sinb * fb1 - yb * (cosb * fb2 - sinb * fb3));
uypn = cond * (pr3 * cosb * fb1 - yb * (sinb * fb2 + cosb * fb3));
sxxps = fb2 + pr2 * (cos2b * fb2 - sin2b * fb3) + yb * (cos2b * fb4 + sin2b * fb5);
sxxpn = fb3 - pr1 * (sin2b * fb2 + cos2b * fb3) + yb * (sin2b * fb4 - cos2b * fb5);
syyps = fb2 - pr2 * (cos2b * fb2 - sin2b * fb3) - yb * (cos2b * fb4 + sin2b * fb5);
syypn = fb3 + pr1 * (sin2b * fb2 + cos2b * fb3) - yb * (sin2b * fb4 - cos2b * fb5);
sxyps = pr2 * (sin2b * fb2 + cos2b * fb3) + yb * (sin2b * fb4 - cos2b * fb5);
sxypn = pr1 * (cos2b * fb2 - sin2b * fb3) - yb * (cos2b * fb4 + sin2b * fb5);
uxs = uxs + msym * uxps;
uxn = uxn + uxpn;
uys = uys + msym * uyps;
uyn = uyn + uypn;
sxxs = sxxs + msym * sxxps;
sxxn = sxxn + sxxpn;
syys = syys + msym * syyps;
syyn = syyn + syypn;
sxys = sxys + msym * sxyps;
sxyn = sxyn + sxypn;
}
void solve(int n, vector<vector<double>>& a, vector<double>& b, vector<double>& x) {
int nb, l;
double xm;
nb = n - 1;
for (int j = 0; j < nb; ++j) {// 1 -> 0; "<=" -> "<"
l = j + 1;
for (int jj = l; jj < n; ++jj) {// "<=" -> "<"
xm = a[jj][j] / a[j][j];
for (int i = j; i < n; ++i) {// "<=" -> "<"
a[jj][i] = a[jj][i] - a[j][i] * xm;
}
b[jj] = b[jj] - b[j] * xm;
}
}
int jj;
double summ;
summ = 0.0;
x[n-1] = b[n-1] / a[n-1][n-1];
for (int j = 0; j < nb; ++j) {// 1 -> 0; "<=" -> "<"
jj = n - j;
l = jj + 1;
summ = 0.0;
for (int i = l; i < n; ++i) {// "<=" -> "<"
summ = summ + a[jj][i] * x[i];
}
x[jj] = (b[jj] - summ) / a[jj][jj];
}
}
int main()
{
int numbs, numos, ksym;
cin >> numbs >> numos >> ksym;
double pr, e, xsym, ysym;
cin >> pr >> e >> xsym >> ysym;
double pxx, pyy, pxy;
cin >> pxx >> pyy >> pxy;
cout << " number of straight-line segments (each counting at least one boundary element) used to define boundaries = " << numbs << endl;
cout << " number of straight-line segments used to specify other locations (i.e., not on a boundary) where results are to be found = " << numos << endl;
cout << " the lines x = xs = " << xsym << " and y = ys = " << ysym << " are lines of symmetry." << endl;
cout << " poissons ratio = " << pr << endl;
cout << " youngs modulus = " << e << endl;
cout << " xx-component of field stress = " << pxx << endl;
cout << " yy-component of field stress = " << pyy << endl;
cout << " xy-component of field stress = " << pxy << endl;
/*cout << "" << << endl;*/
double pi, con, cond, pr1, pr2, pr3;
pi = 4.0 * atan(1.0);
con = 1.0 / (4.0 * pi * (1.0 - pr));
cond = (1.0 + pr) / e;
pr1 = 1.0 - 2.0 * pr;
pr2 = 2.0 * (1.0 - pr);
pr3 = 3.0 - 4.0 * pr;
//define locations, sizes, orientations and boundary conditions of boundary elements.
int numbe, num, kode, m, mn, ms;
vector<int> kod;
double xbeg, ybeg, xend, yend, bvs, bvn, xd, yd, sw;
vector<double> xm, ym, a, sinbet, cosbet, b;
numbe = 0;
for (int n = 1; n <= numbs; ++n) {
cin >> num >> xbeg >> ybeg >> xend >> yend >> kode >> bvs >> bvn;
xd = (xend - xbeg) / num;
yd = (yend - ybeg) / num;
sw = sqrt(xd * xd + yd * yd);
for (int ne = 1; ne <= num; ++ne) {
numbe = numbe + 1;
xm.push_back(xbeg + 0.5 * (2.0 * ne - 1.0) * xd);
ym.push_back(ybeg + 0.5 * (2.0 * ne - 1.0) * yd);
a.push_back(0.5 * sw);
sinbet.push_back(yd / sw);
cosbet.push_back(xd / sw);
kod.push_back(kode);
b.push_back(bvs);
b.push_back(bvn);
}
}
cout << "boundary element data." << endl;
cout << "element kode x (center) y (center) length angle us or sigma-sun or sigma-n" << endl;
double size, angle;
for (int m = 0; m < numbe; ++m) {
size = 2.0 * a[m];
angle = 180.0 * atan2(sinbet[m], cosbet[m]) / pi;
cout << m + 1 << " " << kod[m] << " " << xm[m] << " " << ym[m] << " " << size << " " << angle << " " << b[2 * m] << " " << b[2 * m + 1] << endl;
}
//adjust stress boundary values to account for initial stresses.
int nn, ns;
double cosb, sinb, sigs, sign;
for (int n = 0; n < numbe; ++n) {// 1 -> 0; "<=" -> "<"
nn = 2 * n + 1;// 2 * i
ns = nn - 1;
cosb = cosbet[n];
sinb = sinbet[n];
sigs = (pyy - pxx) * sinb * cosb + pxy * (cosb * cosb - sinb * sinb);
sign = pxx * sinb * sinb - 2. * pxy * sinb * cosb + pyy * cosb * cosb;
b[ns] = b[ns] - sigs;
b[nn] = b[nn] - sign;
b[nn] = b[nn] - sign;
b[ns] = b[ns] - sigs;
}
//compute influence coefficients and set up system of algebraic equations.
int in, is, jn, js;
double xi, yi, cosbi, sinbi, xj, yj, cosbj, sinbj, aj;
double sxxs, sxxn, syys, syyn, sxys, sxyn, uxs, uxn, uys, uyn;
vector<vector<double>> c(12, vector<double>(12, 0));
for (int i = 0; i < numbe; ++i) {// 1 -> 0; "<=" -> "<"
in = 2 * i + 1;// 2 * i
is = in - 1;
xi = xm[i];
yi = ym[i];
cosbi = cosbet[i];
sinbi = sinbet[i];
kode = kod[i];
for (int j = 0; j < numbe; ++j) {// 1 -> 0; "<=" -> "<"
jn = 2 * j + 1;// 2 * j
js = jn - 1;
initl(sxxs, sxxn, syys, syyn, sxys, sxyn, uxs, uxn, uys, uyn);
xj = xm[j];
yj = ym[j];
cosbj = cosbet[j];
sinbj = sinbet[j];
aj = a[j];
coeff(xi, yi, xj, yj, aj, cosbj, sinbj, +1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
if (ksym == 1) {
continue;
}
else if (ksym == 2) {
xj = 2. * xsym - xm[j];
coeff(xi, yi, xj, yj, aj, cosbj, -sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
continue;
}
else if (ksym == 3) {
yj = 2. * ysym - ym[j];
coeff(xi, yi, xj, yj, aj, -cosbj, sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
continue;
}
else if (ksym == 4) {
xj = 2. * xsym - xm[j];
coeff(xi, yi, xj, yj, aj, cosbj, -sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
xj = xm[j];
yj = 2. * ysym - ym[j];
coeff(xi, yi, xj, yj, aj, -cosbj, sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
xj = 2. * xsym - xm[j];
coeff(xi, yi, xj, yj, aj, -cosbj, -sinbj, +1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
}
if (ksym == 1) {
c[is][js] = (syys - sxxs) * sinbi * cosbi + sxys * (cosbi * cosbi - sinbi * sinbi);
c[is][jn] = (syyn - sxxn) * sinbi * cosbi + sxyn * (cosbi * cosbi - sinbi * sinbi);
c[in][js] = sxxs * sinbi * sinbi - 2. * sxys * sinbi * cosbi + syys * cosbi * cosbi;
c[in][jn] = sxxn * sinbi * sinbi - 2. * sxyn * sinbi * cosbi + syyn * cosbi * cosbi;
continue;
}
else if (ksym == 2) {
c[is][js] = uxs * cosbi + uys * sinbi;
c[is][jn] = uxn * cosbi + uyn * sinbi;
c[in][js] = -uxs * sinbi + uys * cosbi;
c[in][jn] = -uxn * sinbi + uyn * cosbi;
continue;
}
else if (ksym == 3) {
c[is][js] = uxs * cosbi + uys * sinbi;
c[is][jn] = uxn * cosbi + uyn * sinbi;
c[in][js] = sxxs * sinbi * sinbi - 2. * sxys * sinbi * cosbi + syys * cosbi * cosbi;
c[in][jn] = sxxn * sinbi * sinbi - 2. * sxyn * sinbi * cosbi + syyn * cosbi * cosbi;
continue;
}
else if (ksym == 4) {
c[is][js] = (syys - sxxs) * sinbi * cosbi + sxys * (cosbi * cosbi - sinbi * sinbi);
c[is][jn] = (syyn - sxxn) * sinbi * cosbi + sxyn * (cosbi * cosbi - sinbi * sinbi);
c[in][js] = -uxs * sinbi + uys * cosbi;
c[in][jn] = -uxn * sinbi + uyn * cosbi;
}
}
}
//solve system of algebraic equations.
int n;
vector<double> p;
n = 2 * numbe;
solve(n, c, b, a);
// compute boundary displacements and stresses.
cout << endl << "displacements and stresses at midpoints of boundary elements." << endl << endl << "element ux uy us un sigxx sigyy sigxy sigma-s sigma-n sigma-t" << endl;
double ux, uy, sigxx, sigyy, sigxy, us, un, sigt;
for (int i = 0; i < numbe; ++i) {// 1 -> 0; "<=" -> "<"
xi = xm[i];
yi = ym[i];
cosbi = cosbet[i];
sinbi = sinbet[i];
ux = 0.0;
uy = 0.0;
sigxx = pxx;
sigyy = pyy;
sigxy = pxy;
for (int j = 0; j < numbe; ++j) {// 1 -> 0; "<=" -> "<"
jn = 2 * j;
js = jn - 1;
initl(sxxs, sxxn, syys, syyn, sxys, sxyn, uxs, uxn, uys, uyn);
xj = xm[j];
yj = ym[j];
aj = a[j];
cosbj = cosbet[j];
sinbj = sinbet[j];
coeff(xi, yi, xj, yj, aj, cosbj, sinbj, +1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
if (ksym == 1) {
continue;
}
if (ksym == 2) {
xj = 2.0 * xsym - xm[j];
coeff(xi, yi, xj, yj, aj, cosbj, -sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
continue;
}
if (ksym == 3) {
yj = 2.0 * ysym - ym[j];
coeff(xi, yi, xj, yj, aj, -cosbj, sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
continue;
}
if (ksym == 4) {
xj = 2.0 * xsym - xm[j];
coeff(xi, yi, xj, yj, aj, cosbj, -sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
xj = xm[j];
yj = 2.0 * ysym - ym[j];
coeff(xi, yi, xj, yj, aj, -cosbj, sinbj, -1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
xj = 2.0 * xsym - xm[j];
coeff(xi, yi, xj, yj, aj, -cosbj, -sinbj, +1, pi, con, cond, pr1, pr2, pr3, uxs, uxn, uys, uyn, sxxs, sxxn, syys, syyn, sxys, sxyn);
}
ux = ux + uxs * p[js] + uxn * p[jn];
uy = uy + uys * p[js] + uyn * p[jn];
sigxx = sigxx + sxxs * p[js] + sxxn * p[jn];
sigyy = sigyy + syys * p[js] + syyn * p[jn];
sigxy = sigxy + sxys * p[js] + sxyn * p[jn];
}
us = ux * cosbi + uy * sinbi;
un = -ux * sinbi + uy * cosbi;
sigs = (sigyy - sigxx) * sinbi * cosbi + sigxy * (cosbi * cosbi - sinbi * sinbi);
sign = sigxx * sinbi * sinbi - 2.0 * sigxy * sinbi * cosbi + sigyy * cosbi * cosbi;
sigt = sigxx * cosbi * cosbi + 2.0 * sigxy * cosbi * sinbi + sigyy * sinbi * sinbi;
cout << i << " " << ux << " " << uy << " " << us << " " << un << " " << sigxx << " " << sigyy << " " << sigxy << " " << sigs << " " << sign << " " << sigt << endl;
}
}
Fortran：
PROGRAM TWOFS
!
COMMON/S1/PI,PR,PR1,PR2,PR3,CON,COND
COMMON/S2/SXXS,SXXN,SYYS,SYYN,SXYS,SXYN,UXS,UXN,UYS,UYN
COMMON/S3/C(100,100),B(100),P(100)
!
DIMENSION XM(50),YM(50),A(50),COSBET(50),SINBET(50),KOD(50)
DIMENSION TITLE(20)
!
READ (5,1) (TITLE(I),I=1,20)
WRITE (6,2) (TITLE(I),I=1,20)
READ (5,3) NUMBS,NUMOS,KSYM
READ (5,4) PR,E,XSYM,YSYM
READ (5,5) PXX,PYY,PXY
WRITE (6,6) NUMBS,NUMOS
GO TO (80,85,90,95), KSYM
80 WRITE (6,7)
GO TO 100
85 WRITE (6,8) XSYM
GO TO 100
90 WRITE (6,9) YSYM
GO TO 100
95 WRITE (6,10) XSYM,YSYM
!
100 CONTINUE
WRITE (6,11) PR,E
WRITE (6,12) PXX,PYY,PXY
!
PI=4.*ATAN(1.)
CON=1./(4.*PI*(1.-PR))
COND=(1.+PR)/E
PR1=1.-2.*PR
PR2=2.*(1.-PR)
PR3=3.-4.*PR
!
! DEFINE LOCATIONS,SIZES,ORIENTATIONS AND BOUNDARY CONDITIONS OF BOUNDARY ELEMENTS
!
NUMBE=0
DO 110 N=1,NUMBS
READ (5,14) NUM,XBEG,YBEG,XEND,YEND,KODE,BVS,BVN
XD=(XEND-XBEG)/NUM
YD=(YEND-YBEG)/NUM
SW=SQRT(XD*XD+YD*YD)
!
DO 110 NE=1,NUM
NUMBE=NUMBE+1
M=NUMBE
XM(M)=XBEG+0.5*(2.*NE-1.)*XD
YM(M)=YBEG+0.5*(2.*NE-1.)*YD
A(M)=0.5*SW
SINBET(M)=YD/SW
COSBET(M)=XD/SW
KOD(M)=KODE
MN=2*M
MS=MN-1
B(MS)=BVS
110 B(MN)=BVN
WRITE (6,13)
DO 115 M=1,NUMBE
SIZE=2.*A(M)
ANGLE=180.*ATAN2(SINBET(M),COSBET(M))/PI
WRITE (6,15) M,KOD(M),XM(M),YM(M),SIZE,ANGLE,B(2*M-1),B(2*M)
115 CONTINUE
!
! ADJUST STRESS BOUNDARY VALUES TO ACCOUNT FOR INITIAL STRESSES.
!
DO 150 N=1,NUMBE
NN=2*N
NS=NN-1
COSB=COSBET(N)
SINB=SINBET(N)
SIGS=(PYY-PXX)*SINB*COSB+PXY*(COSB*COSB-SINB*SINB)
SIG_N=PXX*SINB*SINB-2.*PXY*SINB*COSB+PYY*COSB*COSB
GO TO (120,150,130,140),KOD(N)
120 B(NS)=B(NS)-SIGS
B(NN)=B(NN)-SIG_N
GO TO 150
130 B(NN)=B(NN)-SIG_N
GO TO 150
140 B(NS)=B(NS)-SIGS
150 CONTINUE
!
! COMPUTE INFLUENCE COEFFICIENTS AND SET UP SYSTEM OF ALGEBRAIC EQUATIONS.
!
DO 300 I=1,NUMBE
IN=2*I
IS=IN-1
XI=XM(I)
YI=YM(I)
COSBI=COSBET(I)
SINBI=SINBET(I)
KODE=KOD(I)
!
DO 300 J=1,NUMBE
JN=2*J
JS=JN-1
CALL INITL
XJ=XM(J)
YJ=YM(J)
COSBJ=COSBET(J)
SINBJ=SINBET(J)
AJ=A(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,COSBJ,SINBJ,+1)
GO TO (240,210,220,230),KSYM
!
210 XJ=2.*XSYM-XM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,COSBJ,-SINBJ,-1)
GO TO 240
!
220 YJ=2.*YSYM-YM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,-COSBJ,SINBJ,-1)
GO TO 240
!
230 XJ=2.*XSYM-XM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,COSBJ,-SINBJ,-1)
XJ=XM(J)
YJ=2.*YSYM-YM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,-COSBJ,SINBJ,-1)
XJ=2.*XSYM-XM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,-COSBJ,-SINBJ,+1)
!
240 CONTINUE
GO TO (250,260,270,280),KODE
!
250 C(IS,JS)=(SYYS-SXXS)*SINBI*COSBI+SXYS*(COSBI*COSBI-SINBI*SINBI)
C(IS,JN)=(SYYN-SXXN)*SINBI*COSBI+SXYN*(COSBI*COSBI-SINBI*SINBI)
C(IN,JS)=SXXS*SINBI*SINBI-2.*SXYS*SINBI*COSBI+SYYS*COSBI*COSBI
C(IN,JN)=SXXN*SINBI*SINBI-2.*SXYN*SINBI*COSBI+SYYN*COSBI*COSBI
GO TO 300
!
260 C(IS,JS)=UXS*COSBI+UYS*SINBI
C(IS,JN)=UXN*COSBI+UYN*SINBI
C(IN,JS)=-UXS*SINBI+UYS*COSBI
C(IN,JN)=-UXN*SINBI+UYN*COSBI
GO TO 300
!
270 C(IS,JS)=UXS*COSBI+UYS*SINBI
C(IS,JN)=UXN*COSBI+UYN*SINBI
C(IN,JS)=SXXS*SINBI*SINBI-2.*SXYS*SINBI*COSBI+SYYS*COSBI*COSBI
C(IN,JN)=SXXN*SINBI*SINBI-2.*SXYN*SINBI*COSBI+SYYN*COSBI*COSBI
GO TO 300
!
280 C(IS,JS)=(SYYS-SXXS)*SINBI*COSBI+SXYS*(COSBI*COSBI-SINBI*SINBI)
C(IS,JN)=(SYYN-SXXN)*SINBI*COSBI+SXYN*(COSBI*COSBI-SINBI*SINBI)
C(IN,JS)=-UXS*SINBI+UYS*COSBI
C(IN,JN)=-UXN*SINBI+UYN*COSBI
!
300 CONTINUE
!
! SOLVE SYSTEM OF ALGEBRAIC EQUATIONS.
!
N=2*NUMBE
CALL SOLVE(N)
!
! COMPUTE BOUNDARY DISPLACEMENTS AND STRESSES.
!
WRITE (6,16)
DO 600 I=1,NUMBE
XI=XM(I)
YI=YM(I)
COSBI=COSBET(I)
SINBI=SINBET(I)
!
UX=0.
UY=0.
SIGXX=PXX
SIGYY=PYY
SIGXY=PXY
!
DO 570 J=1,NUMBE
JN=2*J
JS=JN-1
CALL INITL
XJ=XM(J)
YJ=YM(J)
AJ=A(J)
COSBJ=COSBET(J)
SINBJ=SINBET(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,COSBJ,SINBJ,+1)
GO TO (540,510,520,530),KSYM
!
510 XJ=2.*XSYM-XM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,COSBJ,-SINBJ,-1)
GO TO 540
!
520 YJ=2.*YSYM-YM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,-COSBJ,SINBJ,-1)
GO TO 540
!
530 XJ=2.*XSYM-XM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,COSBJ,-SINBJ,-1)
XJ=XM(J)
YJ=2.*YSYM-YM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,-COSBJ,SINBJ,-1)
XJ=2.*XSYM-XM(J)
CALL COEFF(XI,YI,XJ,YJ,AJ,-COSBJ,-SINBJ,+1)
!
540 CONTINUE
!
UX=UX+UXS*P(JS)+UXN*P(JN)
UY=UY+UYS*P(JS)+UYN*P(JN)
SIGXX=SIGXX+SXXS*P(JS)+SXXN*P(JN)
SIGYY=SIGYY+SYYS*P(JS)+SYYN*P(JN)
SIGXY=SIGXY+SXYS*P(JS)+SXYN*P(JN)
!
570 CONTINUE
!
US=UX*COSBI+UY*SINBI
UN=-UX*SINBI+UY*COSBI
SIGS=(SIGYY-SIGXX)*SINBI*COSBI+SIGXY*(COSBI*COSBI-SINBI*SINBI)
SIG_N=SIGXX*SINBI*SINBI-2.*SIGXY*SINBI*COSBI+SIGYY*COSBI*COSBI
SIGT=SIGXX*COSBI*COSBI+2.*SIGXY*COSBI*SINBI+SIGYY*SINBI*SINBI
!
WRITE (6,17) I,UX,UY,US,UN,SIGXX,SIGYY,SIGXY,SIGS,SIG_N,SIGT
!
600 CONTINUE
!
! COMPUTE DISPLACEMENTS AND STRESSES AT SPECIFIED POINTS IN BOD
!
IF (NUMOS.LE.0) GO TO 900
WRITE (6,18)
NPOINT=0
DO 890 N=1,NUMOS
READ (5,19) XBEG,YBEG,XEND,YEND,NUMPB
NUMP=NUMPB+1
DELX=(XEND-XBEG)/NUMP
DELY=(YEND-YBEG)/NUMP
IF (NUMPB.GT.0) NUMP=NUMP+1
IF (DELX**2+DELY**2.EQ.0.) NUMP=1
!
DO 890 NI=1,NUMP
XP=XBEG+(NI-1)*DELX
YP=YBEG+(NI-1)*DELY
!
UX=0.
UY=0.
SIGXX=PXX
SIGYY=PYY
SIGXY=PXY
!
DO 880 J=1,NUMBE
JN=2*J
JS=JN-1
CALL INITL
XJ=XM(J)
YJ=YM(J)
AJ=A(J)
!
IF (SQRT((XP-XJ)**2+(YP-YJ)**2).LT.2.*AJ) GO TO 890
!
COSBJ=COSBET(J)
SINBJ=SINBET(J)
CALL COEFF(XP,YP,XJ,YJ,AJ,COSBJ,SINBJ,+1)
GO TO (840,810,820,830),KSYM
!
810 XJ=2.*XSYM-XM(J)
CALL COEFF(XP,YP,XJ,YJ,AJ,COSBJ,-SINBJ,-1)
GO TO 840
!
820 YJ=2.*YSYM-YM(J)
CALL COEFF(XP,YP,XJ,YJ,AJ,-COSBJ,SINBJ,-1)
GO TO 840
!
830 XJ=2.*XSYM-XM(J)
CALL COEFF(XP,YP,XJ,YJ,AJ,COSBJ,-SINBJ,-1)
XJ=XM(J)
YJ=2.*YSYM-YM(J)
CALL COEFF(XP,YP,XJ,YJ,AJ,-COSBJ,SINBJ,-1)
XJ=2.*XSYM-XM(J)
CALL COEFF(XP,YP,XJ,YJ,AJ,-COSBJ,-SINBJ,+1)
!
840 CONTINUE
!
UX=UX+UXS*P(JS)+UXN*P(JN)
UY=UY+UYS*P(JS)+UYN*P(JN)
SIGXX=SIGXX+SXXS*P(JS)+SXXN*P(JN)
SIGYY=SIGYY+SYYS*P(JS)+SYYN*P(JN)
SIGXY=SIGXY+SXYS*P(JS)+SXYN*P(JN)
!
880 CONTINUE
!
NPOINT=NPOINT+1
WRITE (6,20) NPOINT,XP,YP,UX,UY,SIGXX,SIGYY,SIGXY
!
890 CONTINUE
!
900 CONTINUE
!
! FORMAT STATEMENTS.
!
1 FORMAT (20A4)
2 FORMAT (1H1,/,25X,20A4,/)
3 FORMAT (3I4)
4 FORMAT (F6.2,E11.4,2F12.4)
5 FORMAT (3E11.4)
6 FORMAT (/,109H NUMBE_R OF STRAIGHT-LINE SEGMENTS (EACH COUNTING A&
T LEAST ONE BOUNDARY ELEMENT) USED T_O DEFINE BOUNDARIES =,I3,//,12&
6H NUMBE_R OF STRAIGHT-LINE SEGMENTS USED T_O SPECIFY OTHER LOCATION&
S (IE, NO_T ON A BOUNDARY) WHER_E RESULTS ARE T_O BE FOUND =,I3)
7 FORMAT (/,31H NO SYMMETRY CONDITIONS IMPOSED)
8 FORMAT (/,18H THE LINE X = XS =, F12.4,23H I_S A LINE OF SYMMETRY)
9 FORMAT (/,18H THE LINE Y = YS =, F12.4,23H I_S A LINE OF SYMMETRY)
10 FORMAT (/,19H THE LINES X = XS =,F12.4,13H AND Y = YS =, F12.4,&
22H ARE LINES OF SYMMETRY)
11 FORMAT (/,17H POISSONS RATIO =,F6.2,//,17H YOUNGS MODULUS =,E11.&
4)
12 FORMAT (/,31H XX-COMPONENT OF FIELD STRESS =,E11.4,//,31H YY-COMPO&
NENT OF FIELD STRESS =,E11.4,//,31H XY-COMPONENT OF FIELD STRESS =&
,E11.4)
13 FORMAT (1H1,/,27H BOUNDARY ELEMENT DAT_A,//,96H ELEMENT K&
ODE X (CENTER) Y (CENTER) LENGTH ANGLE US OR SIGMA-S&
UN OR SIGMA-N,/)
14 FORMAT (I4,4F12.4,I4,2E11.4)
15 FORMAT (2I9,3F12.4,F12.2,2E15.4)
16 FORMAT (1H1,/,65H DISPLACEMENTS AND STRESSES AT MIDPOINTS OF B&
OUNDARY ELEMENTS,//,40H ELEMENT UX UY US,&
60H UN SIGXX SIGYY SIGXY SIGMA-S SIGMA-N,&
10H SIGMA-T,/)
17 FORMAT (I10,4F10.6,6F10.1)
18 FORMAT (1H1,/,64H DISPLACEMENTS AND STRESSES AT SPECIFIED POIN&
TS I_N THE BODY,//,93H POINT X CO-ORD Y CO-ORD U&
X UY SIGXX SIGYY SIGXY,/)
19 FORMAT (4F12.4,I4)
20 FORMAT (I9,2F12.4,2F12.6,3F12.1)
!
END
!
SUBROUTINE INITL
!
COMMON/S2/SXXS,SXXN,SYYS,SYYN,SXYS,SXYN,UXS,UXN,UYS,UYN
!
SXXS=0.
SXXN=0.
SYYS=0.
SYYN=0.
SXYS=0.
SXYN=0.
!
UXS=0.
UXN=0.
UYS=0.
UYN=0.
!
RETURN
END
!
SUBROUTINE COEFF(X,Y,CX,CY,A,COSB,SINB,MSYM)
!
COMMON/S1/PI,PR,PR1,PR2,PR3,CON,COND
COMMON/S2/SXXS,SXXN,SYYS,SYYN,SXYS,SXYN,UXS,UXN,UYS,UYN
!
COS2B=COSB*COSB-SINB*SINB
SIN2B=2.*SINB*COSB
!
XB=(X-CX)*COSB+(Y-CY)*SINB
YB=-(X-CX)*SINB+(Y-CY)*COSB
!
R1S=(XB-A)*(XB-A)+YB*YB
R2S=(XB+A)*(XB+A)+YB*YB
FL1=0.5*ALOG(R1S)
FL2=0.5*ALOG(R2S)
FB2=CON*(FL1-FL2)
IF (YB.NE.0.) GO TO 10
FB3=0.
IF (ABS(XB).LT.A) FB3=CON*PI
GO TO 20
10 FB3=-CON*(ATAN((XB+A)/YB)-ATAN((XB-A)/YB))
20 FB1=YB*FB3+CON*((XB-A)*FL1-(XB+A)*FL2)
FB4=CON*(YB/R1S-YB/R2S)
FB5=CON*((XB-A)/R1S-(XB+A)/R2S)
!
UXPS=COND*(PR3*COSB*FB1+YB*(SINB*FB2+COSB*FB3))
UXPN=COND*(-PR3*SINB*FB1-YB*(COSB*FB2-SINB*FB3))
UYPS=COND*(PR3*SINB*FB1-YB*(COSB*FB2-SINB*FB3))
UYPN=COND*(PR3*COSB*FB1-YB*(SINB*FB2+COSB*FB3))
!
SXXPS=FB2+PR2*(COS2B*FB2-SIN2B*FB3)+YB*(COS2B*FB4+SIN2B*FB5)
SXXPN=FB3-PR1*(SIN2B*FB2+COS2B*FB3)+YB*(SIN2B*FB4-COS2B*FB5)
SYYPS=FB2-PR2*(COS2B*FB2-SIN2B*FB3)-YB*(COS2B*FB4+SIN2B*FB5)
SYYPN=FB3+PR1*(SIN2B*FB2+COS2B*FB3)-YB*(SIN2B*FB4-COS2B*FB5)
SXYPS=PR2*(SIN2B*FB2+COS2B*FB3)+YB*(SIN2B*FB4-COS2B*FB5)
SXYPN=PR1*(COS2B*FB2-SIN2B*FB3)-YB*(COS2B*FB4+SIN2B*FB5)
!
UXS=UXS+MSYM*UXPS
UXN=UXN+UXPN
UYS=UYS+MSYM*UYPS
UYN=UYN+UYPN
!
SXXS=SXXS+MSYM*SXXPS
SXXN=SXXN+SXXPN
SYYS=SYYS+MSYM*SYYPS
SYYN=SYYN+SYYPN
SXYS=SXYS+MSYM*SXYPS
SXYN=SXYN+SXYPN
!
RETURN
END
!
SUBROUTINE SOLVE(N)
!
COMMON/S3/A(100,100),B(100),X(100)
!
NB=N-1
DO 20 J=1,NB
L=J+1
DO 20 JJ=L,N
XM=A(JJ,J)/A(J,J)
DO 10 I=J,N
10 A(JJ,I)=A(JJ,I)-A(J,I)*XM
20 B(JJ)=B(JJ)-B(J)*XM
!
X(N)=B(N)/A(N,N)
DO 40 J=1,NB
JJ=N-J
L=JJ+1
SUMM=0.
DO 30 I=L,N
30 SUMM=SUMM+A(JJ,I)*X(I)
40 X(JJ)=(B(JJ)-SUMM)/A(JJ,JJ)
!
RETURN
END
第 94 行是：X[N-1] = B[N-1] / A[N-1][N-1];
输入数据：
6,2,4
0.2, 7.0000E+04,0,0
1.0000E+02,0,0
1,1.0,0.0,0.9659,0.2588,1,0.0,0.0
1,0.9659,0.2588,0.8660,0.5000,1,0.0,0.0
1,0.8660,0.5000,0.7071,0.7071,1,0.0,0.0
1,0.7071,0.7071,0.5000,0.8660,1,0.0,0.0
1,0.5000,0.8660,0.2588,0.9659,1,0.0,0.0
1,0.2588,0.9659,0.0,1.0,1,0.0,0.0
1.0,0.0,4.0,0.0,5
0.0,1.0,0.0,4.0,5
我在等：
NUMBE_R OF STRAIGHT-LINE SEGMENTS (EACH COUNTING AT LEAST ONE BOUNDARY ELEMENT) USED T_O DEFINE BOUNDARIES = 6
NUMBE_R OF STRAIGHT-LINE SEGMENTS USED T_O SPECIFY OTHER LOCATIONS (IE, NO_T ON A BOUNDARY) WHER_E RESULTS ARE T_O BE FOUND = 2
THE LINES X = XS = 0.0000 AND Y = YS = 0.0000 ARE LINES OF SYMMETRY
POISSONS RATIO = 0.20
YOUNGS MODULUS = 0.7000E+05
XX-COMPONENT OF FIELD STRESS = 0.1000E+03
YY-COMPONENT OF FIELD STRESS = 0.0000E+00
XY-COMPONENT OF FIELD STRESS = 0.0000E+00
BOUNDARY ELEMENT DAT_A
ELEMENT KODE X (CENTER) Y (CENTER) LENGTH ANGLE US OR SIGMA-SUN OR SIGMA-N
1 1 0.9829 0.1294 0.2610 97.51 0.0000E+00 0.0000E+00
2 1 0.9160 0.3794 0.2611 112.50 0.0000E+00 0.0000E+00
3 1 0.7865 0.6035 0.2610 127.50 0.0000E+00 0.0000E+00
4 1 0.6035 0.7865 0.2610 142.50 0.0000E+00 0.0000E+00
5 1 0.3794 0.9160 0.2611 157.50 0.0000E+00 0.0000E+00
6 1 0.1294 0.9829 0.2610 172.49 0.0000E+00 0.0000E+00
DISPLACEMENTS AND STRESSES AT MIDPOINTS OF BOUNDARY ELEMENTS
ELEMENT UX UY US UN SIGXX SIGYY SIGXY SIGMA-S SIGMA-N SIGMA-T
1 0.002847 -0.000149 -0.000519 -0.002803 -1.6 -91.5 12.1 0.0 -0.0 -93.1
2 0.002662 -0.000433 -0.001418 -0.002294 -6.1 -35.3 14.6 0.0 -0.0 -41.4
3 0.002299 -0.000678 -0.001937 -0.001411 17.9 30.4 -23.3 0.0 0.0 48.3
4 0.001776 -0.000868 -0.001937 -0.000392 95.5 56.2 -73.3 -0.0 0.0 151.7
5 0.001123 -0.000995 -0.001418 0.000490 206.0 35.3 -85.3 -0.0 -0.0 241.3
6 0.000384 -0.001059 -0.000519 0.000999 288.1 5.0 -38.0 0.0 0.0 293.1
DISPLACEMENTS AND STRESSES AT SPECIFIED POINTS I_N THE BODY
POINT X CO-ORD Y CO-ORD UX UY SIGXX SIGYY SIGXY
1 1.5000 0.0000 0.002246 -0.000000 14.8 -8.0 -0.0
2 2.0000 0.0000 0.001772 -0.000000 44.4 3.2 -0.0
3 2.5000 0.0000 0.001451 -0.000000 62.1 4.3 0.0
4 3.0000 0.0000 0.001223 -0.000000 72.9 3.9 -0.0
5 3.5000 0.0000 0.001056 -0.000000 79.7 3.2 -0.0
6 4.0000 0.0000 0.000929 -0.000000 84.3 2.7 -0.0
7 0.0000 1.5000 0.000000 -0.001050 154.5 38.6 -0.0
8 0.0000 2.0000 -0.000000 -0.000876 123.0 29.4 -0.0
9 0.0000 2.5000 -0.000000 -0.000733 112.4 21.1 -0.0
10 0.0000 3.0000 -0.000000 -0.000626 107.8 15.5 0.0
11 0.0000 3.5000 -0.000000 -0.000544 105.3 11.8 -0.0
12 0.0000 4.0000 -0.000000 -0.000480 103.9 9.2 -0.0```
好的，这里有一个 C++ 版本，可以让一切顺利进行。它提供与 Fortran 版本相同的数字输出。有一定的重复，所以可以改进。
我在处理基于 0 的 (C++) 数组与基于 1 的 (Fortran) 数组时，只是采取了懦夫的方式，使数组比所需的元素大一个元素。每个数组中的 Element[0] 未使用。修复此问题可能会带来输入文件不适用于 Fortran 和 C++ 版本的风险。我还使用了固定大小的数组（尽管使用向量）；这需要改进，但输入文件的格式并不容易。如有必要，请添加检查以确保不超出数组边界。
Fortran 代码非常古老（我认为是 Fortran 66）。该语言中不再存在 Hollerith 常量和计算的 goto。首先将其转换为现代 Fortran 可能会更好 - 世界在前进。我什至编译 Fortran 样本都遇到了很大的困难。我要么根除计算转到部分，要么将它们更改为 C++ 中的 switch 语句。
我仅通过删除标题行并将逗号更改为空格来修改输入文件；对于任何语言来说，尝试读取固定字段宽度都是一场噩梦。 请使用代码下方的输入文件。
C++代码
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <vector>
#include <cmath>
using namespace std;
//---------- Some useful derived types -----------
struct Stress
{
double xx;
double yy;
double xy;
Stress( double a = 0.0, double b = 0.0, double c = 0.0 ) : xx(a), yy(b), xy(c) {};
};
Stress operator + ( const Stress &a, const Stress &b ){ return Stress( a.xx + b.xx, a.yy + b.yy, a.xy + b.xy ); }
Stress operator * ( double a, const Stress &b ){ return Stress( a * b.xx, a * b.yy, a * b.xy ); }
Stress operator * ( const Stress &a, double b ){ return b * a; }
istream &operator >> ( istream &in, Stress &a ){ return in >> a.xx >> a.yy >> a.xy; }
ostream &operator << ( ostream &out, const Stress &a ){ return out << a.xx << ", " << a.yy << ", " << a.xy; }
struct Point{
double x;
double y;
Point( double a = 0.0, double b = 0.0 ) : x(a), y(b) {};
};
Point operator + ( const Point &a, const Point &b ){ return Point( a.x + b.x, a.y + b.y ); }
Point operator * ( double a, const Point &b ){ return Point( a * b.x, a * b.y ); }
Point operator * ( const Point &a, double b ){ return b * a; }
struct TestCase
{
double pr, pr1, pr2, pr3;
double con, cond;
};
//---------- Some constants ----------------------
const double PI = 4.0 * atan( 1.0 );
const int NMAX = 50;
//---------- Function declarations ---------------
void coeff( double x, double y, double cx, double cy, double a, double cosb,double sinb, int msym,
const TestCase &T, Stress &Ss, Stress &Sn, Point &Us, Point &Un );
vector<double> solve( int n, vector<vector<double>> A, vector<double> b );
//------------------------------------------------
int main()
{
string title;
int numbs, numos, numpb;
int ksym, kode;
int num;
double PR, E;
double xsym, ysym;
double xbeg, ybeg, xend, yend;
double bvs, bvn;
Stress Pf;
vector<vector<double>> c( 1 + 2 * NMAX, vector<double>( 1 + 2 * NMAX ) );
vector<double> b( 1 + 2 * NMAX );
vector<double> xm( 1 + NMAX ), ym( 1 + NMAX ), a( 1 + NMAX ), cosbet( 1 + NMAX ), sinbet( 1 + NMAX );
vector<int> kod( 1 + NMAX );
// Input
ifstream in( "input.txt" );
// getline( in, title ); // Stopped reading the title; it isn't in the input file
in >> numbs >> numos >> ksym;
in >> PR >> E >> xsym >> ysym;
in >> Pf;
cout << title << '\n';
cout << "Number of straight-line segments (each counting at least one boundary element) used to define boundaries = " << numbs << '\n';
cout << "Number of straight-line segments used to specify other locations (ie, not on a boundary) where results are to be found = " << numos << '\n';
if ( ksym == 2 || ksym == 4 ) cout << "The line x = " << xsym << " is a line of symmetry\n";
if ( ksym == 3 || ksym == 4 ) cout << "The line y = " << ysym << " is a line of symmetry\n";
cout << "Poisson's ratio = " << PR << " Young's modulus = " << E << '\n';
cout << "Field stress: xx, yy, xy = " << Pf << '\n';
TestCase T;
T.con = 1 / ( 4 * PI * ( 1 - PR ) );
T.cond = ( 1 + PR ) / E;
T.pr1 = 1 - 2 * PR;
T.pr2 = 2 * ( 1 - PR );
T.pr3 = 3 - 4 * PR;
// Define locations, sizes, orientations and boundary conditions of boundary elements
int numbe = 0;
for ( int n = 1; n <= numbs; n++ )
{
in >> num >> xbeg >> ybeg >> xend >> yend >> kode >> bvs >> bvn;
double xd = ( xend - xbeg ) / num;
double yd = ( yend - ybeg ) / num;
double sw = sqrt( xd * xd + yd * yd );
for ( int ne = 1; ne <= num; ne++ )
{
numbe++;
int m = numbe;
xm[m] = xbeg + 0.5 * ( 2 * ne - 1 ) * xd;
ym[m] = ybeg + 0.5 * ( 2 * ne - 1 ) * yd;
a[m] = 0.5 * sw;
sinbet[m] = yd / sw;
cosbet[m] = xd / sw;
kod[m] = kode;
int mn = 2 * m;
int ms = mn - 1;
b[ms] = bvs;
b[mn] = bvn;
}
}
#define FMT << " " << setw( 15 ) << setprecision( 4 ) << scientific <<
cout << "Boundary element data\n";
cout FMT "element" FMT "kode" FMT "x(center)" FMT "y(center)" FMT "length" FMT "angle" FMT "us or sigma-s" FMT "un or sigma-n" << '\n';
for ( int m = 1; m <= numbe; m++ )
{
double size = 2 * a[m];
double angle = 180 * atan2( sinbet[m], cosbet[m] ) / PI;
cout FMT m FMT kod[m] FMT xm[m] FMT ym[m] FMT size FMT angle FMT b[2*m-1] FMT b[2*m] << '\n';
}
// Adjust stress boundary values to account for initial stresses.
for ( int n = 1; n <= numbe; n++ )
{
int nn = 2 * n;
int ns = nn - 1;
double cosb = cosbet[n];
double sinb = sinbet[n];
double sigmas = ( Pf.yy - Pf.xx ) * sinb * cosb + Pf.xy * ( cosb * cosb - sinb * sinb );
double sigman = Pf.xx * sinb * sinb - 2 * Pf.xy * sinb * cosb + Pf.yy * cosb * cosb;
if ( kod[n] == 1 || kod[n] == 4 ) b[ns] -= sigmas;
if ( kod[n] == 1 || kod[n] == 3 ) b[nn] -= sigman;
}
// Compute influence coefficients and set up system of algebraic equations.
for ( int i = 1; i <= numbe; i++ )
{
int in = 2 * i;
int is = in - 1;
double xi = xm[i], yi = ym[i], cosbi = cosbet[i], sinbi = sinbet[i];
int kode = kod[i];
for ( int j = 1; j <= numbe; j++ )
{
int jn = 2 * j;
int js = jn - 1;
double xj = xm[j], yj = ym[j], cosbj = cosbet[j], sinbj = sinbet[j], aj = a[j];
Stress Ss, Sn;
Point Us, Un;
coeff( xi, yi, xj, yj, aj, cosbj, sinbj, +1, T, Ss, Sn, Us, Un );
switch( ksym )
{
case 1:
break;
case 2:
xj = 2 * xsym - xm[j];
coeff( xi, yi, xj, yj, aj, cosbj, -sinbj, -1, T, Ss, Sn, Us, Un );
break;
case 3:
yj = 2 * ysym - ym[j];
coeff( xi, yi, xj, yj, aj, -cosbj, sinbj, -1, T, Ss, Sn, Us, Un );
break;
case 4:
xj = 2 * xsym - xm[j];
coeff( xi, yi, xj, yj, aj, cosbj, -sinbj, -1, T, Ss, Sn, Us, Un );
xj = xm[j];
yj = 2 * ysym - ym[j];
coeff( xi, yi, xj, yj, aj, -cosbj, sinbj, -1, T, Ss, Sn, Us, Un );
xj = 2 * xsym - xm[j];
coeff( xi, yi, xj, yj, aj, -cosbj, -sinbj, +1, T, Ss, Sn, Us, Un );
break;
}
switch( kode )
{
case 1:
c[is][js] = ( Ss.yy - Ss.xx ) * sinbi * cosbi + Ss.xy * ( cosbi * cosbi - sinbi * sinbi );
c[is][jn] = ( Sn.yy - Sn.xx ) * sinbi * cosbi + Sn.xy * ( cosbi * cosbi - sinbi * sinbi );
c[in][js] = Ss.xx * sinbi * sinbi - 2 * Ss.xy * sinbi * cosbi + Ss.yy * cosbi * cosbi;
c[in][jn] = Sn.xx * sinbi * sinbi - 2 * Sn.xy * sinbi * cosbi + Sn.yy * cosbi * cosbi;
break;
case 2:
c[is][js] = Us.x * cosbi + Us.y * sinbi;
c[is][jn] = Un.x * cosbi + Un.y * sinbi;
c[in][js] = -Us.x * sinbi + Us.y * cosbi;
c[in][jn] = -Un.x * sinbi + Un.y * cosbi;
break;
case 3:
c[is][js] = Us.x * cosbi + Us.y * sinbi;
c[is][jn] = Un.x * cosbi + Un.y * sinbi;
c[in][js] = Ss.xx * sinbi * sinbi - 2 * Ss.xy * sinbi * cosbi + Ss.yy * cosbi * cosbi;
c[in][jn] = Sn.xx * sinbi * sinbi - 2 * Sn.xy * sinbi * cosbi + Sn.yy * cosbi * cosbi;
break;
case 4:
c[is][js] = ( Ss.yy - Ss.xx ) * sinbi * cosbi + Ss.xy * ( cosbi * cosbi - sinbi * sinbi );
c[is][jn] = ( Sn.yy - Sn.xx ) * sinbi * cosbi + Sn.xy * ( cosbi * cosbi - sinbi * sinbi );
c[in][js] = -Us.x * sinbi + Us.y * cosbi;
c[in][jn] = -Un.x * sinbi + Un.y * cosbi;
break;
}
}
}
// Solve system of algebraic equations.
vector<double> p = solve( 2 * numbe, c, b );
// Compute boundary displacements and stresses.
cout << "Displacements and stresses at midpoints of boundary elements\n";
cout FMT "element" FMT "ux" FMT "uy" FMT "us" FMT "un" FMT "sigxx" FMT "sigyy" FMT "sigxy" FMT "sigma-s" FMT "sigma-n" FMT "sigma-t" << '\n';
for ( int i = 1; i <= numbe; i++ )
{
double xi = xm[i], yi = ym[i], cosbi = cosbet[i], sinbi = sinbet[i];
Stress Sig = Pf;
Point U;
for ( int j = 1; j <= numbe; j++ )
{
int jn = 2 * j;
int js = jn - 1;
double xj = xm[j], yj = ym[j], aj = a[j], cosbj = cosbet[j], sinbj = sinbet[j];
Point Us, Un;
Stress Ss, Sn;
coeff( xi, yi, xj, yj, aj, cosbj, sinbj, + 1, T, Ss, Sn, Us, Un );
switch( ksym )
{
case 1:
break;
case 2:
xj = 2 * xsym - xm[j];
coeff( xi, yi, xj, yj, aj, cosbj, -sinbj, -1, T, Ss, Sn, Us, Un );
break;
case 3:
yj = 2 * ysym - ym[j];
coeff( xi, yi, xj, yj, aj, -cosbj, sinbj, -1, T, Ss, Sn, Us, Un );
break;
case 4:
xj = 2 * xsym - xm[j];
coeff( xi, yi, xj, yj, aj, cosbj, -sinbj, -1, T, Ss, Sn, Us, Un );
xj = xm[j];
yj = 2 * ysym - ym[j];
coeff( xi, yi, xj, yj, aj, -cosbj, sinbj, -1, T, Ss, Sn, Us, Un );
xj = 2 * xsym - xm[j];
coeff( xi, yi, xj, yj, aj, -cosbj, -sinbj, +1, T, Ss, Sn, Us, Un );
break;
}
U = U + p[js] * Us + p[jn] * Un;
Sig = Sig + p[js] * Ss + p[jn] * Sn;
}
double us = U.x * cosbi + U.y * sinbi;
double un = -U.x * sinbi + U.y * cosbi;
double sigmas = ( Sig.yy - Sig.xx ) * sinbi * cosbi + Sig.xy * ( cosbi * cosbi-sinbi * sinbi );
double sigman = Sig.xx * sinbi * sinbi - 2 * Sig.xy * sinbi * cosbi + Sig.yy * cosbi * cosbi;
double sigmat = Sig.xx * cosbi * cosbi + 2 * Sig.xy * cosbi * sinbi + Sig.yy * sinbi * sinbi;
cout FMT i FMT U.x FMT U.y FMT us FMT un FMT Sig.xx FMT Sig.yy FMT Sig.xy FMT sigmas FMT sigman FMT sigmat << '\n';
}
// Compute displacements and stresses at specified points in body
if ( numos > 0 )
{
cout << "Displacements and stresses at specified points in the body\n";
cout FMT "point" FMT "x co-ord" FMT "y co-ord" FMT "u" FMT "uy" FMT "sigxx" FMT "sigyy" FMT "sigxy" << '\n';
int npoint = 0;
for ( int n = 1; n <= numos; n++ )
{
in >> xbeg >> ybeg >> xend >> yend >> numpb;
int nump = numpb + 1;
double delx = ( xend - xbeg ) / nump;
double dely = ( yend - ybeg ) / nump;
if ( numpb > 0 ) nump++;
if ( delx * delx + dely * dely == 0.0 ) nump = 1;
for ( int ni = 1; ni <= nump; ni++ )
{
double xp = xbeg + ( ni - 1 ) * delx;
double yp = ybeg + ( ni - 1 ) * dely;
Point U;
Stress Sig = Pf;
bool ok = true;
for ( int j = 1; j <= numbe; j++ )
{
int jn = 2 * j;
int js = jn - 1;
double xj = xm[j], yj = ym[j], aj = a[j];
double cosbj = cosbet[j], sinbj = sinbet[j];
if ( sqrt( ( xp - xj ) * ( xp - xj ) + ( yp - yj ) * ( yp - yj ) ) < 2 * aj )
{
ok = false;
break;
}
Point Us, Un;
Stress Ss, Sn;
coeff( xp, yp, xj, yj, aj, cosbj, sinbj, +1, T, Ss, Sn, Us, Un );
switch( ksym )
{
case 1:
break;
case 2:
xj = 2 * xsym - xm[j];
coeff( xp, yp, xj, yj, aj, cosbj, -sinbj, -1, T, Ss, Sn, Us, Un );
break;
case 3:
yj = 2 * ysym - ym[j];
coeff( xp, yp, xj, yj, aj, -cosbj, sinbj, -1, T, Ss, Sn, Us, Un );
break;
case 4:
xj = 2 * xsym - xm[j];
coeff( xp, yp, xj, yj, aj, cosbj, -sinbj, -1, T, Ss, Sn, Us, Un );
xj = xm[j];
yj = 2 * ysym - ym[j];
coeff( xp, yp, xj, yj, aj, -cosbj, sinbj, -1, T, Ss, Sn, Us, Un );
xj = 2 * xsym - xm[j];
coeff( xp, yp, xj, yj, aj, -cosbj, -sinbj, +1, T, Ss, Sn, Us, Un );
break;
}
U = U + p[js] * Us + p[jn] * Un;
Sig = Sig + p[js] * Ss + p[jn] * Sn;
}
if ( ok )
{
npoint++;
cout FMT npoint FMT xp FMT yp FMT U.x FMT U.y FMT Sig.xx FMT Sig.yy FMT Sig.xy << '\n';
}
}
}
}
}
//------------------------------------------------
void coeff( double x, double y, double cx, double cy, double a, double cosb,double sinb, int msym,
const TestCase &T, Stress &Ss, Stress &Sn, Point &Us, Point &Un )
{
double cos2b = cosb * cosb - sinb * sinb;
double sin2b = 2 * sinb * cosb;
double xb = ( x - cx ) * cosb + ( y - cy ) * sinb;
double yb = -( x - cx ) * sinb + ( y - cy ) * cosb;
double r1s = ( xb - a ) * ( xb - a ) + yb * yb;
double r2s = ( xb + a ) * ( xb + a ) + yb * yb;
double fl1 = 0.5 * log( r1s );
double fl2 = 0.5 * log( r2s );
double fb2 = T.con * ( fl1 - fl2 );
double fb3 = 0.0;
if ( yb == 0.0 )
{
if ( abs( xb ) < a ) fb3 = T.con * PI;
}
else
{
fb3 = -T.con * ( atan( ( xb + a ) / yb ) - atan( ( xb - a ) / yb ) );
}
double fb1 = yb * fb3 + T.con * ( ( xb - a ) * fl1 - ( xb + a ) * fl2 );
double fb4 = T.con * ( yb / r1s - yb / r2s );
double fb5 = T.con * ( ( xb - a ) / r1s - ( xb + a ) / r2s );
Point Ups, Upn;
Ups.x = T.cond * ( T.pr3 * cosb * fb1 + yb * ( sinb * fb2 + cosb * fb3 ) );
Upn.x = T.cond * ( -T.pr3 * sinb * fb1 - yb * ( cosb * fb2 - sinb * fb3 ) );
Ups.y = T.cond * ( T.pr3 * sinb * fb1 - yb * ( cosb * fb2 - sinb * fb3 ) );
Upn.y = T.cond * ( T.pr3 * cosb * fb1 - yb * ( sinb * fb2 + cosb * fb3 ) );
Stress Sps, Spn;
Sps.xx = fb2 + T.pr2 * ( cos2b * fb2 - sin2b * fb3 ) + yb * ( cos2b * fb4 + sin2b * fb5 );
Spn.xx = fb3 - T.pr1 * ( sin2b * fb2 + cos2b * fb3 ) + yb * ( sin2b * fb4 - cos2b * fb5 );
Sps.yy = fb2 - T.pr2 * ( cos2b * fb2 - sin2b * fb3 ) - yb * ( cos2b * fb4 + sin2b * fb5 );
Spn.yy = fb3 + T.pr1 * ( sin2b * fb2 + cos2b * fb3 ) - yb * ( sin2b * fb4 - cos2b * fb5 );
Sps.xy = T.pr2 * ( sin2b * fb2 + cos2b * fb3 ) + yb * ( sin2b * fb4 - cos2b * fb5 );
Spn.xy = T.pr1 * ( cos2b * fb2 - sin2b * fb3 ) - yb * ( cos2b * fb4 + sin2b * fb5 );
Us = Us + msym * Ups;
Un = Un + Upn;
Ss = Ss + msym * Sps;
Sn = Sn + Spn;
}
//------------------------------------------------
vector<double> solve( int n, vector<vector<double>> A, vector<double> b )
{
for ( int j = 1; j <= n - 1; j++ )
{
for ( int jj = j + 1; jj <= n; jj++ )
{
double xm = A[jj][j] / A[j][j];
for ( int i = j; i <= n; i++ ) A[jj][i] -= A[j][i] * xm;
b[jj] -= b[j] * xm;
}
}
vector<double> x( b.size() );
x[n] = b[n] / A[n][n];
for ( int j = 1; j <= n - 1; j++ )
{
int jj = n - j;
double summ = 0.0;
for ( int i = jj + 1; i <= n; i++ ) summ += A[jj][i] * x[i];
x[jj] = ( b[jj] - summ ) / A[jj][jj];
}
return x;
}
输入文件（input.txt）
6 2 4
0.2 7.0000E+04 0 0
1.0000E+02 0 0
1 1.0 0.0 0.9659 0.2588 1 0.0 0.0
1 0.9659 0.2588 0.8660 0.5000 1 0.0 0.0
1 0.8660 0.5000 0.7071 0.7071 1 0.0 0.0
1 0.7071 0.7071 0.5000 0.8660 1 0.0 0.0
1 0.5000 0.8660 0.2588 0.9659 1 0.0 0.0
1 0.2588 0.9659 0.0 1.0 1 0.0 0.0
1.0 0.0 4.0 0.0 5
0.0 1.0 0.0 4.0 5
输出（屏幕截图）
Number of straight-line segments (each counting at least one boundary element) used to define boundaries = 6
Number of straight-line segments used to specify other locations (ie, not on a boundary) where results are to be found = 2
The line x = 0 is a line of symmetry
The line y = 0 is a line of symmetry
Poisson's ratio = 0.2 Young's modulus = 70000
Field stress: xx, yy, xy = 100, 0, 0
Boundary element data
element kode x(center) y(center) length angle us or sigma-s un or sigma-n
1 1 9.8295e-01 1.2940e-01 2.6104e-01 9.7506e+01 0.0000e+00 0.0000e+00
2 1 9.1595e-01 3.7940e-01 2.6107e-01 1.1250e+02 0.0000e+00 0.0000e+00
3 1 7.8655e-01 6.0355e-01 2.6104e-01 1.2750e+02 0.0000e+00 0.0000e+00
4 1 6.0355e-01 7.8655e-01 2.6104e-01 1.4250e+02 0.0000e+00 0.0000e+00
5 1 3.7940e-01 9.1595e-01 2.6107e-01 1.5750e+02 0.0000e+00 0.0000e+00
6 1 1.2940e-01 9.8295e-01 2.6104e-01 1.7249e+02 0.0000e+00 0.0000e+00
Displacements and stresses at midpoints of boundary elements
element ux uy us un sigxx sigyy sigxy sigma-s sigma-n sigma-t
1 2.8470e-03 -1.4881e-04 -5.1944e-04 -2.8031e-03 -1.5886e+00 -9.1501e+01 1.2056e+01 5.3291e-15 2.0428e-14 -9.3089e+01
2 2.6619e-03 -4.3250e-04 -1.4182e-03 -2.2938e-03 -6.0549e+00 -3.5296e+01 1.4619e+01 -1.7764e-15 -2.6645e-15 -4.1351e+01
3 2.2991e-03 -6.7780e-04 -1.9373e-03 -1.4114e-03 1.7882e+01 3.0376e+01 -2.3307e+01 8.8818e-15 1.0658e-14 4.8258e+01
4 1.7759e-03 -8.6797e-04 -1.9373e-03 -3.9241e-04 9.5517e+01 5.6230e+01 -7.3287e+01 7.1054e-15 2.1316e-14 1.5175e+02
5 1.1228e-03 -9.9532e-04 -1.4182e-03 4.8992e-04 2.0600e+02 3.5338e+01 -8.5321e+01 -1.4211e-14 1.0658e-14 2.4134e+02
6 3.8430e-04 -1.0586e-03 -5.1929e-04 9.9930e-04 2.8809e+02 5.0017e+00 -3.7960e+01 0.0000e+00 -1.7764e-15 2.9310e+02
Displacements and stresses at specified points in the body
point x co-ord y co-ord u uy sigxx sigyy sigxy
1 1.5000e+00 0.0000e+00 2.2461e-03 -5.0975e-20 1.4850e+01 -7.9906e+00 1.2179e-15
2 2.0000e+00 0.0000e+00 1.7725e-03 3.6011e-20 4.4416e+01 3.1933e+00 1.8014e-16
3 2.5000e+00 0.0000e+00 1.4505e-03 -6.2163e-20 6.2148e+01 4.3214e+00 -2.1320e-16
4 3.0000e+00 0.0000e+00 1.2235e-03 -1.0520e-19 7.2855e+01 3.8602e+00 -6.0026e-17
5 3.5000e+00 0.0000e+00 1.0563e-03 1.5125e-19 7.9676e+01 3.2167e+00 -2.2186e-16
6 4.0000e+00 0.0000e+00 9.2859e-04 1.3077e-19 8.4250e+01 2.6519e+00 -3.3073e-16
7 0.0000e+00 1.5000e+00 0.0000e+00 -1.0503e-03 1.5449e+02 3.8647e+01 0.0000e+00
8 0.0000e+00 2.0000e+00 0.0000e+00 -8.7560e-04 1.2296e+02 2.9426e+01 0.0000e+00
9 0.0000e+00 2.5000e+00 0.0000e+00 -7.3302e-04 1.1242e+02 2.1110e+01 0.0000e+00
10 0.0000e+00 3.0000e+00 0.0000e+00 -6.2558e-04 1.0777e+02 1.5519e+01 0.0000e+00
11 0.0000e+00 3.5000e+00 0.0000e+00 -5.4382e-04 1.0532e+02 1.1782e+01 0.0000e+00
12 0.0000e+00 4.0000e+00 0.0000e+00 -4.8017e-04 1.0389e+02 9.2101e+00 0.0000e+00