femmviewdoc.cpp

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

		{
			if(meshelem[i].lbl!=k)
			{
				blocklist[meshelem[i].lbl].IsSelected=TRUE;
				if (!bMultiplyDefinedLabels) 
				{
					CString msg;
					msg ="Some regions in the problem have been defined\n";
					msg+="by more than one block label.  These potentially\n";
					msg+="problematic regions will appear as selected in\n";
					msg+="the initial view.";
					AfxMessageBox(msg);
					bMultiplyDefinedLabels=TRUE;
				}
			}
		}
	}

	// compute fill factor associated with each block label
	for(k=0;k<blocklist.GetSize();k++)
	{
		GetFillFactor(k);
	}

	FirstDraw=TRUE; 

	return TRUE;
}

/*  Original slooow version 
int CFemmviewDoc::InTriangle(double x, double y)
{
	int i;
	double z;
	BOOL InFlag;

	for(i=0;i<meshelem.GetSize();i++)
	{
		z=(meshelem[i].ctr.re-x)*(meshelem[i].ctr.re-x) +
			(meshelem[i].ctr.im-y)*(meshelem[i].ctr.im-y);
		if(z>meshelem[i].rsqr) InFlag=FALSE;
		else InFlag=InTriangleTest(x,y,i); 

		if(InFlag==TRUE) return i;
	}

	return (-1);
}
*/

// Si hang improved version of InTriangle()
int CFemmviewDoc::InTriangle(double x, double y)
{
	CMeshNode searchpoint;
	TriEdge searchtri;
	
	searchpoint.x = x;
	searchpoint.y = y;
	searchtri.ele = 0; // Choose a Elem Arbitrarily.
	
	enum LocateResult lres = Locate(searchpoint, &searchtri, (void *)this);
	if(lres != eOutSide) return searchtri.ele;
	else return (-1);
}

BOOL CFemmviewDoc::GetPointValues(double x, double y, CPointVals &u)
{
	int k;
	k=InTriangle(x,y);
	if (k<0) return FALSE;
	GetPointValues(x,y,k,u);
	return TRUE;
}

BOOL CFemmviewDoc::GetPointValues(double x, double y, int k, CPointVals &u)
{
	int i,j,n[3],lbl;
	double a[3],b[3],c[3],da,ravg;
	
	for(i=0;i<3;i++) n[i]=meshelem[k].p[i];
	a[0]=meshnode[n[1]].x * meshnode[n[2]].y - meshnode[n[2]].x * meshnode[n[1]].y;
	a[1]=meshnode[n[2]].x * meshnode[n[0]].y - meshnode[n[0]].x * meshnode[n[2]].y;
	a[2]=meshnode[n[0]].x * meshnode[n[1]].y - meshnode[n[1]].x * meshnode[n[0]].y;
	b[0]=meshnode[n[1]].y - meshnode[n[2]].y;
	b[1]=meshnode[n[2]].y - meshnode[n[0]].y;
	b[2]=meshnode[n[0]].y - meshnode[n[1]].y;	
	c[0]=meshnode[n[2]].x - meshnode[n[1]].x;
	c[1]=meshnode[n[0]].x - meshnode[n[2]].x;
	c[2]=meshnode[n[1]].x - meshnode[n[0]].x;
	da=(b[0]*c[1]-b[1]*c[0]);
	ravg=LengthConv[LengthUnits]*
		(meshnode[n[0]].x + meshnode[n[1]].x + meshnode[n[2]].x)/3.;
	
	GetPointB(x,y,u.B1,u.B2,meshelem[k]);

	u.Hc=0;
	if(blockproplist[meshelem[k].blk].LamType>2)
		u.ff=blocklist[meshelem[k].lbl].FillFactor;
	else u.ff=-1;

	if (Frequency==0){
		u.A=0;
		if(ProblemType==0){
			for(i=0;i<3;i++)
				u.A.re+=meshnode[n[i]].A.re*(a[i]+b[i]*x+c[i]*y)/(da);
		}
		else{
		/* Old way that I interpolated potential in axi case:
			// interpolation from A from nodal points.
			// note that the potential that's actually stored
			// for axisymmetric problems is 2*Pi*r*A, so divide
			// by nodal r to get 2*Pi*A at the nodes.  Linearly
			// interpolate this, then multiply by r at the point
			// of interest to get back to 2*Pi*r*A.
			for(i=0,rp=0;i<3;i++){
				r=meshnode[n[i]].x;
				rp+=meshnode[n[i]].x*(a[i]+b[i]*x+c[i]*y)/da;
				if (r>1.e-6) u.A.re+=meshnode[n[i]].A.re*
					(a[i]+b[i]*x+c[i]*y)/(r*da);
			}
			u.A.re*=rp;
		*/

			
			// a ``smarter'' interpolation.  One based on A can't
			// represent constant flux density very well.
			// This works, but I should re-write it in a more
			// efficient form--it's doing a lot of the work
			// twice, because the a-b-c stuff that's already
			// been computed is ignored.
			double v[6];
			double R[3];
			double Z[3];
			double p,q;

			for(i=0;i<3;i++){
				R[i]=meshnode[n[i]].x;
				Z[i]=meshnode[n[i]].y;
			}

			// corner nodes
			v[0]=meshnode[n[0]].A.re;
			v[2]=meshnode[n[1]].A.re;
			v[4]=meshnode[n[2]].A.re;

			// construct values for mid-side nodes;
			if ((R[0]<1.e-06) && (R[1]<1.e-06))
			v[1]=(v[0]+v[2])/2.;
			else
			v[1]=(R[1]*(3.*v[0] + v[2]) + R[0]*(v[0] + 3.*v[2]))/
				 (4.*(R[0] + R[1]));

			if ((R[1]<1.e-06) && (R[2]<1.e-06))
			v[3]=(v[2]+v[4])/2.;
			else
			v[3]=(R[2]*(3.*v[2] + v[4]) + R[1]*(v[2] + 3.*v[4]))/
				 (4.*(R[1] + R[2]));

			if ((R[2]<1.e-06) && (R[0]<1.e-06))
			v[5]=(v[4]+v[0])/2.;
			else
			v[5]=(R[0]*(3.*v[4] + v[0]) + R[2]*(v[4] + 3.*v[0]))/
				 (4.*(R[2] + R[0]));

			// compute location in element transformed onto
			// a unit triangle;
			p=(b[1]*x+c[1]*y + a[1])/da;
			q=(b[2]*x+c[2]*y + a[2])/da;

			// now, interpolate to get potential...
			u.A.re = v[0] - p*(3.*v[0] - 4.*v[1] + v[2]) +
				     2.*p*p*(v[0] - 2.*v[1] + v[2]) - 
				     q*(3.*v[0] + v[4] - 4.*v[5]) + 
					 2.*q*q*(v[0] + v[4] - 2.*v[5]) + 
			         4.*p*q*(v[0] - v[1] + v[3] - v[5]);

	/*		// "simple" way to do it...
			// problem is that this mucks up things
			// near the centerline, where things ought
			// to look pretty quadratic.
			for(i=0;i<3;i++)
				u.A.re+=meshnode[n[i]].A.re*(a[i]+b[i]*x+c[i]*y)/(da);
	*/
		}
		u.mu1.im=0; u.mu2.im=0;
		GetMu(u.B1.re,u.B2.re,u.mu1.re,u.mu2.re,k);
		u.H1 = u.B1/(Re(u.mu1)*muo);
		u.H2 = u.B2/(Re(u.mu2)*muo);
		if ((d_ShiftH) && (blockproplist[meshelem[k].blk].H_c!=0))
		{
			u.Hc = blockproplist[meshelem[k].blk].H_c*
				exp(I*PI*blocklist[meshelem[k].lbl].MagDir/180.);
			u.H1=u.H1-Re(u.Hc);
			u.H2=u.H2-Im(u.Hc);
		}
		u.Je=0;
		u.Js=blockproplist[meshelem[k].blk].Jr;
		lbl=meshelem[k].lbl;
		j=blocklist[lbl].InCircuit;
		if(j>=0){
 			if(blocklist[lbl].Case==0){
				if (ProblemType==0)
					u.Js-=Re(blocklist[meshelem[k].lbl].o)*
						  blocklist[lbl].dVolts;
				else{

					int tn;
					double R[3];
					for(tn=0;tn<3;tn++) 
					{
						R[tn]=meshnode[n[tn]].x;
						if (R[tn]<1.e-6) R[tn]=ravg;
						else R[tn]*=LengthConv[LengthUnits];
					}
					for(ravg=0.,tn=0;tn<3;tn++)
						ravg+=(1./R[tn])*(a[tn]+b[tn]*x+c[tn]*y)/(da);
					u.Js-=Re(blocklist[meshelem[k].lbl].o)*
					      blocklist[lbl].dVolts*ravg;
				}
			}
			else u.Js+=blocklist[lbl].J;
		}			
		u.c=Re(blocklist[meshelem[k].lbl].o);
		u.E=blockproplist[meshelem[k].blk].DoEnergy(u.B1.re,u.B2.re);
		// maybe should add in local stored energy for continuum wound
		// coils, but the contribution is so small, and I'm so lazy...

		u.Ph=0;
		u.Pe=0;
		return TRUE;
	}
		
	if(Frequency!=0){
		u.A=0;
		if(ProblemType==0)
		{
			for(i=0;i<3;i++)
				u.A+=meshnode[n[i]].A*(a[i]+b[i]*x+c[i]*y)/(da);
		}
		else{
			CComplex v[6];
			double R[3];
			double Z[3];
			double p,q;

			for(i=0;i<3;i++){
				R[i]=meshnode[n[i]].x;
				Z[i]=meshnode[n[i]].y;
			}

			// corner nodes
			v[0]=meshnode[n[0]].A;
			v[2]=meshnode[n[1]].A;
			v[4]=meshnode[n[2]].A;

			// construct values for mid-side nodes;
			if ((R[0]<1.e-06) && (R[1]<1.e-06))
			v[1]=(v[0]+v[2])/2.;
			else
			v[1]=(R[1]*(3.*v[0] + v[2]) + R[0]*(v[0] + 3.*v[2]))/
				 (4.*(R[0] + R[1]));

			if ((R[1]<1.e-06) && (R[2]<1.e-06))
			v[3]=(v[2]+v[4])/2.;
			else
			v[3]=(R[2]*(3.*v[2] + v[4]) + R[1]*(v[2] + 3.*v[4]))/
				 (4.*(R[1] + R[2]));

			if ((R[2]<1.e-06) && (R[0]<1.e-06))
			v[5]=(v[4]+v[0])/2.;
			else
			v[5]=(R[0]*(3.*v[4] + v[0]) + R[2]*(v[4] + 3.*v[0]))/
				 (4.*(R[2] + R[0]));

			// compute location in element transformed onto
			// a unit triangle;
			p=(b[1]*x+c[1]*y + a[1])/da;
			q=(b[2]*x+c[2]*y + a[2])/da;

			// now, interpolate to get potential...
			u.A = v[0] - p*(3.*v[0] - 4.*v[1] + v[2]) +
				     2.*p*p*(v[0] - 2.*v[1] + v[2]) - 
				     q*(3.*v[0] + v[4] - 4.*v[5]) + 
					 2.*q*q*(v[0] + v[4] - 2.*v[5]) + 
			         4.*p*q*(v[0] - v[1] + v[3] - v[5]);
		}
	
		GetMu(u.B1,u.B2,u.mu1,u.mu2,k);
		u.H1 = u.B1/(u.mu1*muo);
		u.H2 = u.B2/(u.mu2*muo);

		u.Js=blockproplist[meshelem[k].blk].Jr + 
			 I*blockproplist[meshelem[k].blk].Ji;
		lbl=meshelem[k].lbl;
		j=blocklist[lbl].InCircuit;
		if(j>=0){
			if(blocklist[lbl].Case==0){
				if (ProblemType==0)
					u.Js-=blocklist[meshelem[k].lbl].o*blocklist[lbl].dVolts;
				else
				{

					int tn;
					double R[3];
					for(tn=0;tn<3;tn++) 
					{
						R[tn]=meshnode[n[tn]].x;
						if (R[tn]<1.e-6) R[tn]=ravg;
						else R[tn]*=LengthConv[LengthUnits];
					}
					for(ravg=0.,tn=0;tn<3;tn++)
						ravg+=(1./R[tn])*(a[tn]+b[tn]*x+c[tn]*y)/(da);
					u.Js-=blocklist[meshelem[k].lbl].o*
					      blocklist[lbl].dVolts*ravg;
				}
			}
			else u.Js+=blocklist[lbl].J;
		}	

		// report just loss-related part of conductivity.
		if (blockproplist[meshelem[k].blk].Cduct!=0)
			u.c=1./Re(1./(blocklist[meshelem[k].lbl].o));
		else u.c=0;
		
		if (blockproplist[meshelem[k].blk].Lam_d!=0) u.c=0;

		// only add in eddy currents if the region is solid 
		if (abs(blocklist[meshelem[k].lbl].Turns)==1)
			u.Je=-I*Frequency*2.*PI*u.c*u.A;

		if(ProblemType!=0){
			if(x!=0)
				u.Je/=(2.*PI*x*LengthConv[LengthUnits]);
			else u.Je=0;
		}

		CComplex z;
		z=(u.H1*u.B1.Conj()) + (u.H2*u.B2.Conj());
		u.E=0.25*z.re;
		// maybe should add in local stored energy for continuum wound
		// coils, but the contribution is so small, and I'm so lazy...

		u.Ph=Frequency*PI*z.im;
		z=u.Js + u.Je;
		u.Pe=1.e06*(z.Abs()*z.Abs())/(u.c*2.);
		return TRUE;
	}	

	return FALSE;
}

void CFemmviewDoc::GetPointB(double x, double y, CComplex &B1, CComplex &B2,
							 CElement &elm)
{
	// elm is a reference to the element that contains the point of interest.
	int i,n[3];
	double da,a[3],b[3],c[3];

	if(Smooth==FALSE){
		B1=elm.B1;
		B2=elm.B2;
		return;
	}

	for(i=0;i<3;i++) n[i]=elm.p[i];
	a[0]=meshnode[n[1]].x * meshnode[n[2]].y - meshnode[n[2]].x * meshnode[n[1]].y;
	a[1]=meshnode[n[2]].x * meshnode[n[0]].y - meshnode[n[0]].x * meshnode[n[2]].y;
	a[2]=meshnode[n[0]].x * meshnode[n[1]].y - meshnode[n[1]].x * meshnode[n[0]].y;
	b[0]=meshnode[n[1]].y - meshnode[n[2]].y;
	b[1]=meshnode[n[2]].y - meshnode[n[0]].y;
	b[2]=meshnode[n[0]].y - meshnode[n[1]].y;	
	c[0]=meshnode[n[2]].x - meshnode[n[1]].x;
	c[1]=meshnode[n[0]].x - meshnode[n[2]].x;
	c[2]=meshnode[n[1]].x - meshnode[n[0]].x;
	da=(b[0]*c[1]-b[1]*c[0]);

	B1.Set(0,0);
	B2.Set(0,0);
	for(i=0;i<3;i++){
		B1+=(elm.b1[i]*(a[i]+b[i]*x+c[i]*y)/da);
		B2+=(elm.b2[i]*(a[i]+b[i]*x+c[i]*y)/da);
	}
}

void CFemmviewDoc::GetNodalB(CComplex *b1, CComplex *b2,CElement &elm)
{
	// elm is a reference to the element that contains the point of interest.
	CComplex p;
	CComplex tn,bn,bt,btu,btv,u1,u2,v1,v2;
	int i,j,k,l,q,m,pt,nxt;
	double r,R,z;
	CElement *e;
	int flag;

	// find nodal values of flux density via a patch method.
	for(i=0;i<3;i++){

		k=elm.p[i];
		p.Set(meshnode[k].x,meshnode[k].y);
		b1[i].Set(0,0);
		b2[i].Set(0,0);
		for(j=0,m=0;j<NumList[k];j++)
			if(elm.lbl==meshelem[ConList[k][j]].lbl) m++;
			else{
				if(Frequency==0){
					if ((blockproplist[elm.blk].mu_x==
							blockproplist[meshelem[ConList[k][j]].blk].mu_x) &&
						(blockproplist[elm.blk].mu_y==
							blockproplist[meshelem[ConList[k][j]].blk].mu_y) &&
						(blockproplist[elm.blk].H_c==
						 blockproplist[meshelem[ConList[k][j]].blk].H_c) &&
						(blocklist[elm.lbl].MagDir==
						 blocklist[meshelem[ConList[k][j]].lbl].MagDir)) m++;
					else if ((elm.blk==meshelem[ConList[k][j]].blk) &&
						(blocklist[elm.lbl].MagDir==
						 blocklist[meshelem[ConList[k][j]].lbl].MagDir)) m++;
				}
				else if ((blockproplist[elm.blk].mu_fdx==
						  blockproplist[meshelem[ConList[k][j]].blk].mu_fdx) &&
						 (blockproplist[elm.blk].mu_fdy==
						  blockproplist[meshelem[ConList[k][j]].blk].mu_fdy)) m++;
			}
		
		if(m==NumList[k]) // normal smoothing method for points
		{                 // away from any boundaries
			for(j=0,R=0;j<NumList[k];j++)
			{
				m=ConList[k][j];
				z=1./abs(p-Ctr(m));
				R+=z;
				b1[i]+=(z*meshelem[m].B1);
				b2[i]+=(z*meshelem[m].B2);
			}
			b1[i]/=R;
			b2[i]/=R;
		}

		else{
			R=0; v1=0; v2=0;
			
			//scan ccw for an interface...
			e=&elm;
			for(q=0;q<NumList[k];q++)
			{
				//find ccw side of the element;
				for(j=0;j<3;j++) if(e->p[j]==k) pt=j;
				pt--;
				if(pt<0) pt=2;
				pt=e->p[pt];

				//scan to find element adjacent to this side;
				for(j=0,nxt=-1;j<NumList[k];j++){
					if(&meshelem[ConList[k][j]]!=e){
						for(l=0;l<3;l++)
							if(meshelem[ConList[k][j]].p[l]==pt)
							nxt=ConList[k][j];
					}
				}

				if(nxt==-1){
					// a special-case punt
					q=NumList[k];
					b1[i]=(e->B1);
					b2[i]=(e->B2);
					v1=1;v2=1;
				}
				else if(elm.lbl!=meshelem[nxt].lbl){
					// we have found two elements on either side of the interface
					// now, we take contribution from B at the center of the 
					// interface side
					tn.Set(meshnode[pt].x-meshnode[k].x,
						  meshnode[pt].y-meshnode[k].y);
					r=(meshnode[pt].x+meshnode[k].x)*LengthConv[LengthUnits]/2.;
					bn=(meshnode[pt].A-meshnode[k].A)/
					   (abs(tn)*LengthConv[LengthUnits]);
					if(ProblemType==1){
						bn/=(-2.*PI*r);
					}
					z=0.5/abs(tn);
					tn/=abs(tn);

					// for the moment, kludge with bt...
					bt=e->B1*tn.re + e->B2*tn.im;
										
					R+=z;
					b1[i]+=(z*tn.re*bt);
					b2[i]+=(z*tn.im*bt);
					b1[i]+=(z*tn.im*bn);
					b2[i]+=(-z*tn.re*bn);
					v1=tn;
					q=NumList[k];
				}
				else e=&meshelem[nxt];
			}

		//scan cw for an interface...
		if(v2==0) // catches the "special-case punt" where we have
		{         // already set nodal B values....
			e=&elm;
			for(q=0;q<NumList[k];q++)
			{
				//find cw side of the element;
				for(j=0;j<3;j++) if(e->p[j]==k) pt=j;
				pt++;
				if(pt>2) pt=0;
				pt=e->p[pt];

				//scan to find element adjacent to this side;
				for(j=0,nxt=-1;j<NumList[k];j++){
					if(&meshelem[ConList[k][j]]!=e){
						for(l=0;l<3;l++)
							if(meshelem[ConList[k][j]].p[l]==pt)
							nxt=ConList[k][j];
					}
				}
				if (nxt==-1){
					// a special-case punt
					q=NumList[k];
					b1[i]=(e->B1);
					b2[i]=(e->B2);
					v1=1;v2=1;
				}
				else if(elm.lbl!=meshelem[nxt].lbl){
					// we have found two elements on either side of the interface
					// now, we take contribution from B at the center of the 
					// interface side
					tn.Set(meshnode[pt].x-meshnode[k].x,
						  meshnode[pt].y-meshnode[k].y);
					r=(meshnode[pt].x+meshnode[k].x)*LengthConv[LengthUnits]/2.;
					bn=(meshnode[pt].A-meshnode[k].A)/
						(abs(tn)*LengthConv[LengthUnits]);
					if(ProblemType==1){