📄 prob3big.cpp
字号:
if (blockproplist[k].LamType==0){
mu=blockproplist[k].LamFill;
meshele[i].mu1=blockproplist[k].mu_x*mu;
meshele[i].mu2=blockproplist[k].mu_y*mu;
}
if (blockproplist[k].LamType==1){
mu=blockproplist[k].LamFill;
K=blockproplist[k].mu_x;
meshele[i].mu1=K*mu + (1.-mu);
meshele[i].mu2=K/(mu + K*(1.-mu));
}
if (blockproplist[k].LamType==2){
mu=blockproplist[k].LamFill;
K=blockproplist[k].mu_y;
meshele[i].mu1=K*mu + (1.-mu);
meshele[i].mu2=K/(mu + K*(1.-mu));
}
if (blockproplist[k].LamType>2)
{
meshele[i].mu1=1;
meshele[i].mu2=1;
}
}
else{
k=meshele[i].blk;
if ((blockproplist[k].LamType==0) &&
(meshele[i].mu1==meshele[i].mu2)
&&(blockproplist[k].BHpoints>0))
{
// Derive B directly from energy;
v[0]=0;v[1]=0;v[2]=0;
for(j=0;j<3;j++)
for(w=0;w<3;w++)
v[j]+=(Mx[j][w]+My[j][w])*L.V[n[w]];
for(j=0,dv=0;j<3;j++) dv+=L.V[n[j]]*v[j];
dv*=(10000.*c*c/vol);
B=sqrt(fabs(dv));
// find out new mu from saturation curve;
blockproplist[k].GetBHProps(B,mu,dv);
mu=1./(muo*mu);
meshele[i].mu1=mu;
meshele[i].mu2=mu;
for(j=0;j<3;j++){
for(w=0,v[j]=0;w<3;w++)
v[j]+=(Mx[j][w]+My[j][w])*L.V[n[w]];
}
K=-200.*c*c*c*dv/vol;
for(j=0;j<3;j++)
for(w=0;w<3;w++)
Mn[j][w]=K*v[j]*v[w];
}
if ((blockproplist[k].LamType==1) && (blockproplist[k].BHpoints>0)){
// Derive B directly from energy;
t=blockproplist[k].LamFill;
v[0]=0;v[1]=0;v[2]=0;
for(j=0;j<3;j++)
for(w=0;w<3;w++)
v[j]+=(Mx[j][w]+My[j][w]/(t*t))*L.V[n[w]];
for(j=0,dv=0;j<3;j++) dv+=L.V[n[j]]*v[j];
dv*=(10000.*c*c/vol);
B=sqrt(fabs(dv));
// Evaluate BH curve
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/(vol);
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)){
// Derive B directly from energy;
t=blockproplist[k].LamFill;
v[0]=0;v[1]=0;v[2]=0;
for(j=0;j<3;j++)
for(w=0;w<3;w++)
v[j]+=(Mx[j][w]/(t*t)+My[j][w])*L.V[n[w]];
for(j=0,dv=0;j<3;j++) dv+=L.V[n[j]]*v[j];
dv*=(10000.*c*c/vol);
B=sqrt(fabs(dv));
// Evaluate BH curve
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/(vol);
for(j=0;j<3;j++)
for(w=0;w<3;w++)
Mn[j][w]=K*(v[j]*u[w]+v[w]*u[j]);
}
}
// "Warp" the permeability of this element if part of
// the conformally mapped external region
if((labellist[meshele[i].lbl].IsExternal) && (Iter==0) && (sdi_iter==0))
{
double Z=(meshnode[n[0]].y+meshnode[n[1]].y+meshnode[n[2]].y)/3. - extZo;
double kludge=(R*R+Z*Z)*extRi/(extRo*extRo*extRo);
meshele[i].mu1/=kludge;
meshele[i].mu2/=kludge;
}
// 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)
{
r=meshnode[i].x;
L.b[i]+=(0.01*nodeproplist[meshnode[i].bc].Jr*2.*r);
}
// apply fixed boundary conditions at points;
for(i=0;i<NumNodes;i++)
{
if(fabs(meshnode[i].x)<(units[LengthUnits]*1.e-06)) L.SetValue(i,0.);
else
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);
if(x!=0) 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);
if (x!=0) 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.
if (x!=0) 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.
if (x!=0) L.SetValue(meshele[i].p[k],a/c);
}
}
}
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);
// convert answer back to Webers for plotting purposes.
for (i=0;i<NumNodes;i++){
L.b[i]=L.V[i]*c;
L.b[i]*=(meshnode[i].x*0.01*2*PI);
}
free(V_old);
if (SDIflag==TRUE) free(V_sdi);
if(NumCircProps>0){
free(CircInt1);
free(CircInt2);
free(CircInt3);
}
return TRUE;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -