📄 prob3big.cpp
字号:
#include<stdafx.h>
#include<stdio.h>
#include<math.h>
#include "fkn.h"
#include "fknDlg.h"
#include "complex.h"
#include "mesh.h"
#include "spars.h"
#include "FemmeDocCore.h"
#define Log log
BOOL CFemmeDocCore::StaticAxisymmetric(CBigLinProb &L)
{
int i,j,k,s,w,sdi_iter,sdin;
double Me[3][3],Mx[3][3],My[3][3],Mn[3][3];
double l[3],p[3],q[3],g[3],be[3],u[3],v[3],res,lastres,dv,vol;
int n[3]; // numbers of nodes for a particular element;
double a,K,r,t,x,y,B,mu,R,rn[3],a_hat,R_hat,Cduct;
double c=PI*4.e-05;
double units[]={2.54,0.1,1.,100.,0.00254,1.e-04};
double *V_old,*V_sdi,*CircInt1,*CircInt2,*CircInt3;;
int flag,Iter=0,pctr;
BOOL LinearFlag=TRUE;
BOOL SDIflag=FALSE;
res=0;
CElement *El;
V_old=(double *) calloc(NumNodes,sizeof(double));
extRo*=units[LengthUnits];
extRi*=units[LengthUnits];
extZo*=units[LengthUnits];
// check to see if any circuits have been defined and process them;
if (NumCircProps>0)
{
CircInt1=(double *)calloc(NumCircProps,sizeof(double));
CircInt2=(double *)calloc(NumCircProps,sizeof(double));
CircInt3=(double *)calloc(NumCircProps,sizeof(double));
for(i=0;i<NumEls;i++){
if(meshele[i].lbl>=0)
if(labellist[meshele[i].lbl].InCircuit!=-1){
El=&meshele[i];
// get element area;
for(k=0;k<3;k++) n[k]=El->p[k];
p[0]=meshnode[n[1]].y - meshnode[n[2]].y;
p[1]=meshnode[n[2]].y - meshnode[n[0]].y;
p[2]=meshnode[n[0]].y - meshnode[n[1]].y;
q[0]=meshnode[n[2]].x - meshnode[n[1]].x;
q[1]=meshnode[n[0]].x - meshnode[n[2]].x;
q[2]=meshnode[n[1]].x - meshnode[n[0]].x;
a=(p[0]*q[1]-p[1]*q[0])/2.;
r=(meshnode[n[0]].x+meshnode[n[1]].x+meshnode[n[2]].x)/3.;
// if coils are wound, they act like they have
// a zero "bulk" conductivity...
Cduct=blockproplist[El->blk].Cduct;
if (abs(labellist[El->lbl].Turns)!=1) Cduct=0;
// evaluate integrals;
CircInt1[labellist[El->lbl].InCircuit]+=a;
CircInt2[labellist[El->lbl].InCircuit]+=
100.*a*Cduct/r;
CircInt3[labellist[El->lbl].InCircuit]+=
blockproplist[El->blk].Jr*a*100.;
}
}
for(i=0;i<NumCircProps;i++)
{
if (circproplist[i].CircType==0)
{
if(CircInt2[i]==0){
circproplist[i].Case=1;
if(CircInt1[i]==0.) circproplist[i].J=0.;
else circproplist[i].J=0.01*(circproplist[i].Amps_re -
CircInt3[i])/CircInt1[i];
}
else{
circproplist[i].Case=0;
circproplist[i].dV=-0.01*(circproplist[i].Amps_re -
CircInt3[i])/CircInt2[i];
}
}
else{
circproplist[i].Case=0;
circproplist[i].dV=circproplist[i].dVolts_re;
}
}
}
// check to see if there are any SDI boundaries...
// lineproplist[ meshele[i].e[j] ].BdryFormat==0
for(i=0;i<NumLineProps;i++)
if(lineproplist[i].BdryFormat==3) SDIflag=TRUE;
if(SDIflag==TRUE){
// there is an SDI boundary defined; check to see if it is in use
SDIflag=FALSE;
for(i=0;i<NumEls;i++)
for(j=0;j<3;j++)
if (lineproplist[meshele[i].e[j]].BdryFormat==3){
SDIflag=TRUE;
printf("Problem has SDI boundaries\n");
i=NumEls;
j=3;
}
}
if (SDIflag==TRUE)
{
V_sdi=(double *) calloc(NumNodes,sizeof(double));
sdin=2;
}
else sdin=1;
// first, need to define permeability in each block. In nonlinear
// case, this is sort of a hassle. Linear permeability is simply
// copied from the associated block definition, but nonlinear
// permeability must be updated from iteration to iteration...
// build element matrices using the matrices derived in Allaire's book.
do{
for(sdi_iter=0; sdi_iter<sdin; sdi_iter++)
{
TheView->SetDlgItemText(IDC_FRAME1,"Matrix Construction");
TheView->m_prog1.SetPos(0);
pctr=0;
if(Iter>0) L.Wipe();
for(i=0;i<NumEls;i++){
// update ``building matrix'' progress bar...
j=(i*20)/NumEls+1;
if(j>pctr){
j=pctr*5; if (j>100) j=100;
TheView->m_prog1.SetPos(j);
pctr++;
}
// zero out Me, be;
for(j=0;j<3;j++){
for(k=0;k<3;k++){
Me[j][k]=0.;
Mx[j][k]=0;
My[j][k]=0;
Mn[j][k]=0;
}
be[j]=0.;
}
// Determine shape parameters.
// l == element side lengths;
// p corresponds to the `b' parameter in Allaire
// q corresponds to the `c' parameter in Allaire
El=&meshele[i];
for(k=0;k<3;k++){
n[k]=El->p[k];
rn[k]=meshnode[n[k]].x;
}
p[0]=meshnode[n[1]].y - meshnode[n[2]].y;
p[1]=meshnode[n[2]].y - meshnode[n[0]].y;
p[2]=meshnode[n[0]].y - meshnode[n[1]].y;
q[0]=meshnode[n[2]].x - meshnode[n[1]].x;
q[1]=meshnode[n[0]].x - meshnode[n[2]].x;
q[2]=meshnode[n[1]].x - meshnode[n[0]].x;
g[0]=(meshnode[n[2]].x + meshnode[n[1]].x)/2.;
g[1]=(meshnode[n[0]].x + meshnode[n[2]].x)/2.;
g[2]=(meshnode[n[1]].x + meshnode[n[0]].x)/2.;
for(j=0,k=1;j<3;k++,j++){
if (k==3) k=0;
l[j]=sqrt( pow(meshnode[n[k]].x-meshnode[n[j]].x,2.) +
pow(meshnode[n[k]].y-meshnode[n[j]].y,2.) );
}
a=(p[0]*q[1]-p[1]*q[0])/2.;
R=(meshnode[n[0]].x+meshnode[n[1]].x+meshnode[n[2]].x)/3.;
for(j=0,a_hat=0;j<3;j++) a_hat+=(rn[j]*rn[j]*p[j]/(4.*R));
vol=2.*R*a_hat;
for(j=0,flag=0;j<3;j++) if(rn[j]<1.e-06) flag++;
switch(flag)
{
case 2:
R_hat=R;
break;
case 1:
if(rn[0]<1.e-06){
if (fabs(rn[1]-rn[2])<1.e-06) R_hat=rn[2]/2.;
else R_hat=(rn[1] - rn[2])/(2.*log(rn[1]) - 2.*log(rn[2]));
}
if(rn[1]<1.e-06){
if (fabs(rn[2]-rn[0])<1.e-06) R_hat=rn[0]/2.;
else R_hat=(rn[2] - rn[0])/(2.*log(rn[2]) - 2.*log(rn[0]));
}
if(rn[2]<1.e-06){
if (fabs(rn[0]-rn[1])<1.e-06) R_hat=rn[1]/2.;
else R_hat=(rn[0] - rn[1])/(2.*log(rn[0]) - 2.*log(rn[1]));
}
break;
default:
if (fabs(q[0])<1.e-06)
R_hat=(q[1]*q[1])/(2.*(-q[1] + rn[0]*log(rn[0]/rn[2])));
else if (fabs(q[1])<1.e-06)
R_hat=(q[2]*q[2])/(2.*(-q[2] + rn[1]*log(rn[1]/rn[0])));
else if (fabs(q[2])<1.e-06)
R_hat=(q[0]*q[0])/(2.*(-q[0] + rn[2]*log(rn[2]/rn[1])));
else
R_hat=-(q[0]*q[1]*q[2])/
(2.*(q[0]*rn[0]*log(rn[0]) +
q[1]*rn[1]*log(rn[1]) +
q[2]*rn[2]*log(rn[2])));
break;
}
// Mr Contribution
// Derived from flux formulation with c0 + c1 r^2 + c2 z
// interpolation in the element.
K=(-1./(2.*a_hat*R));
for(j=0;j<3;j++)
for(k=j;k<3;k++)
Mx[j][k] += K*p[j]*rn[j]*p[k]*rn[k];
// need this loop to avoid singularities. This just puts something
// on the main diagonal of nodes that are on the r=0 line.
// The program later sets these nodes to zero, but it's good to
// for scaling reasons to grab entries from the neighboring diagonals
// rather than just setting these entries to 1 or something....
for(j=0;j<3;j++)
if (rn[j]<1.e-06) Mx[j][j]+=Mx[0][0]+Mx[1][1]+Mx[2][2];
// Mz Contribution;
// Derived from flux formulation with c0 + c1 r^2 + c2 z
// interpolation in the element.
K=(-1./(2.*a_hat*R_hat));
for(j=0;j<3;j++)
for(k=j;k<3;k++)
My[j][k] += K*(q[j]*rn[j])*(q[k]*rn[k])*
(g[j]/R)*(g[k]/R);
// Fill out rest of entries of Mx and My;
Mx[1][0]=Mx[0][1]; Mx[2][0]=Mx[0][2]; Mx[2][1]=Mx[1][2];
My[1][0]=My[0][1]; My[2][0]=My[0][2]; My[2][1]=My[1][2];
// contributions to Me, be from derivative boundary conditions;
for(j=0;j<3;j++){
if (El->e[j] >= 0)
if (lineproplist[El->e[j]].BdryFormat==2)
{
// conversion factor is 10^(-4) (I think...)
k=j+1; if(k==3) k=0;
r=(meshnode[n[j]].x+meshnode[n[k]].x)/2.;
K=-0.0001*c*2.*r*lineproplist[ El->e[j] ].c0*l[j]/6.;
k=j+1; if(k==3) k=0;
Me[j][j]+=K*2.;
Me[k][k]+=K*2.;
Me[j][k]+=K;
Me[k][j]+=K;
K=(lineproplist[ El->e[j] ].c1*l[j]/2.)*0.0001*2*r;
be[j]+=K;
be[k]+=K;
}
}
// contribution to be from current density in the block
for(j=0;j<3;j++)
{
if(labellist[El->lbl].InCircuit>=0)
{
k=labellist[El->lbl].InCircuit;
if(circproplist[k].Case==1) t=circproplist[k].J.Re();
if(circproplist[k].Case==0)
t=-100.*circproplist[k].dV.Re()*blockproplist[El->blk].Cduct/R;
}
else t=0;
K=-2.*R*(blockproplist[El->blk].Jr+t)*a/3.;
be[j]+=K;
}
// contribution to be from magnetization in the block;
for(j=0;j<3;j++){
t=labellist[El->lbl].MagDir;
k=j+1; if(k==3) k=0;
r=(meshnode[n[j]].x+meshnode[n[k]].x)/2.;
// need to scale so that everything is in proper units...
// conversion is 0.0001
K=-0.0001*r*blockproplist[El->blk].H_c*(
cos(t*PI/180.)*(meshnode[n[k]].x-meshnode[n[j]].x) +
sin(t*PI/180.)*(meshnode[n[k]].y-meshnode[n[j]].y) );
be[j]+=K;
be[k]+=K;
}
// update permeability for the element;
if (Iter==0){
k=meshele[i].blk;
if (blockproplist[k].BHpoints != 0) LinearFlag=FALSE;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -