prob2big.cpp

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

			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))
			{
				for(j=0,B1=0.,B2=0.;j<3;j++){
					B1+=L.V[n[j]]*q[j];
					B2+=L.V[n[j]]*p[j];
				}
				B=c*sqrt(abs(B1*conj(B1))+abs(B2*conj(B2)))/(0.02*a); 
				// correction for lengths in cm of 1/0.02

#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 for successive approx
				K=-200.*c*c*c*dv/a;
				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;
		}

		// 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){
			K=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]+=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].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 (i=0;i<NumNodes+NumCircProps;i++)
			L.V[i]=(V_sdi[i]+L.V[i])/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);

	for (i=0;i<NumNodes;i++) L.b[i]=(L.V[i]*c);	// convert answer back to AMPS
	for (i=0;i<NumCircProps;i++)
		L.b[NumNodes+i]=(I*c*w*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;
}

BOOL CFemmeDocCore::WriteHarmonic2D(CBigComplexLinProb &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	%.17g\n",meshnode[i].x/cf,
			meshnode[i].y/cf,L.b[i].re,L.b[i].im);
	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	%.17g\n",circproplist[i].dV.Re(),
								      circproplist[i].dV.Im());
		if (circproplist[i].Case==1)
			fprintf(fp,"1	%.17g	%.17g\n",circproplist[i].J.Re(),
									  circproplist[i].J.Im());

		if (circproplist[i].Case==2)
			fprintf(fp,"0	%.17g	%.17g\n",L.b[NumNodes+i].Re(),
									  L.b[NumNodes+i].Im());
	}
*/
	// 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	0\n");
		}
		else{
			if (circproplist[i].Case==0)
				fprintf(fp,"0	%.17g	%.17g\n",circproplist[i].dV.Re(),
										  circproplist[i].dV.Im());
			if (circproplist[i].Case==1)
				fprintf(fp,"1	%.17g	%.17g\n",circproplist[i].J.Re(),
										  circproplist[i].J.Im());

			if (circproplist[i].Case==2)
				fprintf(fp,"0	%.17g	%.17g\n",L.b[NumNodes+i].Re(),
										  L.b[NumNodes+i].Im());
		}
	}

	fclose(fp);
	return TRUE;
}