femmedoc.cpp

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

		return FALSE;
	}

	fprintf(fp,"[Format]      =  4.0\n");
	fprintf(fp,"[Frequency]   =  %.17g\n",Frequency);
	fprintf(fp,"[Precision]   =  %.17g\n",Precision);
	fprintf(fp,"[MinAngle]    =  %.17g\n",MinAngle);
	fprintf(fp,"[Depth]       =  %.17g\n",Depth);
	fprintf(fp,"[LengthUnits] =  ",LengthUnits);
	switch(LengthUnits)
	{
		case 1:
			fprintf(fp,"millimeters\n");
			break;
		case 2:
			fprintf(fp,"centimeters\n");
			break;
		case 3:
			fprintf(fp,"meters\n");
			break;
		case 4:
			fprintf(fp,"mils\n");
			break;
		case 5:
			fprintf(fp,"microns\n");
			break;
		default:
			fprintf(fp,"inches\n");
			break;
	}
	
	if (ProblemType==0) fprintf(fp,"[ProblemType] =  planar\n");
	else{
		fprintf(fp,"[ProblemType] =  axisymmetric\n");
		if ((extRo!=0) && (extRi!=0))
		{
			fprintf(fp,"[extZo] = %.17g\n",extZo);
			fprintf(fp,"[extRo] = %.17g\n",extRo);
			fprintf(fp,"[extRi] = %.17g\n",extRi);
		}
	}

	if (Coords==0) fprintf(fp,"[Coordinates] =  cartesian\n");
	else fprintf(fp,"[Coordinates] =  polar\n");

	s=ProblemNote;
	for(i=0;i<ProblemNote.GetLength();i++){
		if (s[i]==13) s.SetAt(i,'\\');
		if (s[i]==10) s.SetAt(i,'n');
	}
	fprintf(fp,"[Comment]     =  \"%s\"\n",s);

	// write out materials properties stuff...
	fprintf(fp,"[PointProps]   = %i\n",nodeproplist.GetSize());
	for(i=0;i<nodeproplist.GetSize();i++)
	{
		fprintf(fp,"  <BeginPoint>\n");
		fprintf(fp,"    <PointName> = \"%s\"\n",nodeproplist[i].PointName);
		fprintf(fp,"    <I_re> = %.17g\n",nodeproplist[i].Jr);
		fprintf(fp,"    <I_im> = %.17g\n",nodeproplist[i].Ji);
		fprintf(fp,"    <A_re> = %.17g\n",nodeproplist[i].Ar);
		fprintf(fp,"    <A_im> = %.17g\n",nodeproplist[i].Ai);
		fprintf(fp,"  <EndPoint>\n");
	}

	fprintf(fp,"[BdryProps]   = %i\n",lineproplist.GetSize());
	for(i=0;i<lineproplist.GetSize();i++)
	{
		fprintf(fp,"  <BeginBdry>\n");	
		fprintf(fp,"    <BdryName> = \"%s\"\n",lineproplist[i].BdryName);
		fprintf(fp,"    <BdryType> = %i\n",lineproplist[i].BdryFormat);
		fprintf(fp,"    <A_0> = %.17g\n",lineproplist[i].A0);
		fprintf(fp,"    <A_1> = %.17g\n",lineproplist[i].A1);
		fprintf(fp,"    <A_2> = %.17g\n",lineproplist[i].A2);
		fprintf(fp,"    <Phi> = %.17g\n",lineproplist[i].phi);
		fprintf(fp,"    <c0> = %.17g\n",lineproplist[i].c0);
		fprintf(fp,"    <c1> = %.17g\n",lineproplist[i].c1);
		fprintf(fp,"    <Mu_ssd> = %.17g\n",lineproplist[i].Mu);
		fprintf(fp,"    <Sigma_ssd> = %.17g\n",lineproplist[i].Sig);
		fprintf(fp,"  <EndBdry>\n");
	}

	fprintf(fp,"[BlockProps]  = %i\n",blockproplist.GetSize());
	for(i=0;i<blockproplist.GetSize();i++)
	{
		fprintf(fp,"  <BeginBlock>\n");
		fprintf(fp,"    <BlockName> = \"%s\"\n",blockproplist[i].BlockName);
		fprintf(fp,"    <Mu_x> = %.17g\n",blockproplist[i].mu_x);
		fprintf(fp,"    <Mu_y> = %.17g\n",blockproplist[i].mu_y);
		fprintf(fp,"    <H_c> = %.17g\n",blockproplist[i].H_c);
		fprintf(fp,"    <H_cAngle> = %.17g\n",blockproplist[i].Theta_m);
		fprintf(fp,"    <J_re> = %.17g\n",blockproplist[i].Jr);
		fprintf(fp,"    <J_im> = %.17g\n",blockproplist[i].Ji);
		fprintf(fp,"    <Sigma> = %.17g\n",blockproplist[i].Cduct);
		fprintf(fp,"    <d_lam> = %.17g\n",blockproplist[i].Lam_d);
		fprintf(fp,"    <Phi_h> = %.17g\n",blockproplist[i].Theta_hn);
		fprintf(fp,"    <Phi_hx> = %.17g\n",blockproplist[i].Theta_hx);
		fprintf(fp,"    <Phi_hy> = %.17g\n",blockproplist[i].Theta_hy);
		fprintf(fp,"    <LamType> = %i\n",blockproplist[i].LamType);
		fprintf(fp,"    <LamFill> = %.17g\n",blockproplist[i].LamFill);
		fprintf(fp,"    <NStrands> = %i\n",blockproplist[i].NStrands);
		fprintf(fp,"    <WireD> = %.17g\n",blockproplist[i].WireD);
		fprintf(fp,"    <BHPoints> = %i\n",blockproplist[i].BHpoints);
		for(j=0;j<blockproplist[i].BHpoints;j++)
			fprintf(fp,"      %.17g	%.17g\n",blockproplist[i].BHdata[j].re,
								  blockproplist[i].BHdata[j].im);
		fprintf(fp,"  <EndBlock>\n");
	}
 
	fprintf(fp,"[CircuitProps]  = %i\n",circproplist.GetSize());
	for(i=0;i<circproplist.GetSize();i++)
	{
		fprintf(fp,"  <BeginCircuit>\n");
		fprintf(fp,"    <CircuitName> = \"%s\"\n",circproplist[i].CircName);
		fprintf(fp,"    <TotalAmps_re> = %.17g\n",circproplist[i].Amps.Re());
		fprintf(fp,"    <TotalAmps_im> = %.17g\n",circproplist[i].Amps.Im());
		fprintf(fp,"    <CircuitType> = %i\n",circproplist[i].CircType);
		fprintf(fp,"  <EndCircuit>\n");
	}

	// write out node list
	fprintf(fp,"[NumPoints] = %i\n",nodelist.GetSize());
	for(i=0;i<nodelist.GetSize();i++)
	{
		for(j=0,t=0;j<nodeproplist.GetSize();j++)
				if(nodeproplist[j].PointName==nodelist[i].BoundaryMarker) t=j+1;
		fprintf(fp,"%.17g	%.17g	%i	%i\n",nodelist[i].x,nodelist[i].y,t,
			nodelist[i].InGroup);
	}

	// write out segment list
	fprintf(fp,"[NumSegments] = %i\n",linelist.GetSize());
	for(i=0;i<linelist.GetSize();i++)
	{
		for(j=0,t=0;j<lineproplist.GetSize();j++)
				if(lineproplist[j].BdryName==linelist[i].BoundaryMarker) t=j+1;
		fprintf(fp,"%i	%i	",linelist[i].n0,linelist[i].n1);
		if(linelist[i].MaxSideLength<0) fprintf(fp,"-1	");
		else fprintf(fp,"%.17g	",linelist[i].MaxSideLength);
		fprintf(fp,"%i	%i	%i\n",t,linelist[i].Hidden,linelist[i].InGroup);
	}

	// write out arc segment list
	fprintf(fp,"[NumArcSegments] = %i\n",arclist.GetSize());
	for(i=0;i<arclist.GetSize();i++)
	{
		for(j=0,t=0;j<lineproplist.GetSize();j++)
				if(lineproplist[j].BdryName==arclist[i].BoundaryMarker) t=j+1;
		fprintf(fp,"%i	%i	%.17g	%.17g	%i	%i	%i\n",arclist[i].n0,arclist[i].n1,
			arclist[i].ArcLength,arclist[i].MaxSideLength,t,
			arclist[i].Hidden,arclist[i].InGroup);
	}

	// write out list of holes;
	for(i=0,j=0;i<blocklist.GetSize();i++)
		if(blocklist[i].BlockType=="<No Mesh>") j++;
	fprintf(fp,"[NumHoles] = %i\n",j);
	for(i=0,k=0;i<blocklist.GetSize();i++)
		if(blocklist[i].BlockType=="<No Mesh>")
		{
			fprintf(fp,"%.17g	%.17g	%i\n",blocklist[i].x,blocklist[i].y,
				blocklist[i].InGroup);
			k++;
		}
	
	// write out regional attributes
	fprintf(fp,"[NumBlockLabels] = %i\n",blocklist.GetSize()-j);
	for(i=0,k=0;i<blocklist.GetSize();i++)
		if(blocklist[i].BlockType!="<No Mesh>")
		{
			fprintf(fp,"%.17g	%.17g	",blocklist[i].x,blocklist[i].y);
			for(j=0,t=0;j<blockproplist.GetSize();j++)
				if(blockproplist[j].BlockName==blocklist[i].BlockType) t=j+1;
			fprintf(fp,"%i	",t);
			if (blocklist[i].MaxArea>0)
				fprintf(fp,"%.17g	",sqrt(4.*blocklist[i].MaxArea/PI));
			else fprintf(fp,"-1	");
			for(j=0,t=0;j<circproplist.GetSize();j++)
				if(circproplist[j].CircName==blocklist[i].InCircuit) t=j+1;
			fprintf(fp,"%i	%.17g	%i	%i	%i\n",t,blocklist[i].MagDir,
				blocklist[i].InGroup,blocklist[i].Turns,blocklist[i].IsExternal);
			k++;
		}
	fclose(fp);

	return TRUE;
}

BOOL CFemmeDoc::LoadMesh()
{
	int i,j,k,q,nl;
	CString pathname,rootname,infile;
	FILE *fp;
	char s[1024];

	// clear out the old mesh...
	meshnode.RemoveAll();
	meshline.RemoveAll();
	greymeshline.RemoveAll();

	pathname=GetPathName();
	if (pathname.GetLength()==0){
		AfxMessageBox("No mesh to display");
		return FALSE;
	}

	rootname=pathname.Left(pathname.ReverseFind('.'));

	//read meshnodes;
	infile=rootname+".node";
	if((fp=fopen(infile,"rt"))==NULL){
		AfxMessageBox("No mesh to display");
		return FALSE;
	}
	fgets(s,1024,fp);
	sscanf(s,"%i",&k);
	meshnode.SetSize(k);
	CNode node;
	for(i=0;i<k;i++){
		fgets(s,1024,fp);
		sscanf(s,"%i	%lf	%lf",&j,&node.x,&node.y);
		meshnode.SetAt(i,node);
	}
	fclose(fp);

	//read meshlines;
	infile=rootname+".edge";
	if((fp=fopen(infile,"rt"))==NULL){
		AfxMessageBox("No mesh to display");
		return FALSE;
	}
	fgets(s,1024,fp);
	sscanf(s,"%i",&k);
	meshline.SetSize(k);
	fclose(fp);

	infile=rootname+".ele";
	if((fp=fopen(infile,"rt"))==NULL){
		AfxMessageBox("No mesh to display");
		return FALSE;
	}
	fgets(s,1024,fp);
	sscanf(s,"%i",&k);

	CPoint segm;
	int n[3],p;
	for(i=0,nl=0;i<k;i++)
	{
		fgets(s,1024,fp);
		sscanf(s,"%i	%i	%i	%i	%i",&q,&n[0],&n[1],&n[2],&j);
		for(q=0;q<3;q++)
		{
			p=q+1; if(p==3) p=0;
			if (n[p]>n[q])
			{ 
				segm.x=n[p]; 
				segm.y=n[q]; 
				if (j!=0) meshline.SetAt(nl++,segm);
				else greymeshline.Add(segm);
			}
		}
	}
	meshline.SetSize(nl);
	fclose(fp);

	// clear out temporary files
	infile=rootname+".ele";		DeleteFile(infile);
	infile=rootname+".node";	DeleteFile(infile);	
	infile=rootname+".edge";	DeleteFile(infile);
	infile=rootname+".pbc";		DeleteFile(infile);
	infile=rootname+".poly";	DeleteFile(infile);

	return TRUE;
}

void CFemmeDoc::UpdateUndo()
{
	
	int i;
	
	undonodelist.RemoveAll();
	undolinelist.RemoveAll();
	undoarclist.RemoveAll();
	undoblocklist.RemoveAll();

	for(i=0;i<nodelist.GetSize();i++) undonodelist.Add(nodelist[i]);
	for(i=0;i<linelist.GetSize();i++) undolinelist.Add(linelist[i]);
	for(i=0;i<arclist.GetSize();i++) undoarclist.Add(arclist[i]);
	for(i=0;i<blocklist.GetSize();i++) undoblocklist.Add(blocklist[i]);


}

void CFemmeDoc::Undo()
{
	int i;

	CArray< CNode, CNode&>             tempnodelist;
	CArray< CSegment, CSegment&>       templinelist;
	CArray< CArcSegment, CArcSegment&> temparclist;
	CArray< CBlockLabel, CBlockLabel&> tempblocklist;

	tempnodelist.RemoveAll();
	templinelist.RemoveAll();
	temparclist.RemoveAll();
	tempblocklist.RemoveAll();

	for(i=0;i<nodelist.GetSize();i++) tempnodelist.Add(nodelist[i]);
	for(i=0;i<linelist.GetSize();i++) templinelist.Add(linelist[i]);
	for(i=0;i<arclist.GetSize();i++) temparclist.Add(arclist[i]);
	for(i=0;i<blocklist.GetSize();i++) tempblocklist.Add(blocklist[i]);

	nodelist.RemoveAll();
	linelist.RemoveAll();
	arclist.RemoveAll();
	blocklist.RemoveAll();

	for(i=0;i<undonodelist.GetSize();i++) nodelist.Add(undonodelist[i]);
	for(i=0;i<undolinelist.GetSize();i++) linelist.Add(undolinelist[i]);
	for(i=0;i<undoarclist.GetSize();i++) arclist.Add(undoarclist[i]);
	for(i=0;i<undoblocklist.GetSize();i++) blocklist.Add(undoblocklist[i]);

	undonodelist.RemoveAll();
	undolinelist.RemoveAll();
	undoarclist.RemoveAll();
	undoblocklist.RemoveAll();

	for(i=0;i<tempnodelist.GetSize();i++) undonodelist.Add(tempnodelist[i]);
	for(i=0;i<templinelist.GetSize();i++) undolinelist.Add(templinelist[i]);
	for(i=0;i<temparclist.GetSize();i++) undoarclist.Add(temparclist[i]);
	for(i=0;i<tempblocklist.GetSize();i++) undoblocklist.Add(tempblocklist[i]);
}

BOOL CFemmeDoc::ScanPreferences()
{
	FILE *fp;
	CString fname;

	fname=BinDir+"femme.cfg";

	fp=fopen(fname,"rt");
	if (fp!=NULL)
	{
		BOOL flag=FALSE;
		char s[1024];
		char q[1024];
		char *v;
	
		// parse the file
		while (fgets(s,1024,fp)!=NULL)
		{
			sscanf(s,"%s",q);

			if( _strnicmp(q,"<Precision>",11)==0)
			{
			  v=StripKey(s);
			  sscanf(v,"%lf",&d_prec);
			  q[0]=NULL;
			}
	
			if( _strnicmp(q,"<MinAngle>",10)==0)
			{
			  v=StripKey(s);
			  sscanf(v,"%lf",&d_minangle);
			  q[0]=NULL;
			}

			if( _strnicmp(q,"<Frequency>",11)==0)
			{
			  v=StripKey(s);
			  sscanf(v,"%lf",&d_freq);
			  q[0]=NULL;
			}

			if( _strnicmp(q,"<Depth>",7)==0)
			{
			  v=StripKey(s);
			  sscanf(v,"%lf",&d_depth);
			  q[0]=NULL;
			}

			if( _strnicmp(q,"<Coordinates>",13)==0)
			{
			  v=StripKey(s);
			  sscanf(v,"%i",&d_coord);
			  q[0]=NULL;
			}
				
			if( _strnicmp(q,"<LengthUnits>",13)==0)
			{
			  v=StripKey(s);
			  sscanf(v,"%i",&d_length);
			  q[0]=NULL;
			}

			if( _strnicmp(q,"<ProblemType>",13)==0)
			{
			  v=StripKey(s);
			  sscanf(v,"%i",&d_type);
			  q[0]=NULL;
			}

		}
		fclose(fp);
			return TRUE;
	}
	
	return FALSE;
}


//--------------------------------------------------------------


void CFemmeView::lnuMakeMesh()
{
	OnMakeMesh();
}

void CFemmeView::lnu_purge_mesh()
{
	OnPurgemesh();
}

void CFemmeView::lnu_show_mesh()
{
	OnShowMesh();
}

void CFemmeView::lnu_analyze(int bShow)
{
	if (bShow) bLinehook=HiddenLua;
	else bLinehook=NormalLua;
	OnMenuAnalyze();
}

void CFemmeView::lua_zoomnatural()
{
	OnZoomNatural();
}

void CFemmeView::lua_zoomout()
{
	OnZoomOut();
}

void CFemmeView::lua_zoomin()
{
	OnZoomIn();
}

BOOL CFemmeDoc::OnCmdMsg(UINT nID, int nCode, void* pExtra, AFX_CMDHANDLERINFO* pHandlerInfo) 
{
	// TODO: Add your specialized code here and/or call the base class
	if (bLinehook!=FALSE) return TRUE;
	return CDocument::OnCmdMsg(nID, nCode, pExtra, pHandlerInfo);
}

void CFemmeDoc::OnEditExterior() 
{
	// TODO: Add your command handler code here
	CExteriorProps dlg;

	dlg.m_Ro=extRo;
	dlg.m_Ri=extRi;
	dlg.m_Zo=extZo;
	if(dlg.DoModal()==IDOK)
	{
		extRo=dlg.m_Ro;
		extRi=dlg.m_Ri;
		extZo=dlg.m_Zo;
	}
}

BOOL CFemmeDoc::CanCreateRadius(int n)
{
	// check to see if a selected point, specified by n, can be
	// converted to a radius.  To be able to be converted to a radius,
	// the point must be an element of either 2 lines, 2 arcs, or
	// 1 line and 1 arc.

	int j,k;

	for(k=0,j=0;k<linelist.GetSize();k++)
		if ((linelist[k].n0==n) || (linelist[k].n1==n)) j++;
	for(k=0;k<arclist.GetSize();k++)
		if ((arclist[k].n0==n) || (arclist[k].n1==n)) j++;

	if (j==2) return TRUE;

	return FALSE;
}

BOOL CFemmeDoc::CreateRadius(int n, double r)
{
	// replace the node indexed by n with a radius of r

	if(r<=0) return FALSE;

	int nseg=0;
	int narc=0;
	int k,arc[2],seg[2];

	for(k=0;k<linelist.GetSize();k++)
		if ((linelist[k].n0==n) || (linelist[k].n1==n)) seg[nseg++]=k;
	for(k=0;k<arclist.GetSize();k++)
		if ((arclist[k].n0==n) || (arclist[k].n1==n)) arc[narc++]=k;
	if ((narc+nseg)!=2) return FALSE;

	// there are three valid cases:
	switch (nseg-narc)
	{
		case 0:  // One arc and one line
		{
			CComplex c,u,p0,p1,q,p[4],v[8],i1[8],i2[8];
			double rc,b,R[4],phi;
			CArcSegment ar;
			int j,m;

			// inherit the boundary condition from the arc so that 
			// we can apply it to the newly created arc later;
			ar.InGroup       =arclist[arc[0]].InGroup;
			ar.BoundaryMarker=arclist[arc[0]].BoundaryMarker;
			
			// get the center and radius of the circle associated with the arc;
			GetCircle(arclist[arc[0]],c,rc);

			// get the locations of the endpoints of the segment;
			p0=nodelist[n].CC();
			if(linelist[seg[0]].n0==n) 
				p1=nodelist[linelist[seg[0]].n1].CC();
			else
				p1=nodelist[linelist[seg[0]].n0].CC();

			u=(p1-p0)/abs(p1-p0);  // unit vector along the line
			q=p0 + u*Re((c-p0)/u); // closest point on line to center of circle
			u=(q-c)/abs(q-c); // unit vector from center to closest point on line;

			p[0]=q+r*u; R[0]=rc+r;
			p[1]=q-r*u; R[1]=rc+r;
			p[2]=q+r*u; R[2]=rc-r;
			p[3]=q-r*u; R[3]=rc-r;
		
			for(j=0,k=0;k<4;k++)
			{
				b=R[k]*R[k]-abs(p[k]-c)*abs(p[k]-c);