femmedoc.cpp

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

// femmeDoc.cpp : implementation of the CFemmeDoc class
//

#include "stdafx.h"
#include "femme.h"
#include "femmeDoc.h"
#include "femmeView.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"
#include "ExteriorProps.h"

extern void * pFemmdoc;
extern lua_State *lua;
extern BOOL bLinehook;
extern CFemmeApp theApp;
extern CString luascriptname;

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CFemmeDoc

IMPLEMENT_DYNCREATE(CFemmeDoc, CDocument)

BEGIN_MESSAGE_MAP(CFemmeDoc, CDocument)
	//{{AFX_MSG_MAP(CFemmeDoc)
	ON_COMMAND(ID_DEFINE_PROBLEM, OnDefineProblem)
	ON_COMMAND(ID_EDIT_MATPROPS, OnEditMatprops)
	ON_COMMAND(ID_EDIT_PTPROPS, OnEditPtprops)
	ON_COMMAND(ID_EDIT_SEGPROPS, OnEditSegprops)
	ON_COMMAND(ID_EDIT_CIRCPROPS, OnEditCircprops)
	ON_COMMAND(ID_EDIT_EXTERIOR, OnEditExterior)
	ON_COMMAND(ID_VIEW_LUACONSOLE, OnViewLuaConsole)
	ON_COMMAND(ID_FILE_OPEN_LUA_SCRIPT, OnFileOpenLuaScript)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

CFemmeDoc::CFemmeDoc()
{
	// set some default values for starting up rendering
	// things properly
	FirstDraw=FALSE;
	NoDraw=FALSE;

	// set up some default document behaviors
	d_prec=1.e-08;
	d_minangle=30.;
	d_freq=0;
	d_depth=1;
	d_coord=0;
	d_length=0;
	d_type=0;
	d_luaconsole=FALSE;

	// Figure out what directory the executables
	// are in, so we can call `triangle' if we need to.
	BinDir=theApp.GetExecutablePath();

	// read document default behaviors from disk
	ScanPreferences();

	// fire up lua
	initalise_lua();
		
	// initialize the data in the document
	OnNewDocument();
}

BOOL CFemmeDoc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// clear out all current lines, nodes, and block labels
	nodelist.RemoveAll();
	linelist.RemoveAll();
	arclist.RemoveAll();
	blocklist.RemoveAll();
	undonodelist.RemoveAll();
	undolinelist.RemoveAll();
	undoarclist.RemoveAll();
	undoblocklist.RemoveAll();
	nodeproplist.RemoveAll();
	lineproplist.RemoveAll();
	blockproplist.RemoveAll();
	circproplist.RemoveAll();
	meshnode.RemoveAll();
	meshline.RemoveAll();
	greymeshline.RemoveAll();

	// set problem attributes to generic ones;
	Frequency=d_freq;
	Precision=d_prec;
	MinAngle=d_minangle;
	Depth=d_depth;
	LengthUnits=d_length;
	ProblemType=d_type;
	Coords=d_coord;
	ProblemNote="Add comments here.";
	extRo=extRi=extZo=0;

	// reset view to default attributes
	CFemmeView *pView;
	POSITION pos;
	pos=GetFirstViewPosition();
	if (pos!=NULL){	
		pView=(CFemmeView *)GetNextView(pos);
//		if (pView!=NULL) pView->OnNewDocument();
	}

	return TRUE;
}

CFemmeDoc::~CFemmeDoc()
{
//	pFemmeDoc=NULL;
}

/////////////////////////////////////////////////////////////////////////////
// CFemmeDoc diagnostics

#ifdef _DEBUG
void CFemmeDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CFemmeDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CFemmeDoc serialization

void CFemmeDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CFemmeDoc commands

void CFemmeDoc::OnDefineProblem() 
{
	probdlg pDlg;
	
	// Send present parameter values to the dialog
	// pDlg.m_rji=TheDoc->vi[0];
	pDlg.probtype       = ProblemType;
	pDlg.m_problem_note = ProblemNote;
	pDlg.m_frequency    = Frequency;
	pDlg.m_precision	= Precision;
	pDlg.m_minangle		= MinAngle;
	pDlg.m_depth		= Depth;
	pDlg.lengthunits	= LengthUnits;

	// Display dialog and collect data
	if(pDlg.DoModal()==IDOK)
	{
		Frequency   = pDlg.m_frequency;
		Precision	= pDlg.m_precision;
		MinAngle	= pDlg.m_minangle;
		ProblemNote = pDlg.m_problem_note;
		ProblemType = pDlg.probtype;
		LengthUnits = pDlg.lengthunits;
		Depth		= pDlg.m_depth;
	}
}

void CFemmeDoc::UnselectAll()
{
	int i;

	for(i=0;i < nodelist.GetSize();i++) nodelist[i].IsSelected=FALSE;
	for(i=0;i < linelist.GetSize();i++) linelist[i].IsSelected=FALSE;
	for(i=0;i < blocklist.GetSize();i++) blocklist[i].IsSelected=FALSE;
	for(i=0;i < arclist.GetSize();i++) arclist[i].IsSelected=FALSE;
}

BOOL CFemmeDoc::AddNode(double x, double y, double d)
{	
	int i,k;
	CNode pt;
	CSegment segm;
	CArcSegment asegm;
	CComplex c,a0,a1,a2;
	double R;

	// test to see if ``too close'' to existing node...
	for (i=0;i<nodelist.GetSize();i++)
		if(nodelist[i].GetDistance(x,y)<d) return FALSE;

	// can't put a node on top of a block label; do same sort of test.
	for (i=0;i<blocklist.GetSize();i++)
		if(blocklist[i].GetDistance(x,y)<d) return FALSE;

	// if all is OK, add point in to the node list...
	pt.x=x; pt.y=y;
	nodelist.Add(pt);

	// test to see if node is on an existing line; if so, 
	// break into two lines;
	k=linelist.GetSize();
	for(i=0;i<k;i++)
	{
		if (fabs(ShortestDistance(x,y,i))<d)
		{
			segm=linelist[i];
			linelist[i].n1=nodelist.GetSize()-1;
			segm.n0=nodelist.GetSize()-1;
			linelist.Add(segm);
		}
	}

	// test to see if node is on an existing arc; if so, 
	// break into two arcs;
	k=arclist.GetSize();
	for(i=0;i<k;i++)
	{
		if (ShortestDistanceFromArc(CComplex(x,y),arclist[i])<d)
		{
			a0.Set(nodelist[arclist[i].n0].x,nodelist[arclist[i].n0].y);
			a1.Set(nodelist[arclist[i].n1].x,nodelist[arclist[i].n1].y);
			a2.Set(x,y);
			GetCircle(arclist[i],c,R);
			asegm=arclist[i];
			arclist[i].n1=nodelist.GetSize()-1;
			arclist[i].ArcLength=arg((a2-c)/(a0-c))*180./PI;
			asegm.n0=nodelist.GetSize()-1;
			asegm.ArcLength=arg((a1-c)/(a2-c))*180./PI;
			arclist.Add(asegm);
		}
	}
	return TRUE;
}

BOOL CFemmeDoc::AddSegment(int n0, int n1, double tol)
{
	return AddSegment(n0,n1,NULL,tol);
}

BOOL CFemmeDoc::AddSegment(int n0, int n1, CSegment *parsegm, double tol)
{
	int i,j,k;
	double xi,yi,t;
	CComplex p[2];
	CSegment segm;
	CArray< CComplex, CComplex&> newnodes;

	newnodes.RemoveAll();

	// don't add if line is degenerate
	if (n0==n1) return FALSE;
	
	// don't add if the line is already in the list;
	for(i=0;i<linelist.GetSize();i++){
		if ((linelist[i].n0==n0) && (linelist[i].n1==n1)) return FALSE;
		if ((linelist[i].n0==n1) && (linelist[i].n1==n0)) return FALSE;
	}

	// add proposed line to the linelist
	 segm.BoundaryMarker="<None>";
	if (parsegm!=NULL) segm=*parsegm;
	segm.IsSelected=FALSE;
	segm.n0=n0; segm.n1=n1; 
	
	// check to see if there are intersections with segments
	for(i=0;i<linelist.GetSize();i++)
		if(GetIntersection(n0,n1,i,&xi,&yi)==TRUE)  newnodes.Add(CComplex(xi,yi));

	// check to see if there are intersections with arcs
	for(i=0;i<arclist.GetSize();i++){
		j=GetLineArcIntersection(segm,arclist[i],p);
		if (j>0) for(k=0;k<j;k++) newnodes.Add(p[k]);
	}

	// add nodes at intersections
	if (tol==0)
	{
		if (nodelist.GetSize()<2) t=1.e-08;
		else{
			CComplex p0,p1;
			p0=nodelist[0].CC();
			p1=p0;
			for(i=1;i<nodelist.GetSize();i++)
			{
				if(nodelist[i].x<p0.re) p0.re=nodelist[i].x;
				if(nodelist[i].x>p1.re) p1.re=nodelist[i].x;
				if(nodelist[i].y<p0.im) p0.im=nodelist[i].y;
				if(nodelist[i].y>p1.im) p1.im=nodelist[i].y;
			}
			t=abs(p1-p0)*1.e-06;
		}
	}
	else t=tol;

	for(i=0;i<newnodes.GetSize();i++) 
		AddNode(newnodes[i].re,newnodes[i].im,t);
	
	// Add proposed line segment
	linelist.Add(segm);

	// check to see if proposed line passes through other points;
    // if so, delete line and create lines that link intermediate points;
    // does this by recursive use of AddSegment;
	double d,dmin;
    UnselectAll();
    if (tol==0) dmin=abs(nodelist[n1].CC()-nodelist[n0].CC())*1.e-05;
	else dmin=tol;

    k=linelist.GetSize()-1;
    for(i=0;i<nodelist.GetSize();i++)
    {
        if( (i!=n0) && (i!=n1) )
        {
            d=ShortestDistance(nodelist[i].x,nodelist[i].y,k);
			if (abs(nodelist[i].CC()-nodelist[n0].CC())<dmin) d=2.*dmin;
			if (abs(nodelist[i].CC()-nodelist[n1].CC())<dmin) d=2.*dmin;
            if (d<dmin){
                linelist[k].ToggleSelect();
                DeleteSelectedSegments();
				if(parsegm==NULL)
				{
					AddSegment(n0,i,dmin);
					AddSegment(i,n1,dmin);
				}
				else{
					AddSegment(n0,i,&segm,dmin);
					AddSegment(i,n1,&segm,dmin);
				}
                i=nodelist.GetSize();
            }
        }
    }

	return TRUE;
}

void CFemmeDoc::GetCircle(CArcSegment &arc,CComplex &c, double &R)
{
	CComplex a0,a1,t;
	double d,tta;

	a0.Set(nodelist[arc.n0].x,nodelist[arc.n0].y);
	a1.Set(nodelist[arc.n1].x,nodelist[arc.n1].y);
	d=abs(a1-a0);			// distance between arc endpoints

	// figure out what the radius of the circle is...
	t=(a1-a0)/d;
	tta=arc.ArcLength*PI/180.;
	R=d/(2.*sin(tta/2.));
	c=a0 + (d/2. + I*sqrt(R*R-d*d/4.))*t; // center of the arc segment's circle...
}

int CFemmeDoc::GetLineArcIntersection(CSegment &seg, CArcSegment &arc, CComplex *p)
{
	CComplex p0,p1,a0,a1,t,v,c;
	double d,l,R,z,tta;
	int i=0;

	p0.Set(nodelist[seg.n0].x,nodelist[seg.n0].y);
	p1.Set(nodelist[seg.n1].x,nodelist[seg.n1].y);
	a0.Set(nodelist[arc.n0].x,nodelist[arc.n0].y);
	a1.Set(nodelist[arc.n1].x,nodelist[arc.n1].y);
	d=abs(a1-a0);			// distance between arc endpoints

	// figure out what the radius of the circle is...
	t=(a1-a0)/d;
	tta=arc.ArcLength*PI/180.;
	R=d/(2.*sin(tta/2.));
	c=a0 + (d/2. + I*sqrt(R*R-d*d/4.))*t; // center of the arc segment's circle...
	
	// figure out the distance between line and circle's center;
	d=abs(p1-p0);
	t=(p1-p0)/d;
	v=(c-p0)/t;
	if (fabs(Im(v))>R) return 0;
	l=sqrt( R*R - Im(v)*Im(v));	// Im(v) is distance between line and center...

	if ((l/R) < 1.e-05){		// case where line is very close to a tangent;
		p[i]=p0 + Re(v)*t;		// make it be a tangent.
		R=Re((p[i]-p0)/t);
		z=arg((p[i]-c)/(a0-c));
		if ((R>0) && (R<d) && (z>0.) && (z<tta)) i++;
		return i;
	}

	p[i]=p0 + (Re(v)+l)*t;		// first possible intersection;
	R=Re((p[i]-p0)/t);
	z=arg((p[i]-c)/(a0-c));
	if ((R>0) && (R<d) && (z>0.) && (z<tta)) i++;

	p[i]=p0 + (Re(v)-l)*t;		// second possible intersection
	R=Re((p[i]-p0)/t);
	z=arg((p[i]-c)/(a0-c));
	if ((R>0) && (R<d) && (z>0.) && (z<tta)) i++;

	// returns the number of valid intersections found;
	// intersections are returned in the array p[];
	return i;
	
}

int CFemmeDoc::GetArcArcIntersection(CArcSegment &arc0, CArcSegment &arc1, CComplex *p)
{
	CComplex a0,a1,t,c0,c1;
	double d,l,R0,R1,z0,z1,c,tta0,tta1;
	int i=0;

	a0.Set(nodelist[arc0.n0].x,nodelist[arc0.n0].y);
	a1.Set(nodelist[arc1.n0].x,nodelist[arc1.n0].y);

	GetCircle(arc1,c1,R1);
	GetCircle(arc0,c0,R0);

	d=abs(c1-c0);			// distance between centers
	
	if ((d>R0+R1) || (d<1.e-08)) return 0;
							// directly eliminate case where there can't
							// be any crossings....

	l=sqrt((R0+R1-d)*(d+R0-R1)*(d-R0+R1)*(d+R0+R1))/(2.*d);
	c=1.+(R0/d)*(R0/d)-(R1/d)*(R1/d);
	t=(c1-c0)/d;
	tta0=arc0.ArcLength*PI/180;
	tta1=arc1.ArcLength*PI/180;

	p[i]=c0 + (c*d/2.+ I*l)*t;		// first possible intersection;
	z0=arg((p[i]-c0)/(a0-c0));
	z1=arg((p[i]-c1)/(a1-c1));
	if ((z0>0.) && (z0<tta0) && (z1>0.) && (z1<tta1)) i++;
	
	if(fabs(d-R0+R1)/(R0+R1)< 1.e-05){
		p[i]=c0+ c*d*t/2.;
		return i;
	}
	
	p[i]=c0 + (c*d/2.-I*l)*t;		// second possible intersection
	z0=arg((p[i]-c0)/(a0-c0));
	z1=arg((p[i]-c1)/(a1-c1));
	if ((z0>0.) && (z0<tta0) && (z1>0.) && (z1<tta1)) i++;

	// returns the number of valid intersections found;
	// intersections are returned in the array p[];
	return i;
}

BOOL CFemmeDoc::AddArcSegment(CArcSegment &asegm, double tol)
{
	int i,j,k;
	CSegment segm;
	CArcSegment newarc;
	CComplex c,p[2];
	CArray< CComplex, CComplex&> newnodes;
	double R,d,dmin,t;

	newnodes.RemoveAll();

	// don't add if line is degenerate
	if (asegm.n0==asegm.n1) return FALSE;

	// don't add if the arc is already in the list;
	for(i=0;i<arclist.GetSize();i++){
		if ((arclist[i].n0==asegm.n0) && (arclist[i].n1==asegm.n1) &&
			(fabs(arclist[i].ArcLength-asegm.ArcLength)<1.e-02)) return FALSE;
		// arcs are ``the same'' if start and end points are the same, and if
		// the arc lengths are relatively close (but a lot farther than
		// machine precision...
	}

	// add proposed arc to the linelist
	asegm.IsSelected=FALSE;
	
	// check to see if there are intersections
	for(i=0;i<linelist.GetSize();i++)
	{
		j=GetLineArcIntersection(linelist[i],asegm,p);
		if(j>0) for(k=0;k<j;k++) newnodes.Add(p[k]);
	}
	for(i=0;i<arclist.GetSize();i++)
	{
		j=GetArcArcIntersection(asegm,arclist[i],p);
		if(j>0) for(k=0;k<j;k++) newnodes.Add(p[k]);
	}

	// add nodes at intersections
	if(tol==0)
	{
		if (nodelist.GetSize()<2) t=1.e-08;
		else{
			CComplex p0,p1;
			p0=nodelist[0].CC();
			p1=p0;
			for(i=1;i<nodelist.GetSize();i++)
			{
				if(nodelist[i].x<p0.re) p0.re=nodelist[i].x;
				if(nodelist[i].x>p1.re) p1.re=nodelist[i].x;
				if(nodelist[i].y<p0.im) p0.im=nodelist[i].y;
				if(nodelist[i].y>p1.im) p1.im=nodelist[i].y;
			}
			t=abs(p1-p0)*1.e-06;
		}
	}
	else t=tol;

	for(i=0;i<newnodes.GetSize();i++) 
		AddNode(newnodes[i].re,newnodes[i].im,t);

	// add proposed arc segment;
	arclist.Add(asegm);

	// check to see if proposed arc passes through other points;
    // if so, delete arc and create arcs that link intermediate points;
    // does this by recursive use of AddArcSegment;
	
    UnselectAll();
	GetCircle(asegm,c,R);
    if (tol==0) dmin=fabs(R*PI*asegm.ArcLength/180.)*1.e-05;
	else dmin=tol;
    k=arclist.GetSize()-1;
			
    for(i=0;i<nodelist.GetSize();i++)
    {
        if( (i!=asegm.n0) && (i!=asegm.n1) )
        {
			d=ShortestDistanceFromArc(CComplex(nodelist[i].x,nodelist[i].y),arclist[k]);
			
		//	MsgBox("d=%g dmin=%g",d,dmin);
			// what is the purpose of this test?
		//	if (abs(nodelist[i].CC()-nodelist[asegm.n0].CC())<2.*dmin) d=2.*dmin;
		//	if (abs(nodelist[i].CC()-nodelist[asegm.n1].CC())<2.*dmin) d=2.*dmin;