There is code in Fortran that I am trying to translate into C++. Apparently, I misunderstood the Fortran syntax somewhere and that's why I'm getting a problem on line 94. I looked at the values of the components on line 94 and realized that, with a high probability, I incorrectly calculated the two-dimensional array A. But I double-checked its calculations many times and cannot find any errors. Please tell me what needs to be done for the code to work correctly? (I used debuger, I’m inclined to believe that I messed something up with the difference in syntax in the languages)
С++:
#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:
!
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
Line 94 is: X[N-1] = B[N-1] / A[N-1][N-1];
Input data:
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
I am waiting:
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```
OK, here's a C++ version which will get the ball rolling. It gives the same numerical output as the Fortran version. There's a certain amount of repetition, so it could be improved.
I've dealt with the 0-based (C++) versus 1-based (Fortran) arrays simply by taking the coward's way out and making the arrays one element larger than they need to be. Element[0] in each array isn't used. Fixing this runs the risk of having input files that don't work with both Fortran and C++ versions. I've also used fixed-size arrays (despite using vectors); this needs improving, but the format of the input file doesn't make that easy. Add a check that you don't exceed array bounds if necessary.
The Fortran code is very old (Fortran 66, I think). Hollerith constants and computed gotos no longer exist in the language. It might have been better to convert it to modern Fortran first - the world moves on. I had great difficulty even compiling the Fortran sample. I've either eradicated the computed-goto sections or changed them for switch statements in C++.
I have modified the input file only by removing the title line and changing the commas to spaces; trying to read fixed field widths is a nightmare in any language. Please use the input file below the code.
C++ Code
#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 file (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
Output (screen capture)
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