writepoly.cpp

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

				segm=linelist[s0];
				for(j=0;j<k;j++)
				{
					a2=a0+(a1-a0)*((double) (j+1))/((double) k);
					b2=b0+(b1-b0)*((double) (j+1))/((double) k);
					node0.x=a2.re; node0.y=a2.im;
					node1.x=b2.re; node1.y=b2.im;
					if(j==0){
						l=nodelst.GetSize();
						nodelst.Add(node0);
						segm.n0=linelist[s0].n0;
						segm.n1=l;
						linelst.Add(segm);
						pt.x=l;

						l=nodelst.GetSize();
						nodelst.Add(node1);
						segm.n0=linelist[s1].n0;
						segm.n1=l;
						linelst.Add(segm);
						pt.y=l;

						pt.t=pbclst[n].BdryFormat;
						ptlst.Add(pt);
					}
					else if(j==(k-1))
					{
						// last subdivision--no ptlst
						// entry associated with this one.
						l=nodelst.GetSize()-2;
						segm.n0=l;
						segm.n1=linelist[s0].n1;
						linelst.Add(segm);

						l=nodelst.GetSize()-1;
						segm.n0=l;
						segm.n1=linelist[s1].n1;
						linelst.Add(segm);
					}
					else{
						l=nodelst.GetSize();
						
						nodelst.Add(node0);
						nodelst.Add(node1);

						segm.n0=l-2;
						segm.n1=l;
						linelst.Add(segm);

						segm.n0=l-1;
						segm.n1=l+1;
						linelst.Add(segm);
						
						pt.x=l;
						pt.y=l+1;
						pt.t=pbclst[n].BdryFormat;
						ptlst.Add(pt);
					}
				}
			}
		}
		else{  // if this pbc is an arc segment...
		
			int s0,s1;
			int p0[2],p1[2];
			CNode node0,node1;
			CComplex bgn0,bgn1,c0,c1,d0,d1;
			double r0,r1;
			
			s0=pbclst[n].seg[0];
			s1=pbclst[n].seg[1];
			arclist[s0].IsSelected=TRUE;
			arclist[s1].IsSelected=TRUE;

			k=(int) ceil(arclist[s0].ArcLength/arclist[s0].MaxSideLength);
			segm.BoundaryMarker=arclist[s0].BoundaryMarker;
			GetCircle(arclist[s0],c0,r0);
			GetCircle(arclist[s1],c1,r1);

			if (arclist[s0].NormalDirection==0){ 
				bgn0=nodelist[arclist[s0].n0].CC();
				d0=exp(I*arclist[s0].ArcLength*PI/(((double) k)*180.));
				p0[0]=arclist[s0].n0;
				p0[1]=arclist[s0].n1;
			}
			else{
				bgn0=nodelist[arclist[s0].n1].CC();
				d0=exp(-I*arclist[s0].ArcLength*PI/(((double) k)*180.));
				p0[0]=arclist[s0].n1;
				p0[1]=arclist[s0].n0;
			}
	
			if (arclist[s1].NormalDirection!=0){ 
				bgn1=nodelist[arclist[s1].n0].CC();
				d1=exp(I*arclist[s1].ArcLength*PI/(((double) k)*180.));
				p1[0]=arclist[s1].n0;
				p1[1]=arclist[s1].n1;
			}
			else{
				bgn1=nodelist[arclist[s1].n1].CC();
				d1=exp(-I*arclist[s1].ArcLength*PI/(((double) k)*180.));
				p1[0]=arclist[s1].n1;
				p1[1]=arclist[s1].n0;
			}

			// add arc segment end points to the list;
			pt.x=p0[0]; pt.y=p1[0]; pt.t=pbclst[n].BdryFormat;
			ptlst.Add(pt);
			pt.x=p0[1]; pt.y=p1[1]; pt.t=pbclst[n].BdryFormat;
			ptlst.Add(pt);
	
			if (k==1){

				// catch the case in which the line
				// doesn't get subdivided.
				segm.n0=p0[0]; segm.n1=p0[1];
				linelst.Add(segm);
				segm.n0=p1[0]; segm.n1=p1[1];
				linelst.Add(segm);
			}
			else{
				for(j=0;j<k;j++)
				{
					bgn0=(bgn0-c0)*d0+c0;
					node0.x=bgn0.re; node0.y=bgn0.im;

					bgn1=(bgn1-c1)*d1+c1;
					node1.x=bgn1.re; node1.y=bgn1.im;

					if(j==0){
						l=nodelst.GetSize();
						nodelst.Add(node0);
						segm.n0=p0[0];
						segm.n1=l;
						linelst.Add(segm);
						pt.x=l;

						l=nodelst.GetSize();
						nodelst.Add(node1);
						segm.n0=p1[0];
						segm.n1=l;
						linelst.Add(segm);
						pt.y=l;

						pt.t=pbclst[n].BdryFormat;
						ptlst.Add(pt);
					}
					else if(j==(k-1))
					{
						// last subdivision--no ptlst
						// entry associated with this one.
						l=nodelst.GetSize()-2;
						segm.n0=l;
						segm.n1=p0[1];
						linelst.Add(segm);

						l=nodelst.GetSize()-1;
						segm.n0=l;
						segm.n1=p1[1];
						linelst.Add(segm);
					}
					else{
						l=nodelst.GetSize();
						
						nodelst.Add(node0);
						nodelst.Add(node1);

						segm.n0=l-2;
						segm.n1=l;
						linelst.Add(segm);

						segm.n0=l-1;
						segm.n1=l+1;
						linelst.Add(segm);
						
						pt.x=l;
						pt.y=l+1;
						pt.t=pbclst[n].BdryFormat;
						ptlst.Add(pt);
					}
				}
				
			}
		}
	}
	
	// Then, do the rest of the lines and arcs in the
	// "normal" way and write .poly file.

	// discretize input segments
	for(i=0;i<linelist.GetSize();i++)
	if(linelist[i].IsSelected==FALSE){
		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++)
	if(arclist[i].IsSelected==FALSE){
		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;
	pn = GetPathName();
	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");
		Undo();  UnselectAll();
		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);

	// Make sure to prune out any duplications in the ptlst	
	for(k=0;k<ptlst.GetSize();k++) ptlst[k].Order();
	k=0;
	while((k+1) < ptlst.GetSize())
	{
		j=k+1;
		while(j < ptlst.GetSize())
		{
			if((ptlst[k].x==ptlst[j].x) && (ptlst[k].y==ptlst[j].y))
				ptlst.RemoveAt(j);
			else j++;
		}
		k++;
	}

	// used to have a check to eliminate the case where a point
	// and its companion are the same point--actually, this shouldn't
	// be a problem just to let the algorithm deal with this
	// as usual.

	// One last error check--each point must have only one companion point.
	// however, it would be possible to screw up in the definition of the BCs
	// so that this isn't the case.  Look through the points to try and catch
	// this one.
/*
	// let's let this check go away for a minute...

	for(k=0,n=FALSE;(k+1)<ptlst.GetSize();k++)
	{
		for(j=k+1;j<ptlst.GetSize();j++)
		{
			if(ptlst[k].x==ptlst[j].x) n=TRUE;
			if(ptlst[k].y==ptlst[j].y) n=TRUE;
			if(ptlst[k].x==ptlst[j].y) n=TRUE;
			if(ptlst[k].y==ptlst[j].x) n=TRUE;
		}
	}
	if (n==TRUE){
		AfxMessageBox("Nonphysical (anti)periodic boundary assignments");
		Undo();  UnselectAll();
		return FALSE;
	}
*/
	// write out a pbc file containing a list of linked nodes
	plyname=pn.Left(pn.ReverseFind('.')) + ".pbc";
	if ((fp=fopen(plyname,"wt"))==NULL){
		AfxMessageBox("Couldn't write to specified .pbc file");
		Undo();  UnselectAll();
		return FALSE;
	}
	fprintf(fp,"%i\n",ptlst.GetSize());
	for(k=0;k<ptlst.GetSize();k++) 
		fprintf(fp,"%i	%i	%i	%i\n",k,ptlst[k].x,ptlst[k].y,ptlst[k].t);
	fclose(fp);

	// call triangle with -Y flag.

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

	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");
		Undo();  UnselectAll();
		return FALSE;
	}
	
	GetExitCodeProcess(
		ProcessInfo.hProcess,	// handle to the process 
		&ExitCode 				// address to receive termination status 
	);

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

	UnselectAll();

	// Now restore boundary segment discretizations that have
	// been mucked up in the process...
	for(i=0;i<linelist.GetSize();i++)
		linelist[i].MaxSideLength=undolinelist[i].MaxSideLength;
	
	// and save the latest version of the document to make sure
	// any changes to arc discretization get propagated into
	// the solution description....
	OnSaveDocument(pn);

	return TRUE;
}