📄 hexafron.c
字号:
}
gn[i][j] = .125 * xi[i][j] * c;
}
}
/* --- formation of jacobian tj */
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
tj[i][j] = 0;
for (k = 0; k < 8; k++) {
kn = abs(noc[8*n+k]);
tj[i][j] = tj[i][j] + gn[i][k] * x[3*(kn-1)+j];
}
}
}
/* --- determinant of the jacobian */
dj1 = tj[0][0] * (tj[1][1] * tj[2][2] - tj[2][1] * tj[1][2]);
dj2 = tj[0][1] * (tj[1][2] * tj[2][0] - tj[2][2] * tj[1][0]);
dj3 = tj[0][2] * (tj[1][0] * tj[2][1] - tj[2][0] * tj[1][1]);
*dj = dj1 + dj2 + dj3;
/* --- inverse of the jacobian aj() */
aj[0][0] = (tj[1][1] * tj[2][2] - tj[1][2] * tj[2][1]) / *dj;
aj[0][1] = (tj[2][1] * tj[0][2] - tj[2][2] * tj[0][1]) / *dj;
aj[0][2] = (tj[0][1] * tj[1][2] - tj[0][2] * tj[1][1]) / *dj;
aj[1][0] = (tj[1][2] * tj[2][0] - tj[1][0] * tj[2][2]) / *dj;
aj[1][1] = (tj[0][0] * tj[2][2] - tj[0][2] * tj[2][0]) / *dj;
aj[1][2] = (tj[0][2] * tj[1][0] - tj[0][0] * tj[1][2]) / *dj;
aj[2][0] = (tj[1][0] * tj[2][1] - tj[1][1] * tj[2][0]) / *dj;
aj[2][1] = (tj[0][1] * tj[2][0] - tj[0][0] * tj[2][1]) / *dj;
aj[2][2] = (tj[0][0] * tj[1][1] - tj[0][1] * tj[1][0]) / *dj;
/* --- h() matrix relates local derivatives of u to local */
/* displacements q */
for (i = 0; i < 9; i++) {
for (j = 0; j < 24; j++) {
h[i][j] = 0;
}
}
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
ir = 3 * i + j ;
for (k = 0; k < 8; k++) {
ic = 3 * k + i;
h[ir][ic] = gn[j][k];
}
}
}
/* --- g() matrix relates strains to local derivatives of u */
for (i = 0; i < 6; i++) {
for (j = 0; j < 9; j++) {
g[i][j] = 0;
}
}
g[0][0] = aj[0][0];
g[0][1] = aj[0][1];
g[0][2] = aj[0][2];
g[1][3] = aj[1][0];
g[1][4] = aj[1][1];
g[1][5] = aj[1][2];
g[2][6] = aj[2][0];
g[2][7] = aj[2][1];
g[2][8] = aj[2][2];
g[3][3] = aj[2][0];
g[3][4] = aj[2][1];
g[3][5] = aj[2][2];
g[3][6] = aj[1][0];
g[3][7] = aj[1][1];
g[3][8] = aj[1][2];
g[4][0] = aj[2][0];
g[4][1] = aj[2][1];
g[4][2] = aj[2][2];
g[4][6] = aj[0][0];
g[4][7] = aj[0][1];
g[4][8] = aj[0][2];
g[5][0] = aj[1][0];
g[5][1] = aj[1][1];
g[5][2] = aj[1][2];
g[5][3] = aj[0][0];
g[5][4] = aj[0][1];
g[5][5] = aj[0][2];
/* --- b() matrix relates strains to q */
for (i = 0; i < 6; i++) {
for (j = 0; j < 24; j++) {
b[i][j] = 0;
for (k = 0; k < 9; k++) {
b[i][j] = b[i][j] + g[i][k] * h[k][j];
}
}
}
/* --- db() matrix relates stresses to q */
for (i = 0; i < 6; i++) {
for (j = 0; j < 24; j++ ) {
db[i][j] = 0;
for (k = 0; k < 6; k++) {
db[i][j] = db[i][j] + d[i][k] * b[k][j];
}
}
}
return(0);
}
prefront(nn,ne,nen,ndn,nq,nmpc,noc,mpc,ibl)
int nn,ne,nen,ndn,nq,nmpc,*ibl;
int *noc,*mpc;
{
int i,j,k,n,ifron,*ide,i1,ia,ineg;
/* ----- mark last appearance of node / make it negative in noc() */
/* last appearance is first appearance for reverse element order */
for (i = 1; i <= nn; i++) {
for (j = ne-1; j >= 0; j--) {
for (k = 0; k < nen; k++) {
if (i == noc[nen*j+k])
goto label1;
}
}
label1:
noc[nen*j+k] = -i;
}
/* ===== block size determination */
ide = (int *) calloc(nq+1, sizeof(int));
for (i = 1; i <= nq; i++) {
ide[i] = 0;
}
for (i = 0; i < nmpc; i++) {
for (j = 0; j < 2; j++) {
ide[mpc[2*i+j]] = 1;
}
}
ifron = 0;
for (i = 1; i <= nq; i++) {
ifron = ifron + ide[i];
}
*ibl = ifron;
for (n = 0; n < ne; n++) {
ineg = 0;
for (i = 0; i < nen; i++) {
i1 = noc[nen*n+i];
ia = ndn * (abs(i1) - 1);
for (j = 0; j < ndn; j++) {
ia = ia + 1;
if (ide[ia] == 0) {
ifron = ifron + 1;
ide[ia] = 1;
}
}
if (i1 < 0)
ineg = ineg + 1;
}
if (*ibl < ifron)
*ibl = ifron;
ifron = ifron - ndn * ineg;
}
printf( "block size = %d\n", *ibl);
return(0);
}
mpcfron(indx,isbl,mpc,nmpc,nfron,s,f,ibl,beta,cnst)
int *indx,*isbl,*mpc,nmpc,*nfron,ibl;
double *s,*f,*beta,cnst;
{
int i,j,i1,j1,ifl,k,k1,i2;
/* ----- modifications for multipoint constraints by penalty method */
for (i = 0; i < nmpc; i++) {
i1 = mpc[2*i];
ifl = 0;
for (j = 1; j <= *nfron; j++) {
j1 = indx[j];
if (i1 == isbl[j1]) {
ifl = 1;
break;
}
}
if (ifl == 0) {
*nfron = *nfron + 1;
j1 = indx[*nfron];
isbl[j1] = i1;
}
i2 = mpc[2*i+1];
ifl = 0;
for (k = 1; k <= *nfron; k++) {
k1 = indx[k];
if (k1 == isbl[k1]) {
ifl = 1;
break;
}
}
if (ifl == 0) {
*nfron = *nfron + 1;
k1 = indx[*nfron];
isbl[k1] = i2;
}
/* ----- stiffness modification */
j1 = j1 - 1;
k1 = k1 - 1;
s[ibl*j1+j1] = s[ibl*j1+j1] + cnst * beta[3*i] * beta[3*i];
s[ibl*k1+k1] = s[ibl*k1+k1] + cnst * beta[3*i+1] * beta[3*i+1];
s[ibl*j1+k1] = s[ibl*j1+k1] + cnst * beta[3*i] * beta[3*i+1];
s[ibl*k1+j1] = s[ibl*j1+k1];
/* ----- force modification */
f[i1-1] = f[i1-1] + cnst * beta[3*i+2] * beta[3*i];
f[i2-1] = f[i2-1] + cnst * beta[3*i+2] * beta[3*i+1];
}
return(0);
}
front(n,noc,nen,ndn,nd,icount,indx,isbl,ibl,s,f,nfron,ntogo,ndcnt,se,nu,cnst,u)
int n,nen,ndn,nd,*noc,*indx,*isbl,ibl,*nfron,*ntogo,*ndcnt,*nu,*icount;
double *s,*f,se[][24],cnst,*u;
{
struct data record;
double pivot,c;
int i,i1,ia,is1,idj,idf,ie1,ifl,ii,ix,j,j1;
int itemp,ipv,ipg,ig,iebl[24],ntg1,iba;
/* ----- frontal method assembly and elimination ----- */
/* ------------- assembly of element n ------------- */
for (i = 0; i < nen; i++) {
i1 = noc[nen*n+i];
ia = abs(i1);
is1 = 1;
if(i1 < 0)
is1 = -1;
idf = ndn * (ia - 1);
ie1 = ndn * i;
for (j = 0; j < ndn; j++) {
idf = idf + 1;
ie1 = ie1 + 1;
ifl = 0;
if (*nfron > *ntogo) {
for (ii = *ntogo+1; ii <= *nfron; ii++) {
ix = indx[ii];
if (idf == isbl[ix]) {
ifl = 1;
break;
}
}
}
if (ifl == 0) {
*nfron = *nfron + 1;
ii = *nfron;
ix = indx[ii];
}
isbl[ix] = idf;
iebl[ie1] = ix;
if (is1 == -1) {
*ntogo = *ntogo + 1;
itemp = indx[*ntogo];
indx[*ntogo] = indx[ii];
indx[ii] = itemp;
}
}
}
for (i = 0; i < 24; i++) {
i1 = iebl[i+1]-1;
for (j = 0; j < 24; j++) {
j1 = iebl[j+1]-1;
s[ibl*i1+j1] = s[ibl*i1+j1] + se[i][j];
}
}
/* ------------------------------------------------------------------ */
if (*ndcnt < nd) {
/* ----- modification for displacement bcs / penalty approach ----- */
for (i = 1; i <= *ntogo; i++) {
i1 = indx[i];
ig = isbl[i1];
for (j = 0; j < nd; j++) {
if (ig == nu[j]) {
i1 = i1 - 1;
s[ibl*i1+i1] = s[ibl*i1+i1] + cnst;
f[ig-1] = f[ig-1] + cnst * u[j];
*ndcnt = *ndcnt + 1; /* counter for check */
break;
}
}
}
}
/* ------------ elimination of completed variables --------------- */
ntg1 = *ntogo;
for (ii = 1; ii <= ntg1; ii++) {
ipv = indx[1];
ipg = isbl[ipv];
pivot = s[(ibl+1)*(ipv-1)];
/* ----- write separator "0" and pivot value to disk ----- */
*icount = *icount + 1;
record.variable = 0;
record.coefft = pivot;
fwrite(&record, sizeof(record),1,fptr);
s[(ibl+1)*(ipv-1)] = 0;
for (i = 2; i <= *nfron; i++) {
i1 = indx[i];
ig = isbl[i1];
if (s[ibl*(i1-1)+ipv-1] != 0) {
c = s[ibl*(i1-1)+ipv-1] / pivot;
s[ibl*(i1-1)+ipv-1] = 0;
for (j = 2; j <= *nfron; j++) {
j1 = indx[j];
if (s[ibl*(ipv-1)+j1-1] != 0)
s[ibl*(i1-1)+j1-1] = s[ibl*(i1-1)+j1-1] - c * s[ibl*(ipv-1)+j1-1];
}
f[ig-1] = f[ig-1] - c * f[ipg-1];
}
}
for (j = 2; j <= *nfron; j++) {
/* ----- write variable# and reduced coeff/pivot to disk ----- */
j1 = indx[j];
if (s[ibl*(ipv-1)+j1-1] != 0) {
*icount = *icount + 1;
iba = isbl[j1];
record.variable = iba;
record.coefft = s[ibl*(ipv-1)+j1-1]/pivot;
fwrite(&record, sizeof(record),1,fptr);
s[ibl*(ipv-1)+j1-1] = 0;
}
}
*icount = *icount + 1;
/* ----- write eliminated variable# and rhs/pivot to disk ----- */
record.variable = ipg;
record.coefft = f[ipg-1]/pivot;
fwrite(&record, sizeof(record),1,fptr);
f[ipg-1] = 0;
/* ----- (*ntogo) into (1); (*nfron) into (*ntogo) */
/* ----- ipv into (*nfron) and reduce front & *ntogo sizes by 1 */
if (*ntogo > 1)
indx[1] = indx[*ntogo];
indx[*ntogo] = indx[*nfron];
indx[*nfron] = ipv;
*nfron = *nfron - 1;
*ntogo = *ntogo - 1;
}
return(0);
}
backsub(icount,f)
long int icount;
double *f;
{
struct data record;
long int offset;
int n1,n2;
/* ===== backsubstitution */;
while (icount > 0) {
offset = (icount-1) * sizeof(record);
fseek(fptr,offset,0);
fread(&record,sizeof(record),1,fptr);
icount = icount -1;
n1 = record.variable;
f[n1-1] = record.coefft;
while (icount > 0) {
offset = (icount-1) * sizeof(record);
fseek(fptr,offset,0);
fread(&record,sizeof(record),1,fptr);
icount = icount -1;
n2 = record.variable;
if (n2 == 0)
break;
f[n1-1] = f[n1-1] - record.coefft * f[n2-1];
}
}
return(0);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -