femmviewdoc.cpp

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

						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);
					v2=tn;	
					q=NumList[k];
				}
				else e=&meshelem[nxt];
			}
			b1[i]/=R;
			b2[i]/=R;
		}

			// check to see if angle of corner is too sharp to apply
			// this rule; really only does right if the interface is flat;
			flag=FALSE;	
			// if there is only one edge, approx is ok;
			if ((abs(v1)<0.9) || (abs(v2)<0.9)) flag=TRUE;
			// if the interfaces make less than a 10 degree angle, things are ok;
			if ( (-v1.re*v2.re-v1.im*v2.im) > 0.985) flag=TRUE;
			
			// Otherwise, punt...
			if(flag==FALSE){
				bn=0;
				k=elm.p[i];
				for(j=0;j<NumList[k];j++)
				{
					if(elm.lbl==meshelem[ConList[k][j]].lbl){
						m=ConList[k][j];
						bt.re=sqrt(meshelem[m].B1.re*meshelem[m].B1.re +
							meshelem[m].B2.re*meshelem[m].B2.re);
						bt.im=sqrt(meshelem[m].B1.im*meshelem[m].B1.im +
							meshelem[m].B2.im*meshelem[m].B2.im);
						if(bt.re>bn.re) bn.re=bt.re;
						if(bt.im>bn.im) bn.im=bt.im;
					}
				}
			
				R=sqrt(elm.B1.re*elm.B1.re + elm.B2.re*elm.B2.re);
				if(R!=0)
				{
					b1[i].re=bn.re/R * elm.B1.re;
					b2[i].re=bn.re/R * elm.B2.re;
				}
				else{
					b1[i].re=0;
					b2[i].re=0;
				}

				R=sqrt(elm.B1.im*elm.B1.im + elm.B2.im*elm.B2.im);
				if(R!=0)
				{
					b1[i].im=bn.im/R * elm.B1.im;
					b2[i].im=bn.im/R * elm.B2.im;
				}
				else{
					b1[i].im=0;
					b2[i].im=0;
				}
			}
		}

		// check to see if the point has a point current; if so, just
		// use element average values;
		if (nodeproplist.GetSize()!=0)
		for(j=0;j<nodelist.GetSize();j++)
		{
			if (abs(p-(nodelist[j].x+nodelist[j].y*I))<1.e-08)
			if(nodelist[j].BoundaryMarker>=0)
			{	
				if ((nodeproplist[nodelist[j].BoundaryMarker].Jr!=0) ||
					(nodeproplist[nodelist[j].BoundaryMarker].Ji!=0))
				{
					b1[i]=elm.B1;
					b2[i]=elm.B2;	
				}
			}
		}
		

		//check for special case of node on r=0 axisymmetric; set Br=0;
		if ((fabs(p.re)<1.e-06) && (ProblemType==1)) b1[i].Set(0.,0);
	}
}

void CFemmviewDoc::GetElementB(CElement &elm)
{
	int i,n[3];
	double b[3],c[3],da;

	for(i=0;i<3;i++) n[i]=elm.p[i];

	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]);

	if (ProblemType==0){
		elm.B1=0;
		elm.B2=0;
		for(i=0;i<3;i++)
		{
			elm.B1+=meshnode[n[i]].A*c[i]/(da*LengthConv[LengthUnits]);
			elm.B2-=meshnode[n[i]].A*b[i]/(da*LengthConv[LengthUnits]);
		}
		return;
	}
	else{
		CComplex v[6],dp,dq;
		double R[3],Z[3],r;

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

		// 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]));

		// derivatives w.r.t. p and q:
		dp=(-v[0] + v[2] + 4.*v[3] - 4.*v[5])/3.;
		dq=(-v[0] - 4.*v[1] + 4.*v[3] + v[4])/3.;

		// now, compute flux.
		da*=2.*PI*r*LengthConv[LengthUnits]*LengthConv[LengthUnits];
		elm.B1=-(c[1]*dp+c[2]*dq)/da;
		elm.B2= (b[1]*dp+b[2]*dq)/da;

		return;
	}
}


void CFemmviewDoc::OnReload() 
{
	// TODO: Add your command handler code here
	CString pname = GetPathName();
	if(pname.GetLength()>0)
	{
		OnNewDocument();
		SetPathName(pname,FALSE);
		OnOpenDocument(pname);
	}
}

int CFemmviewDoc::ClosestNode(double x, double y)
{
	int i,j;
	double d0,d1;

	if(nodelist.GetSize()==0) return -1;

	j=0;
	d0=nodelist[0].GetDistance(x,y);
	for(i=0;i<nodelist.GetSize();i++){
		d1=nodelist[i].GetDistance(x,y);
		if(d1<d0){
			d0=d1;
			j=i;
		}
	}

	return j;
}

void CFemmviewDoc::GetLineValues(CXYPlot &p,int PlotType,int NumPlotPoints)
{
	double *q,z,u,dz;
	CComplex pt,n,t;
	int i,j,k,m,elm;
	CPointVals v;
	BOOL flag;
	
	q=(double *)calloc(contour.GetSize(),sizeof(double));
	for(i=1,z=0.;i<contour.GetSize();i++)
	{
		z+=abs(contour[i]-contour[i-1]);
		q[i]=z;
	}
	dz=z/(NumPlotPoints-1);
	
	/*
		m_XYPlotType.AddString("Potential");
		m_XYPlotType.AddString("|B|        (Magnitude of flux density)");
		m_XYPlotType.AddString("B . n      (Normal flux density)");
		m_XYPlotType.AddString("B . t      (Tangential flux density)");
		m_XYPlotType.AddString("|H|        (Magnitude of field intensity)");
		m_XYPlotType.AddString("H . n      (Normal field intensity)");
		m_XYPlotType.AddString("H . t      (Tangential field intensity)");
		m_XYPlotType.AddString("J_eddy     
	*/

	if(Frequency==0){
		switch (PlotType)
		{
			case 0:
				p.Create(NumPlotPoints,2);
				if (ProblemType==0) strcpy(p.lbls[1],"Potential, Wb/m");
				else strcpy(p.lbls[1],"Flux, Wb");
				break;
			case 1:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"|B|, Tesla");
				break;
			case 2:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"B.n, Tesla");
				break;
			case 3:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"B.t, Tesla");
				break;
			case 4:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"|H|, Amp/m");
				break;
			case 5:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"H.n, Amp/m");
				break;
			case 6:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"H.t, Amp/m");
				break;
			default:
				p.Create(NumPlotPoints,2);
				break;
		}
	}
	else{
		switch (PlotType)
		{
			case 0:
				p.Create(NumPlotPoints,4);
				if(ProblemType==0){
					strcpy(p.lbls[1],"|A|, Wb/m");
					strcpy(p.lbls[2],"Re[A], Wb/m");
					strcpy(p.lbls[3],"Im[A], Wb/m");
				}
				else{
					strcpy(p.lbls[1],"|Flux|, Wb");
					strcpy(p.lbls[2],"Re[Flux], Wb");
					strcpy(p.lbls[3],"Im[Flux], Wb");
				}
				break;
			case 1:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"|B|, Tesla");
				break;
			case 2:
				p.Create(NumPlotPoints,4);
				strcpy(p.lbls[1],"|B.n|, Tesla");
				strcpy(p.lbls[2],"Re[B.n], Tesla");
				strcpy(p.lbls[3],"Im[B.n], Tesla");
				break;
			case 3:
				p.Create(NumPlotPoints,4);
				strcpy(p.lbls[1],"|B.t|, Tesla");
				strcpy(p.lbls[2],"Re[B.t], Tesla");
				strcpy(p.lbls[3],"Im[B.t], Tesla");
				break;
			case 4:
				p.Create(NumPlotPoints,2);
				strcpy(p.lbls[1],"|H|, Amp/m");
				break;
			case 5:
				p.Create(NumPlotPoints,4);
				if(ProblemType==0){
					strcpy(p.lbls[1],"|H.n|, Amp/m");
					strcpy(p.lbls[2],"Re[H.n], Amp/m");
					strcpy(p.lbls[3],"Im[H.n], Amp/m");
				}
				break;
			case 6:
				p.Create(NumPlotPoints,4);
				strcpy(p.lbls[1],"|H.t|, Amp/m");
				strcpy(p.lbls[2],"Re[H.t], Amp/m");
				strcpy(p.lbls[3],"Im[H.t], Amp/m");
				break;
			case 7:
				p.Create(NumPlotPoints,4);
				strcpy(p.lbls[1],"|Je|, MA/m^2");
				strcpy(p.lbls[2],"Re[Je], MA/m^2");
				strcpy(p.lbls[3],"Im[Je], MA/m^2");
				break;
			case 8:
				p.Create(NumPlotPoints,4);
				strcpy(p.lbls[1],"|Js+Je|, MA/m^2");
				strcpy(p.lbls[2],"Re[Js+Je], MA/m^2");
				strcpy(p.lbls[3],"Im[Js+Je], MA/m^2");
				break;
			default:
				p.Create(NumPlotPoints,2);
				break;
		}
	}
	
	switch(LengthUnits){
		case 1:
			strcpy(p.lbls[0],"Length, mm");
			break;
		case 2:
			strcpy(p.lbls[0],"Length, cm");
			break;
		case 3:
			strcpy(p.lbls[0],"Length, m");
			break;
		case 4:
			strcpy(p.lbls[0],"Length, mils");
			break;
		case 5:
			strcpy(p.lbls[0],"Length, um");
			break;
		default:
			strcpy(p.lbls[0],"Length, inches");
			break;
	}

	for(i=0,k=1,z=0,elm=-1;i<NumPlotPoints;i++,z+=dz)
	{
		while((z>q[k]) && (k<(contour.GetSize()-1))) k++;
		u=(z-q[k-1])/(q[k]-q[k-1]);
		pt=contour[k-1]+u*(contour[k]-contour[k-1]);
		t=contour[k]-contour[k-1];
		t/=abs(t);
		n = I*t;
		pt+=(n*1.e-06);
		
		if (elm<0) elm=InTriangle(pt.re,pt.im);
		else if (InTriangleTest(pt.re,pt.im,elm)==FALSE)
		{
			flag=FALSE;
			for(j=0;j<3;j++)
				for(m=0;m<NumList[meshelem[elm].p[j]];m++)
				{
					elm=ConList[meshelem[elm].p[j]][m];
					if (InTriangleTest(pt.re,pt.im,elm)==TRUE)
					{
						flag=TRUE;
						m=100;
						j=3;
					}
				}
			if (flag==FALSE) elm=InTriangle(pt.re,pt.im);
		}
		if(elm>=0)
			flag=GetPointValues(pt.re,pt.im,elm,v);
		else flag=FALSE;

		p.M[i][0]=z;
		if ((Frequency==0) && (flag!=FALSE)){
			switch (PlotType){
				case 0:
					p.M[i][1]=v.A.re;
					break;
				case 1:
					p.M[i][1]=sqrt(v.B1.Abs()*v.B1.Abs() + v.B2.Abs()*v.B2.Abs());
					break;
				case 2:
					p.M[i][1]= n.re*v.B1.re + n.im*v.B2.re;
					break;
				case 3:
					p.M[i][1]= t.re*v.B1.re + t.im*v.B2.re;
					break;
				case 4:
					p.M[i][1]=sqrt(v.H1.Abs()*v.H1.Abs() + v.H2.Abs()*v.H2.Abs());
					break;
				case 5:
					p.M[i][1]= n.re*v.H1.re + n.im*v.H2.re;
					break;
				case 6:
					p.M[i][1]= t.re*v.H1.re + t.im*v.H2.re;
					break;
				default:
					p.M[i][1]=0;
					break;
			}
		}
		else if (flag!=FALSE){
			switch (PlotType){
				case 0:
					p.M[i][1]=v.A.Abs();
					p.M[i][2]=v.A.re;
					p.M[i][3]=v.A.im;
					break;
				case 1:
					p.M[i][1]=sqrt(v.B1.Abs()*v.B1.Abs() + v.B2.Abs()*v.B2.Abs());
					break;
				case 2:
					p.M[i][2]= n.re*v.B1.re + n.im*v.B2.re;
					p.M[i][3]= n.re*v.B1.im + n.im*v.B2.im;
					p.M[i][1]=sqrt(p.M[i][2]*p.M[i][2] +
						           p.M[i][3]*p.M[i][3]);
					break;
				case 3:
					p.M[i][2]= t.re*v.B1.re + t.im*v.B2.re;
					p.M[i][3]= t.re*v.B1.im + t.im*v.B2.im;
					p.M[i][1]=sqrt(p.M[i][2]*p.M[i][2] +
						           p.M[i][3]*p.M[i][3]);
					break;
				case 4:
					p.M[i][1]=sqrt(v.H1.Abs()*v.H1.Abs() + v.H2.Abs()*v.H2.Abs());
					break;
				case 5:
					p.M[i][2]= n.re*v.H1.re + n.im*v.H2.re;
					p.M[i][3]= n.re*v.H1.im + n.im*v.H2.im;
					p.M[i][1]=sqrt(p.M[i][2]*p.M[i][2] +
						           p.M[i][3]*p.M[i][3]);
					break;
				case 6:
					p.M[i][2]= t.re*v.H1.re + t.im*v.H2.re;
					p.M[i][3]= t.re*v.H1.im + t.im*v.H2.im;
					p.M[i][1]=sqrt(p.M[i][2]*p.M[i][2] +
						           p.M[i][3]*p.M[i][3]);
					break;
				case 7:
					p.M[i][2]= v.Je.re;
					p.M[i][3]= v.Je.im;
					p.M[i][1]= abs(v.Je);
					break;
				case 8:
					p.M[i][2]= v.Je.re+v.Js.re;
					p.M[i][3]= v.Je.im+v.Js.im;
					p.M[i][1]= abs(v.Je+v.Js);
					break;
				default:
					p.M[i][1]=0;
					break;
			}
		}
	}
		
	free(q);
}

BOOL CFemmviewDoc::InTriangleTest(double x, double y, int i)
{
	int j,k;
	double z;
	BOOL InFlag;

	if(i<0) return FALSE;

	for(j=0,InFlag=TRUE;((j<3) && (InFlag==TRUE));j++)
	{
		k=j+1; if(k==3) k=0;
		z=(meshnode[meshelem[i].p[k]].x-meshnode[meshelem[i].p[j]].x)*
		  (y-meshnode[meshelem[i].p[j]].y) -
		  (meshnode[meshelem[i].p[k]].y-meshnode[meshelem[i].p[j]].y)*
		  (x-meshnode[meshelem[i].p[j]].x);
		if(z<0) InFlag=FALSE;
	}

	return InFlag;
}

CPointVals::CPointVals()
{
	A=0;			// vector potential
	B1=0; B2=0;		// flux density
	mu1=1; mu2=1;	// permeability
	H1=0; H2=0;		// field intensity
	Je=0; Js=0;		// eddy current and source current densities
	c=0;			// conductivity
	E=0;			// energy stored in the magnetic field
	Ph=0;			// power dissipated by hysteresis
	Pe=0;
	ff=1;
}

CComplex CFemmviewDoc::Ctr(int i)
{
	CComplex p,c;
	int j;
	
	for(j=0,c=0;j<3;j++){
		p.Set(meshnode[meshelem[i].p[j]].x/3.,meshnode[meshelem[i].p[j]].y/3.);
		c+=p;
	}

	return c;
}

double CFemmviewDoc::ElmArea(int i)
{
	int j,n[3];
	double b0,b1,c0,c1;

	for(j=0;j<3;j++) n[j]=meshelem[i].p[j];

	b0=meshnode[n[1]].y - meshnode[n[2]].y;
	b1=meshnode[n[2]].y - meshnode[n[0]].y;
	c0=meshnode[n[2]].x - meshnode[n[1]].x;
	c1=meshnode[n[0]].x - meshnode[n[2]].x;
	return (b0*c1-b1*c0)/2.;

}

double CFemmviewDoc::ElmArea(CElement *elm)
{
	int j,n[3];
	double b0,b1,c0,c1;

	for(j=0;j<3;j++) n[j]=elm->p[j];