writepoly.cpp

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

			}
		}
	}


	// 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");
		Undo();  UnselectAll();
		return FALSE;
	}
	
	// write out node list
	fprintf(fp,"%i	2	0	1\n",nodelst.GetSize());
	for(i=0;i<nodelst.GetSize();i++)
	{
		fprintf(fp,"%i	%.17g	%.17g	%i\n",
			     i,nodelst[i].x,nodelst[i].y,0);
	}

	// write out segment list
	fprintf(fp,"%i	1\n",linelst.GetSize());
	for(i=0;i<linelst.GetSize();i++)
	{
		t=-(linelst[i].IsSelected+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);

	//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");
		Undo();  UnselectAll();
		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");
		Undo();  UnselectAll();
		return FALSE;
	}

	// So far, so good.  Now, read back in the .edge file
	// to make sure the points in the segments and arc
	// segments are ordered in a consistent way so that
	// the (anti)periodic boundary conditions can be applied.

	//read meshlines;
	plyname=pn.Left(pn.ReverseFind('.')) + ".edge";
	if((fp=fopen(plyname,"rt"))==NULL){
		fclose(fp);
		AfxMessageBox("Couldn't open .edge file");
		Undo();  UnselectAll();
		return FALSE;
	}
	fgets(instring,1024,fp);
	sscanf(instring,"%i",&k);
	UnselectAll();	// abuse IsSelected again to keep a
					// tally of how many subsegments each
					// entity is sliced into.
	
	ptlst.SetSize(linelist.GetSize()+arclist.GetSize());
	for(i=0;i<ptlst.GetSize();i++) ptlst[i].t=0;

	for(i=0;i<k;i++)
	{
		fgets(instring,1024,fp);
		sscanf(instring,"%i	%i	%i	%i",&l,&n0,&n1,&j);
		if(j!=0)
		{
			j=-(j+2); // convert back to the `right' numbering

			// store a reference line that we can use to
			// determine whether or not this is a
			// boundary segment w/out re-running triangle.
			if (ptlst[j].t==0)
			{
				ptlst[j].t=1;
				if(n0<n1){
					ptlst[j].x=n0;
					ptlst[j].y=n1;
				}
				else{
					ptlst[j].x=n1;
					ptlst[j].y=n0;
				}
			}

			if(j<linelist.GetSize())
			{
				// deal with segments
				linelist[j].IsSelected++;
				
				if((linelist[j].n0==n1) || (linelist[j].n1==n0))
				{
					t=linelist[j].n0;
					linelist[j].n0=linelist[j].n1;
					linelist[j].n1=t;
				}
			}
			else{
				// deal with arc segments;
				// Can't just flip the point order with
				// impunity in the arc segments, so we flip
				// a marker which denotes which side the
				// normal is on.
			
				j=j-linelist.GetSize();
				arclist[j].IsSelected++;
				if((arclist[j].n0==n1) || (arclist[j].n1==n0))
					arclist[j].NormalDirection=FALSE;
				if((arclist[j].n0==n0) || (arclist[j].n1==n1))
					arclist[j].NormalDirection=TRUE;
			}
		}
	}
	fclose(fp);

	// figure out which segments / arcsegments are on the
	// boundary and force an appropriate mesh density on 
	// these based on how many divisions are in the first
	// trial meshing of the domain.

	// paw through the element list to find out how many
	// elements each reference segment appears in.  If a 
	// segment is on the boundary, it ought to appear in just
	// one element.  Otherwise, it appears in two.
	plyname=pn.Left(pn.ReverseFind('.')) + ".ele";
	if((fp=fopen(plyname,"rt"))==NULL){
		fclose(fp);
		AfxMessageBox("Couldn't open .ele file");
		Undo();  UnselectAll();
		return FALSE;
	}
	fgets(instring,1024,fp);
	sscanf(instring,"%i",&k);
	for(i=0;i<k;i++)
	{
		fgets(instring,1024,fp);
		sscanf(instring,"%i	%i	%i	%i",&j,&n0,&n1,&n2);
		// Sort out the three nodes...
		if (n0>n1) { n=n0; n0=n1; n1=n; }
		if (n1>n2) { n=n1; n1=n2; n2=n; }
		if (n0>n1) { n=n0; n0=n1; n1=n; }

		// now, check to see if any of the test segments
		// are sides of this node...
		for(j=0;j<ptlst.GetSize();j++)
		{
			if ((n0==ptlst[j].x) && (n1==ptlst[j].y)) ptlst[j].t--;
			if ((n0==ptlst[j].x) && (n2==ptlst[j].y)) ptlst[j].t--;
			if ((n1==ptlst[j].x) && (n2==ptlst[j].y)) ptlst[j].t--;
		}
	}
	fclose(fp);

	// impose "new" mesh constraints on bdry arcs and segments....
	for(i=0;i<linelist.GetSize();i++)
	{
		if (ptlst[i].t==0) linelist[i].MaxSideLength=
			LineLength(i)/((double) linelist[i].IsSelected);
	}
	for(i=0;i<arclist.GetSize();i++)
	{
		if (ptlst[i+linelist.GetSize()].t==0){
			// alter maxsidelength, but do it in such
			// a way that it carries only 4 significant
			// digits.  There's no use in carrying double
			// precision here, because it looks crappy
			// when you open up the arc segment to see
			// its properties.
			char kludge[32];
			arclist[i].MaxSideLength=
			arclist[i].ArcLength/((double) arclist[i].IsSelected);
			sprintf(kludge,"%.4g",arclist[i].MaxSideLength);
			sscanf(kludge,"%lf",&arclist[i].MaxSideLength);
		}
	}
	ptlst.RemoveAll();	

		// want to impose explicit discretization only on
		// the boundary arcs and segments.  After the meshing
		// is done, spacing on boundary segments should be
		// restored to the value that was there before meshing
		// was called, but the arc segments should keep the
		// "new" MaxSideLength--this is used in other places
		// and must always be consistent with the the mesh.

	
	// Now, do a shitload of checking to make sure that
	// the PBCs haven't been defined by the user
	// in a messed up way.

	// First, search through defined bc's for periodic ones;
	for(i=0;i<lineproplist.GetSize();i++)
	{
		if ((lineproplist[i].BdryFormat==4) ||
			(lineproplist[i].BdryFormat==5)){
			pbc.BdryName=lineproplist[i].BdryName;
			pbc.BdryFormat=lineproplist[i].BdryFormat-4; // 0 for pbc, 1 for apbc
			pbclst.Add(pbc);
		}
	}

	for(i=0;i<linelist.GetSize();i++)
	{
		for(j=0;j<pbclst.GetSize();j++)
		{
			if (pbclst[j].BdryName==linelist[i].BoundaryMarker)
			{
				// A pbc or apbc can only be applied to 2 segs
				// at a time.  If it is applied to multiple segs
				// at the same time, flag it and kick it out.
				if (pbclst[j].nseg==2)
				{
					AfxMessageBox("An (anti)periodic BC is assigned to more than two segments");
					Undo();  UnselectAll();
					return FALSE;
				}
				pbclst[j].seg[pbclst[j].nseg]=i;
				pbclst[j].nseg++;
			}
		}
	}
	
	for(i=0;i<arclist.GetSize();i++)
	{
		for(j=0;j<pbclst.GetSize();j++)
		{
			if (pbclst[j].BdryName==arclist[i].BoundaryMarker)
			{
				// A pbc or apbc can only be applied to 2 arcs
				// at a time.  If it is applied to multiple arcs
				// at the same time, flag it and kick it out.
				if (pbclst[j].narc==2)
				{
					AfxMessageBox("An (anti)periodic BC is assigned to more than two arcs");
					Undo();  UnselectAll();
					return FALSE;
				}
				pbclst[j].seg[pbclst[j].narc]=i;
				pbclst[j].narc++;
			}
		}
	}

	j=0;
	while(j<pbclst.GetSize())
	{
		// check for a bc that is a mix of arcs and segments.
		// this is an error, and it should get flagged.
		if ((pbclst[j].nseg>0) && (pbclst[j].narc>0))
		{
			AfxMessageBox("Can't mix arcs and segments for (anti)periodic BCs");
			Undo();  UnselectAll();
			return FALSE;
		}

		
		// remove any periodic BC's that aren't actually in play
		if((pbclst[j].nseg<2) && (pbclst[j].narc<2)) pbclst.RemoveAt(j);
		else j++;
	}

	for(j=0;j<pbclst.GetSize();j++)
	{
		// check to see if adjoining entries are applied
		// to objects of compatible size/shape, and
		// reconcile meshing on the objects.
		
		// for segments:
		if(pbclst[j].nseg>0){
			
			// make sure that lines are pretty much the same length
			if(fabs(LineLength(pbclst[j].seg[0])
		           -LineLength(pbclst[j].seg[1]))>1.e-06)
			{
				AfxMessageBox("(anti)periodic BCs applied to dissimilar segments");
				Undo();  UnselectAll();
				return FALSE;
			}
			
			// make sure that both lines have the same spacing
			double len1,len2,len;
			len1=linelist[pbclst[j].seg[0]].MaxSideLength;
			len2=linelist[pbclst[j].seg[1]].MaxSideLength;
			
			if(len1<=0) len1=len2;
			if(len2<=0) len2=len1;
			len=min(len1,len2);

			linelist[pbclst[j].seg[0]].MaxSideLength=len;
			linelist[pbclst[j].seg[1]].MaxSideLength=len;
		}

		// for arc segments:
		if(pbclst[j].narc>0){
			
			// make sure that arcs are pretty much the 
			// same arc length
			if(fabs(arclist[pbclst[j].seg[0]].ArcLength
		           -arclist[pbclst[j].seg[1]].ArcLength)>1.e-06)
			{
				AfxMessageBox("(anti)periodic BCs applied to dissimilar arc segments");
				Undo();  UnselectAll();
				return FALSE;
			}

			// make sure that both lines have the same spacing
			double len1,len2,len;
			len1=arclist[pbclst[j].seg[0]].MaxSideLength;
			len2=arclist[pbclst[j].seg[1]].MaxSideLength;
	
			len=min(len1,len2);

			arclist[pbclst[j].seg[0]].MaxSideLength=len;
			arclist[pbclst[j].seg[1]].MaxSideLength=len;
		}
	}

	// write out new poly and write out adjacent
	// boundary nodes in a separate .pbc file.
	
	// kludge things a bit and use IsSelected to denote
	// whether or not a line or arc has already been processed.
	UnselectAll();
	nodelst.RemoveAll();
	linelst.RemoveAll();

	// first, add in existing nodes
	for(n=0;n<nodelist.GetSize();n++) nodelst.Add(nodelist[n]);

	for(n=0;n<pbclst.GetSize();n++)
	{
		if (pbclst[n].nseg!=0) // if this pbc is a line segment...
		{
			int s0,s1;
			CNode node0,node1;

			s0=pbclst[n].seg[0];
			s1=pbclst[n].seg[1];
			linelist[s0].IsSelected=TRUE;
			linelist[s1].IsSelected=TRUE;
			
			// make is so that first point on first line
			// maps to first point on second line...
			t=linelist[s1].n1;
			linelist[s1].n1=linelist[s1].n0;
			linelist[s1].n0=t;

			// store number of sub-segments in k
			if (linelist[s0].MaxSideLength==-1) k=1;
			else{
				a0=nodelist[linelist[s0].n0].CC();
				a1=nodelist[linelist[s0].n1].CC();
				b0=nodelist[linelist[s1].n0].CC();
				b1=nodelist[linelist[s1].n1].CC();
				z=abs(a1-a0);
				k=(int) ceil(z/linelist[s0].MaxSideLength);
			}

			// add segment end points to the list;
			pt.x=linelist[s0].n0;
			pt.y=linelist[s1].n0;
			pt.t=pbclst[n].BdryFormat;
			ptlst.Add(pt);
			pt.x=linelist[s0].n1;
			pt.y=linelist[s1].n1;
			pt.t=pbclst[n].BdryFormat;
			ptlst.Add(pt);
	
			if (k==1){
				// catch the case in which the line
				// doesn't get subdivided.
				linelst.Add(linelist[s0]);
				linelst.Add(linelist[s1]);
			}
			else{