📄 prob1big.cpp
字号:
for(j=0,B1=0.,B2=0.;j<3;j++){
B1+=L.V[n[j]]*q[j];
B2+=L.V[n[j]]*p[j]/t;
}
B=c*sqrt(B1*B1+B2*B2)/(0.02*a);
blockproplist[k].GetBHProps(B,mu,dv);
mu=1./(muo*mu);
meshele[i].mu1=mu*t;
meshele[i].mu2=mu/(t+mu*(1.-t));
for(j=0;j<3;j++){
for(w=0,v[j]=0,u[j]=0;w<3;w++)
{
v[j]+=(My[j][w]/t+Mx[j][w])*L.V[n[w]];
u[j]+=(My[j][w]/t + t*Mx[j][w])*L.V[n[w]];
}
}
K=-100.*c*c*c*dv/(a);
for(j=0;j<3;j++)
for(w=0;w<3;w++)
Mn[j][w]=K*(v[j]*u[w]+v[w]*u[j]);
}
if ((blockproplist[k].LamType==2) && (blockproplist[k].BHpoints>0)){
t=blockproplist[k].LamFill;
for(j=0,B1=0.,B2=0.;j<3;j++){
B1+=(L.V[n[j]]*q[j])/t;
B2+=L.V[n[j]]*p[j];
}
B=c*sqrt(B1*B1+B2*B2)/(0.02*a);
blockproplist[k].GetBHProps(B,mu,dv);
mu=1./(muo*mu);
meshele[i].mu2=mu*t;
meshele[i].mu1=mu/(t+mu*(1.-t));
for(j=0;j<3;j++){
for(w=0,v[j]=0,u[j]=0;w<3;w++)
{
v[j]+=(Mx[j][w]/t + My[j][w])*L.V[n[w]];
u[j]+=(Mx[j][w]/t + t*My[j][w])*L.V[n[w]];
}
}
K=-100.*c*c*c*dv/(a);
for(j=0;j<3;j++)
for(w=0;w<3;w++)
Mn[j][w]=K*(v[j]*u[w]+v[w]*u[j]);
}
}
// combine block matrices into global matrices;
for(j=0;j<3;j++)
for(k=0;k<3;k++){
Me[j][k]+= (Mx[j][k]/Re(El->mu2) + My[j][k]/Re(El->mu1) + Mn[j][k]);
be[j]+=Mn[j][k]*L.V[n[k]];
}
for (j=0;j<3;j++){
for (k=j;k<3;k++)
L.Put(L.Get(n[j],n[k])-Me[j][k],n[j],n[k]);
L.b[n[j]]-=be[j];
}
}
// add in contribution from point currents;
for(i=0;i<NumNodes;i++)
if(meshnode[i].bc>=0)
L.b[i]+=(0.01*nodeproplist[meshnode[i].bc].Jr);
// apply fixed boundary conditions at points;
for(i=0;i<NumNodes;i++)
if(meshnode[i].bc >=0)
if((nodeproplist[meshnode[i].bc].Jr==0) &&
(nodeproplist[meshnode[i].bc].Ji==0) && (sdi_iter==0))
L.SetValue(i,nodeproplist[meshnode[i].bc].Ar / c);
// apply fixed boundary conditions along segments;
for(i=0;i<NumEls;i++)
for(j=0;j<3;j++){
k=j+1; if(k==3) k=0;
if(meshele[i].e[j]>=0)
if(lineproplist[ meshele[i].e[j] ].BdryFormat==0)
{
if(Coords==0){
// first point on the side;
x=meshnode[meshele[i].p[j]].x;
y=meshnode[meshele[i].p[j]].y;
x/=units[LengthUnits]; y/=units[LengthUnits];
s=meshele[i].e[j];
a=lineproplist[s].A0 + x*lineproplist[s].A1 +
y*lineproplist[s].A2;
// just take ``real'' component.
a*=cos(lineproplist[s].phi*DEG);
L.SetValue(meshele[i].p[j],a/c);
// second point on the side;
x=meshnode[meshele[i].p[k]].x;
y=meshnode[meshele[i].p[k]].y;
x/=units[LengthUnits]; y/=units[LengthUnits];
s=meshele[i].e[j];
a=lineproplist[s].A0 + x*lineproplist[s].A1 +
y*lineproplist[s].A2;
// just take``real'' component.
a*=cos(lineproplist[s].phi*DEG);
L.SetValue(meshele[i].p[k],a/c);
}
else{
// first point on the side;
x=meshnode[meshele[i].p[j]].x;
y=meshnode[meshele[i].p[j]].y;
r=sqrt(x*x+y*y);
if((x==0)&&(y==0)) t=0;
else t=atan2(y,x)/DEG;
r/=units[LengthUnits];
s=meshele[i].e[j];
a=lineproplist[s].A0 + r*lineproplist[s].A1 +
t*lineproplist[s].A2;
a*=cos(lineproplist[s].phi*DEG); // just take ``real'' component.
L.SetValue(meshele[i].p[j],a/c);
// second point on the side;
x=meshnode[meshele[i].p[k]].x;
y=meshnode[meshele[i].p[k]].y;
r=sqrt(x*x+y*y);
if((x==0) && (y==0)) t=0;
else t=atan2(y,x)/DEG;
r/=units[LengthUnits];
s=meshele[i].e[j];
a=lineproplist[s].A0 + r*lineproplist[s].A1 +
t*lineproplist[s].A2;
a*=cos(lineproplist[s].phi*DEG); // just take ``real'' component.
L.SetValue(meshele[i].p[k],a/c);
}
}
}
// Apply SDI boundary condition;
if ((SDIflag==TRUE) && (sdi_iter==1)) for(i=0;i<NumEls;i++)
for(j=0;j<3;j++){
k=j+1; if(k==3) k=0;
if(meshele[i].e[j]>=0)
if(lineproplist[ meshele[i].e[j] ].BdryFormat==3)
{
L.SetValue(meshele[i].p[j],0.);
L.SetValue(meshele[i].p[k],0.);
}
}
// Apply any periodicity/antiperiodicity boundary conditions that we have
for(k=0;k<NumPBCs;k++)
{
if (pbclist[k].t==0) L.Periodicity(pbclist[k].x,pbclist[k].y);
if (pbclist[k].t==1) L.AntiPeriodicity(pbclist[k].x,pbclist[k].y);
}
// solve the problem;
if (SDIflag==FALSE) for(j=0;j<NumNodes;j++) V_old[j]=L.V[j];
else{
if(sdi_iter==0)
for(j=0;j<NumNodes;j++) V_sdi[j]=L.V[j];
else
for(j=0;j<NumNodes;j++){
V_old[j]=V_sdi[j];
V_sdi[j]=L.V[j];
}
}
if (L.PCGSolve(Iter+sdi_iter)==FALSE) return FALSE;
if(sdi_iter==1)
for(j=0;j<NumNodes;j++) L.V[j]=(V_sdi[j]+L.V[j])/2.;
} // end of SDI iteration loop;
if (LinearFlag==FALSE){
for(j=0,x=0,y=0;j<NumNodes;j++){
x+=(L.V[j]-V_old[j])*(L.V[j]-V_old[j]);
y+=(L.V[j]*L.V[j]);
}
if (y==0) LinearFlag=TRUE;
else{
lastres=res;
res=sqrt(x/y);
}
// relaxation if we need it
if(Iter>5)
{
if ((res>lastres) && (Relax>0.125)) Relax/=2.;
else Relax+= 0.1 * (1. - Relax);
for(j=0;j<NumNodes;j++) L.V[j]=Relax*L.V[j]+(1.0-Relax)*V_old[j];
}
// report some results
char outstr[256];
sprintf(outstr,"Newton Iteration(%i) Relax=%.4g",Iter,Relax);
TheView->SetDlgItemText(IDC_FRAME2,outstr);
j=(int) (100.*log10(res)/(log10(Precision)+2.));
if (j>100) j=100;
TheView->m_prog2.SetPos(j);
}
// nonlinear iteration has to have a looser tolerance
// than the linear solver--otherwise, things can't ever
// converge. Arbitrarily choose 100*tolerance.
if((res<100.*Precision) && (Iter>0)) LinearFlag=TRUE;
Iter++;
} while(LinearFlag==FALSE);
for(i=0;i<NumNodes;i++) L.b[i]=L.V[i]*c; // convert answer to Amps
free(V_old);
if(SDIflag==TRUE) free(V_sdi);
if(NumCircProps>0){
free(CircInt1);
free(CircInt2);
free(CircInt3);
}
return TRUE;
}
//=========================================================================
//=========================================================================
BOOL CFemmeDocCore::WriteStatic2D(CBigLinProb &L)
{
// write solution to disk;
char c[1024];
FILE *fp,*fz;
int i,k;
double cf;
double unitconv[]={2.54,0.1,1.,100.,0.00254,1.e-04};
// first, echo input .fem file to the .ans file;
sprintf(c,"%s.fem",PathName);
if((fz=fopen(c,"rt"))==NULL){
MsgBox("Couldn't open %s.fem\n",PathName);
return FALSE;
}
sprintf(c,"%s.ans",PathName);
if((fp=fopen(c,"wt"))==NULL){
MsgBox("Couldn't write to %s.ans\n",PathName);
return FALSE;
}
while(fgets(c,1024,fz)!=NULL) fputs(c,fp);
fclose(fz);
// then print out node, line, and element information
fprintf(fp,"[Solution]\n");
cf=unitconv[LengthUnits];
fprintf(fp,"%i\n",NumNodes);
for(i=0;i<NumNodes;i++)
fprintf(fp,"%.17g %.17g %.17g\n",meshnode[i].x/cf,
meshnode[i].y/cf,L.b[i]);
fprintf(fp,"%i\n",NumEls);
for(i=0;i<NumEls;i++)
fprintf(fp,"%i %i %i %i\n",
meshele[i].p[0],meshele[i].p[1],meshele[i].p[2],meshele[i].lbl);
/*
// print out circuit info
fprintf(fp,"%i\n",NumCircPropsOrig);
for(i=0;i<NumCircPropsOrig;i++){
if (circproplist[i].Case==0)
fprintf(fp,"0 %.17g\n",circproplist[i].dV.Re());
if (circproplist[i].Case==1)
fprintf(fp,"1 %.17g\n",circproplist[i].J.Re());
}
*/
// print out circuit info on a blocklabel by blocklabel basis;
fprintf(fp,"%i\n",NumBlockLabels);
for(k=0;k<NumBlockLabels;k++){
i=labellist[k].InCircuit;
if(i<0) // if block not associated with any particular circuit
{
// print out some "dummy" propeties that say that
// there is a fixed additional current density,
// but that that additional current density is zero.
fprintf(fp,"1 0\n");
}
else{
if (circproplist[i].Case==0)
fprintf(fp,"0 %.17g\n",circproplist[i].dV.Re());
if (circproplist[i].Case==1)
fprintf(fp,"1 %.17g\n",circproplist[i].J.Re());
}
}
fclose(fp);
return TRUE;
}
//=========================================================================
//=========================================================================
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -