writepoly.cpp

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

// implementation of various incarnations of calls
// to triangle from the CFemmeDoc class

#include "stdafx.h"
#include "femme.h"
#include "femmeDoc.h"
#include "probdlg.h"
#include "PtProp.h"
#include "OpBlkDlg.h"
#include "OpNodeDlg.h"
#include "OpSegDlg.h"
#include "OpArcSegDlg.h"
#include "OpGrp.h"
#include "ArcDlg.h"

extern BOOL bLinehook;
extern lua_State *lua;
extern CFemmeApp theApp;

double CFemmeDoc::LineLength(int i)
{
	return abs(nodelist[linelist[i].n0].CC()-
		   nodelist[linelist[i].n1].CC());
}


BOOL CFemmeDoc::HasPeriodicBC()
{
	BOOL flag=FALSE;
	int i,j,k;

	for(i=0;i<lineproplist.GetSize();i++)
	{
		if ((lineproplist[i].BdryFormat==4) ||
			(lineproplist[i].BdryFormat==5)) 
			flag=TRUE;
	}
	// if flag is false, there can't be any lines
	// with periodic BC's, because no periodic boundary
	// conditions have been defined.
	if (flag==FALSE) return FALSE;

	//now, if there are some periodic boundary conditions,
	//we have to check to see if any have actually been
	//applied to the model
	flag=FALSE; // reset flag
	
	// first, test the segments
	for(i=0;i<linelist.GetSize();i++)
	{
		for(j=0,k=-1;j<lineproplist.GetSize();j++)
		{
			if(lineproplist[j].BdryName==
			   linelist[i].BoundaryMarker)
			{
				k=j;
				break;
			}
		}
		if(k>=0){
			if ((lineproplist[k].BdryFormat==4) ||
				(lineproplist[k].BdryFormat==5))
			{
				flag=TRUE;
				break;
			}
		}
	}

	if (flag==TRUE) return TRUE;

	// If we've gotten this far, we still need to check the
	// arc segments.
	for(i=0;i<arclist.GetSize();i++)
	{
		for(j=0,k=-1;j<lineproplist.GetSize();j++)
		{
			if(lineproplist[j].BdryName==
			   arclist[i].BoundaryMarker)
			{
				k=j;
				break;
			}
		}
		if(k>=0){
			if ((lineproplist[k].BdryFormat==4) ||
				(lineproplist[k].BdryFormat==5))
			{
				flag=TRUE;
				break;
			}
		}
	}

	// Finally, we're done. The value of flag now reflects
	// the judgement on whether or not we have periodic
	// and/or antiperiodic boundaries.
	return flag;
}



// What we do in the normal case is OnWritePoly
BOOL CFemmeDoc::OnWritePoly() 
{
	FILE *fp;
	int i,j,k,l,t;
	double z,R;
	CComplex a0,a1,a2,c;
	CString s;
	CArray< CNode, CNode&>             nodelst;
	CArray< CSegment, CSegment&>       linelst;
	CArray< CArcSegment, CArcSegment&> arclst;
	CArray< CBlockLabel, CBlockLabel&> blocklst;
	CNode node;
	CSegment segm;

	nodelst.RemoveAll();
	linelst.RemoveAll();

	// copy node list as it is;
	for(i=0;i<nodelist.GetSize();i++) nodelst.Add(nodelist[i]);

	// discretize input segments
	for(i=0;i<linelist.GetSize();i++)
	{
		if (linelist[i].MaxSideLength==-1) k=1;
		else{
			a0.Set(nodelist[linelist[i].n0].x,nodelist[linelist[i].n0].y);
			a1.Set(nodelist[linelist[i].n1].x,nodelist[linelist[i].n1].y);
			z=abs(a1-a0);
			k=(int) ceil(z/linelist[i].MaxSideLength);
		}

		if (k==1) linelst.Add(linelist[i]);
		else{
			segm=linelist[i];
			for(j=0;j<k;j++)
			{
				a2=a0+(a1-a0)*((double) (j+1))/((double) k);
				node.x=a2.re; node.y=a2.im;
				if(j==0){
					l=nodelst.GetSize();
					nodelst.Add(node);
					segm.n0=linelist[i].n0;
					segm.n1=l;
					linelst.Add(segm);
				}
				else if(j==(k-1))
				{
					l=nodelst.GetSize()-1;
					segm.n0=l;
					segm.n1=linelist[i].n1;
					linelst.Add(segm);
				}
				else{
					l=nodelst.GetSize();
					nodelst.Add(node);
					segm.n0=l-1;
					segm.n1=l;
					linelst.Add(segm);
				}
			}
		}
	}

	// discretize input arc segments
	for(i=0;i<arclist.GetSize();i++)
	{
		a2.Set(nodelist[arclist[i].n0].x,nodelist[arclist[i].n0].y);
		k=(int) ceil(arclist[i].ArcLength/arclist[i].MaxSideLength);
		segm.BoundaryMarker=arclist[i].BoundaryMarker;
		GetCircle(arclist[i],c,R);
		a1=exp(I*arclist[i].ArcLength*PI/(((double) k)*180.));

		if(k==1){
			segm.n0=arclist[i].n0;
			segm.n1=arclist[i].n1;
			linelst.Add(segm);
		}
		else for(j=0;j<k;j++)
		{
			a2=(a2-c)*a1+c;
			node.x=a2.re; node.y=a2.im;
			if(j==0){
				l=nodelst.GetSize();
				nodelst.Add(node);
				segm.n0=arclist[i].n0;
				segm.n1=l;
				linelst.Add(segm);
			}
			else if(j==(k-1))
			{
				l=nodelst.GetSize()-1;
				segm.n0=l;
				segm.n1=arclist[i].n1;
				linelst.Add(segm);
			}
			else{
				l=nodelst.GetSize();
				nodelst.Add(node);
				segm.n0=l-1;
				segm.n1=l;
				linelst.Add(segm);
			}
		}
	}


	// create correct output filename;
	CString pn = GetPathName();
	CString plyname=pn.Left(pn.ReverseFind('.')) + ".poly";
	
	// check to see if we are ready to write a datafile;
	
	if ((fp=fopen(plyname,"wt"))==NULL){
		AfxMessageBox("Couldn't write to specified .poly file");
		return FALSE;
	}
	
	// write out node list
	fprintf(fp,"%i	2	0	1\n",nodelst.GetSize());
	for(i=0;i<nodelst.GetSize();i++)
	{
		for(j=0,t=0;j<nodeproplist.GetSize();j++)
				if(nodeproplist[j].PointName==nodelst[i].BoundaryMarker) t=j+2;
		fprintf(fp,"%i	%.17g	%.17g	%i\n",i,nodelst[i].x,nodelst[i].y,t);
	}

	// write out segment list
	fprintf(fp,"%i	1\n",linelst.GetSize());
	for(i=0;i<linelst.GetSize();i++)
	{
		for(j=0,t=0;j<lineproplist.GetSize();j++)
				if(lineproplist[j].BdryName==linelst[i].BoundaryMarker) t=-(j+2);
		fprintf(fp,"%i	%i	%i	%i\n",i,linelst[i].n0,linelst[i].n1,t);
	}

	// write out list of holes;
	for(i=0,j=0;i<blocklist.GetSize();i++)
		if(blocklist[i].BlockType=="<No Mesh>") j++;
	fprintf(fp,"%i\n",j);
	for(i=0,k=0;i<blocklist.GetSize();i++)
		if(blocklist[i].BlockType=="<No Mesh>")
		{
			fprintf(fp,"%i	%.17g	%.17g\n",k,blocklist[i].x,blocklist[i].y);
			k++;
		}
	
	// write out regional attributes
	fprintf(fp,"%i\n",blocklist.GetSize()-j);
	for(i=0,k=0;i<blocklist.GetSize();i++)
		if(blocklist[i].BlockType!="<No Mesh>")
		{
			fprintf(fp,"%i	%.17g	%.17g	",k,blocklist[i].x,blocklist[i].y);
			fprintf(fp,"%i	",k+1);
			if (blocklist[i].MaxArea>0)
				fprintf(fp,"%.17g\n",blocklist[i].MaxArea);
			else fprintf(fp,"-1\n");
			k++;
		}
	fclose(fp);

	// write out a trivial pbc file
	plyname=pn.Left(pn.ReverseFind('.')) + ".pbc";
	if ((fp=fopen(plyname,"wt"))==NULL){
		AfxMessageBox("Couldn't write to specified .pbc file");
		return FALSE;
	}
	fprintf(fp,"0\n");
	fclose(fp);

	//call triangle

	CString rootname="\"" + pn.Left(pn.ReverseFind('.')) + "\"";
	char CommandLine[512];
	sprintf(CommandLine,"\"%striangle.exe\" -p -P -q30 -e -A -a -z -Q -I %s",
		BinDir,rootname);


	STARTUPINFO StartupInfo = {0};
	PROCESS_INFORMATION ProcessInfo;
	StartupInfo.cb = sizeof(STARTUPINFO);
	StartupInfo.dwFlags = STARTF_USESHOWWINDOW;
	StartupInfo.wShowWindow = SW_MINIMIZE;
	if (CreateProcess(NULL,CommandLine, NULL, NULL, FALSE,
		0, NULL, NULL, &StartupInfo, &ProcessInfo)){

		if(bLinehook==FALSE) WaitForSingleObject(ProcessInfo.hProcess, INFINITE);
		else{
			DWORD ExitCode;
			hProc=ProcessInfo.hProcess;
			do{
				GetExitCodeProcess(ProcessInfo.hProcess,&ExitCode);
				line_hook(lua,NULL);	
			} while(ExitCode==STILL_ACTIVE);
			hProc=NULL;
		}
	
	}
	else
	{
		AfxMessageBox("Couldn't spawn triangle.exe");
		return FALSE;
	}
	
	DWORD ExitCode;
	GetExitCodeProcess(
		ProcessInfo.hProcess,	// handle to the process 
		&ExitCode 				// address to receive termination status 
	);

	if (ExitCode!=0)
	{
		AfxMessageBox("Call to triangle was unsuccessful");
		return FALSE;
	}



	return TRUE;
}


// Call triangle to order segments on the boundary properly
BOOL CFemmeDoc::FunnyOnWritePoly() 
{
	FILE *fp;
	int i,j,k,l,t,n,n0,n1,n2;
	double z,R;
	CComplex a0,a1,a2,c;
	CComplex b0,b1,b2;
	char instring[1024];
	CString s;
	CArray< CNode, CNode&>             nodelst;
	CArray< CSegment, CSegment&>       linelst;
	CArray< CArcSegment, CArcSegment&> arclst;
	CArray< CBlockLabel, CBlockLabel&> blocklst;
	CArray< CPeriodicBoundary, CPeriodicBoundary&> pbclst;
	CArray< CCommonPoint, CCommonPoint& >ptlst;
	CNode node;
	CSegment segm;
	CPeriodicBoundary pbc;
	CCommonPoint pt;

	nodelst.RemoveAll();
	linelst.RemoveAll();
	pbclst.RemoveAll();
	ptlst.RemoveAll();

	UpdateUndo();

	// copy node list as it is;
	for(i=0;i<nodelist.GetSize();i++) nodelst.Add(nodelist[i]);

	// discretize input segments
	for(i=0;i<linelist.GetSize();i++)
	{
		// abuse the IsSelected flag to carry a notation
		// of which line or arc in the input geometry a
		// particular segment is associated with
		segm=linelist[i];
		segm.IsSelected=i;

		if (linelist[i].MaxSideLength==-1) k=1;
		else{
			a0.Set(nodelist[linelist[i].n0].x,nodelist[linelist[i].n0].y);
			a1.Set(nodelist[linelist[i].n1].x,nodelist[linelist[i].n1].y);
			z=abs(a1-a0);
			k=(int) ceil(z/linelist[i].MaxSideLength);
		}

		if (k==1) linelst.Add(segm);
		else{
			for(j=0;j<k;j++)
			{
				a2=a0+(a1-a0)*((double) (j+1))/((double) k);
				node.x=a2.re; node.y=a2.im;
				if(j==0){
					l=nodelst.GetSize();
					nodelst.Add(node);
					segm.n0=linelist[i].n0;
					segm.n1=l;
					linelst.Add(segm);
				}
				else if(j==(k-1))
				{
					l=nodelst.GetSize()-1;
					segm.n0=l;
					segm.n1=linelist[i].n1;
					linelst.Add(segm);
				}
				else{
					l=nodelst.GetSize();
					nodelst.Add(node);
					segm.n0=l-1;
					segm.n1=l;
					linelst.Add(segm);
				}
			}
		}
	}

	// discretize input arc segments
	for(i=0;i<arclist.GetSize();i++)
	{
		segm.IsSelected=i+linelist.GetSize();
		a2.Set(nodelist[arclist[i].n0].x,nodelist[arclist[i].n0].y);
		k=(int) ceil(arclist[i].ArcLength/arclist[i].MaxSideLength);
		segm.BoundaryMarker=arclist[i].BoundaryMarker;
		GetCircle(arclist[i],c,R);
		a1=exp(I*arclist[i].ArcLength*PI/(((double) k)*180.));

		if(k==1){
			segm.n0=arclist[i].n0;
			segm.n1=arclist[i].n1;
			linelst.Add(segm);
		}
		else for(j=0;j<k;j++)
		{
			a2=(a2-c)*a1+c;
			node.x=a2.re; node.y=a2.im;
			if(j==0){
				l=nodelst.GetSize();
				nodelst.Add(node);
				segm.n0=arclist[i].n0;
				segm.n1=l;
				linelst.Add(segm);
			}
			else if(j==(k-1))
			{
				l=nodelst.GetSize()-1;
				segm.n0=l;
				segm.n1=arclist[i].n1;
				linelst.Add(segm);
			}
			else{
				l=nodelst.GetSize();
				nodelst.Add(node);
				segm.n0=l-1;
				segm.n1=l;
				linelst.Add(segm);