prob4big.cpp

一个2D电磁场FEM计算的VC++源程序

					K = -0.0001*c*2.*r*lineproplist[ El->e[j] ].c0*l[j]/6.;
					Me[j][j]+=2*K;
					Me[k][k]+=2*K;
					Me[j][k]+=K;
					Me[k][j]+=K;	
	
					K = (lineproplist[ El->e[j] ].c1*l[j]/2.)*2.*r*0.0001;
					be[j]+=K;
					be[k]+=K;
				}
			
				if (lineproplist[El->e[j]].BdryFormat==1)
				{
					ds=sqrt(2./(0.4*PI*w*lineproplist[El->e[j]].Sig*
						lineproplist[El->e[j]].Mu));
					K=deg45/(-ds*lineproplist[El->e[j]].Mu*100.);
					K*=(2.*r*l[j]/6.);
					Me[j][j]+=2*K;
					Me[k][k]+=2*K;
					Me[j][k]+=K;
					Me[k][j]+=K;
				}

			}
		}

		// contribution to be from current density in the block
		for(j=0;j<3;j++){
			Jv=0;
			if(labellist[El->lbl].InCircuit>=0)
			{
				k=labellist[El->lbl].InCircuit;
				if(circproplist[k].Case==1) Jv=circproplist[k].J;
				if(circproplist[k].Case==0)
					Jv=-100.*circproplist[k].dV*
					    blockproplist[El->blk].Cduct/R;
			}

			K=-2.*R*(blockproplist[El->blk].Jr+I*blockproplist[El->blk].Ji+Jv)*a/3.;
			be[j]+=K;

			if(labellist[El->lbl].InCircuit>=0){
				k=labellist[El->lbl].InCircuit;
				if(circproplist[k].Case==2)
					L.b[NumNodes+k]+=K/R;
			}
		}

		// do Case 2 circuit stuff for element
		if(labellist[El->lbl].InCircuit>=0){
			k=labellist[El->lbl].InCircuit;
			if(circproplist[k].Case==2){
				K=-2.*I*a*w*blockproplist[meshele[i].blk].Cduct*c;
				for(j=0;j<3;j++)
					L.Put(L.Get(n[j],NumNodes+k)+K/3.,n[j],NumNodes+k);
				L.Put(L.Get(NumNodes+k,NumNodes+k)+K/R,NumNodes+k,NumNodes+k);
			}
		}

/////////////////////////
//
//  Nonlinear Stuff
//
/////////////////////////

		// update permeability for the element;
		if (Iter==0){ 
			k=meshele[i].blk;
			if (blockproplist[k].BHpoints != 0) LinearFlag=FALSE;
			meshele[i].mu1=Mu[k][0];
			meshele[i].mu2=Mu[k][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(ww=0;ww<3;ww++)
						v[j]+=(Mx[j][ww]+My[j][ww])*L.V[n[ww]];
				for(j=0,dv=0;j<3;j++) dv+=conj(L.V[n[j]])*v[j];
				dv*=(10000.*c*c/vol); 
				B=sqrt(abs(dv));

#ifdef NEWTON
				// 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(ww=0,v[j]=0;ww<3;ww++)
						v[j]+=(Mx[j][ww]+My[j][ww])*L.V[n[ww]];
				}

				// Newton-like iteration
				// comment out this block for successive approximation
				K=-200.*c*c*c*dv/vol;
				for(j=0;j<3;j++)
					for(ww=0;ww<3;ww++)
						Mn[j][ww]=K*Re(v[j]*conj(v[ww])); 
#else
			   // find out new mu from saturation curve;
               murel=1./(muo*blockproplist[k].Get_v(B));
               muinc=1./(muo*blockproplist[k].GetdHdB(B));

               // successive approximation;
       //      K=muinc;                            // total incremental
       //      K=murel;                            // total updated
               K=2.*murel*muinc/(murel+muinc);     // averaged
               meshele[i].mu1=K;
               meshele[i].mu2=K;
               K=-(1./murel - 1/K);
               for(j=0;j<3;j++)
                   for(ww=0;ww<3;ww++)
                       Mn[j][ww]=K*(Mx[j][ww]+My[j][ww]); 

#endif
			}	
		}

		// Apply correction for elements subject to prox effects
		if((blockproplist[meshele[i].blk].LamType>2) && (Iter==0) && (sdi_iter==0))
		{
			meshele[i].mu1=labellist[meshele[i].lbl].ProximityMu;
			meshele[i].mu2=labellist[meshele[i].lbl].ProximityMu;
		}

		// "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++){
#ifdef NEWTON
				Me[j][k]+= (Mx[j][k]/(El->mu2) + My[j][k]/(El->mu1) + Mn[j][k]);
#else
				Me[j][k]+= (Mx[j][k]/(El->mu2) + My[j][k]/(El->mu1));
#endif			
				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;
			K = (2.*r*0.01)*(nodeproplist[meshnode[i].bc].Jr +
				  I*nodeproplist[meshnode[i].bc].Ji);
			L.b[i]-=K;
		}

	// add in total current constraints for circuits;
	for(i=0;i<NumCircProps;i++)
		if (circproplist[i].Case==2){
			L.b[NumNodes+i]+=2.*0.01*(circproplist[i].Amps_re +
				                   I*circproplist[i].Amps_im);
		}
	
	// apply fixed boundary conditions at points;
	for(i=0;i<NumNodes;i++)
		if(meshnode[i].x<(units[LengthUnits]*1.e-06)){
			K=0;
			L.SetValue(i,K);
		}
		else if(meshnode[i].bc >=0)
			if((nodeproplist[meshnode[i].bc].Jr==0) &&
			   (nodeproplist[meshnode[i].bc].Ji==0) && (sdi_iter==0))
			{
				K =  (nodeproplist[meshnode[i].bc].Ar 
			      + I*nodeproplist[meshnode[i].bc].Ai) / c;
				L.SetValue(i,K);
			}

	// 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;
					K=(a/c)*exp(I*lineproplist[s].phi*DEG);
					L.SetValue(meshele[i].p[j],K);
					
					// 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;
					K=(a/c)*exp(I*lineproplist[s].phi*DEG);
					L.SetValue(meshele[i].p[k],K);
				}
				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;
					K=(a/c)*exp(I*lineproplist[s].phi*DEG);
					L.SetValue(meshele[i].p[j],K);
					
					// 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;
					K=(a/c)*exp(I*lineproplist[s].phi*DEG);
					L.SetValue(meshele[i].p[k],K);
				}

			}
		}
	
	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.*I);
					L.SetValue(meshele[i].p[k],0.*I);
			}
		}

	// "fix" diagonal entries associated with circuits that have 
	// applied current or voltage that is known a priori
	// so that solver doesn't throw a "singular" flag
	for(j=0;j<NumCircProps;j++)
		if (circproplist[j].Case<2)	L.Put(L.Get(0,0),NumNodes+j,NumNodes+j);
	
	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+NumCircProps;j++) V_old[j]=L.V[j];
	else{
		if(sdi_iter==0)
			for(j=0;j<NumNodes+NumCircProps;j++) V_sdi[j]=L.V[j];
		else
			for(j=0;j<NumNodes+NumCircProps;j++){
				V_old[j]=V_sdi[j];
				V_sdi[j]=L.V[j];
			}
	}

	if (L.PBCGSolveMod(Iter+sdi_iter)==FALSE) return FALSE;
	
	if(sdi_iter==1) 
		for(j=0;j<NumNodes+NumCircProps;j++) L.V[j]=(V_sdi[j]+L.V[j])/2.;

  } //end of SDI loop;

  	if (LinearFlag==FALSE){
		
		for(j=0,x=0,y=0;j<NumNodes;j++){
			x+=Re((L.V[j]-V_old[j])*conj(L.V[j]-V_old[j]));
			y+=Re(L.V[j]*conj(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.1)) Relax/=2.;
            else Relax+= 0.1 * (1. - Relax);
       
            for(j=0;j<NumNodes+NumCircProps;j++) L.V[j]=Relax*L.V[j]+(1.0-Relax)*V_old[j];
        }

       
        // report some results
		char outstr[256];
#ifdef NEWTON
		sprintf(outstr,"Newton Iteration(%i) Relax=%.4g",Iter,Relax);
#else
		sprintf(outstr,"Successive Approx(%i) Relax=%.4g",Iter,Relax);
#endif
        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 (i=0;i<NumNodes;i++) L.b[i]=L.V[i]*c*2.*PI*meshnode[i].x*0.01;
	for (i=0;i<NumCircProps;i++)
		L.b[NumNodes+i]=(I*w*c*0.01*L.V[NumNodes+i]);
	

	// free up space allocated in this routine
	for(k=0;k<NumBlockProps;k++) free(Mu[k]);
	free(Mu);
	free(V_old);
	if (SDIflag==TRUE) free(V_sdi);
	if(NumCircProps>0){
		free(CircInt1);
		free(CircInt2);
		free(CircInt3);
	}
	return TRUE;
}