junk.cpp

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

#include "stdafx.h"
#include <afx.h>
#include <afxtempl.h>
#include "complex.h"
#include "problem.h"
#include "femmview.h"
#include "xyplot.h"
#include "MainFrm.h"
#include "femmviewDoc.h"
#include "femmviewView.h"
#include "ptloc.h"
#include "lua.h"

#include "spars.h"

// #define SIMPLE

#ifdef SIMPLE

BOOL CFemmviewDoc::MakeMask()
{
	// Good placeholder mask generator
	// This actually gives a valid mask that butts right up against the 
	// part of interest.  This should actually give the same result as the
	// virtual work formulation that you get by differentiating element
	// shape functions.

	int i,j;

	for(i=0;i<meshnode.GetSize();i++) meshnode[i].msk=0;

	for(i=0;i<meshelem.GetSize();i++)
	{
		if(blocklist[meshelem[i].lbl].IsSelected)
		{
			for(j=0;j<3;j++) meshnode[meshelem[i].p[j]].msk=1;
		}
	}

	return TRUE;
}

#else

BOOL CFemmviewDoc::MakeMask()
{
	int i,j,k,d;
	CBigLinProb L;
	double Me[3][3],be[3];		// element matrix;
	double p[3],q[3];			// element shape parameters;
	double a;					// element area;
	int n[3];					// numbers of nodes for a particular element;
	int *matflag;
	int *lblflag;
	
	int NumNodes=meshnode.GetSize();
	int NumEls=meshelem.GetSize();
	int bw;

	BOOL bOnAxis=FALSE;

	static int plus1mod3[3] = {1, 2, 0};
	static int minus1mod3[3] = {2, 0, 1};

	// figure out bandwidth--helps speed somethings up;
	for(i=0,bw=0;i<NumEls;i++)
	{
		for(j=0;j<3;j++)
		{
			k=j+1; if (k==3) k=0;
			d=abs(meshelem[i].p[j]-meshelem[i].p[k]);
			if (d>bw) bw=d;
		}
	}
	bw++;
	
	L.Create(NumNodes,bw);
	
	// Sort through materials to see if they denote air;
	matflag=(int*)calloc(blockproplist.GetSize(),sizeof(int));
	lblflag=(int*)calloc(blocklist.GetSize(),sizeof(int));
	for(i=0;i<blockproplist.GetSize();i++)
	{
		// k==0 for air, k==1 for other than air
		k=0;
		if((blockproplist[i].mu_x!=1) || (blockproplist[i].mu_y!=1)) k=1;
		if(blockproplist[i].BHpoints!=0) k=1;
		if(blockproplist[i].LamType!=0) k=1;		
		if(blockproplist[i].H_c!=0) k=1;
		if((blockproplist[i].Jr!=0) || (blockproplist[i].Ji!=0)) k=1;
		if(blockproplist[i].Cduct!=0) k=1;
		if((blockproplist[i].Theta_hn!=0) || 
		   (blockproplist[i].Theta_hx!=0) || 
		   (blockproplist[i].Theta_hy!=0)) k=1;
		matflag[i]=k;
	}

	// Now, sort through the labels to see which ones correspond to air blocks.
	for(i=0;i<blocklist.GetSize();i++)
	{
		lblflag[i]=matflag[blocklist[i].BlockType];
		if(blocklist[i].InCircuit>=0) lblflag[i]=1;
	}

	// Figure out which nodes are exterior edges and set them to zero;
	for(i=0;i<NumEls;i++)
	{
		for(j=0;j<3;j++)
		{
			if ((meshelem[i].n[j] & 0x80000000) != FALSE)
			{
				L.V[meshelem[i].p[plus1mod3[j]]] =0;
				L.V[meshelem[i].p[minus1mod3[j]]]=0;
			}
		}
	} 

	// Set all nodes in a selected block equal to 1;
	// Set all nodes in unselected, non-air blocks to 0;
	for(i=0;i<NumEls;i++)
	{
		if(blocklist[meshelem[i].lbl].IsSelected)
		{
			for(j=0;j<3;j++){
				L.V[meshelem[i].p[j]]=1;
			}
		}
		else if(lblflag[meshelem[i].lbl]!=0)
		{
			for(j=0;j<3;j++) L.V[meshelem[i].p[j]]=0;
		}
	}

	// build up element matrices;
	for(i=0;i<NumEls;i++){

		// zero out Me;
		for(j=0;j<3;j++){
			for(k=0;k<3;k++) Me[j][k]=0;
			be[j]=0;
		}

		// Determine shape parameters.
		// l == element side lengths;
		// p corresponds to the `b' parameter in Allaire
		// q corresponds to the `c' parameter in Allaire	
	
		for(k=0;k<3;k++) n[k] = meshelem[i].p[k];
		p[0]=meshnode[n[1]].y - meshnode[n[2]].y;
		p[1]=meshnode[n[2]].y - meshnode[n[0]].y;
		p[2]=meshnode[n[0]].y - meshnode[n[1]].y;	
		q[0]=meshnode[n[2]].x - meshnode[n[1]].x;
		q[1]=meshnode[n[0]].x - meshnode[n[2]].x;
		q[2]=meshnode[n[1]].x - meshnode[n[0]].x;
		
		a = (p[0]*q[1]-p[1]*q[0])/2.;

		// x-contribution; only need to do main diagonal and above;
		for(j=0;j<3;j++)
			for(k=j;k<3;k++)
				Me[j][k]+=(p[j]*p[k]/a);
	
		// y-contribution; only need to do main diagonal and above;
		for(j=0;j<3;j++)
			for(k=j;k<3;k++)
				Me[j][k]+=(q[j]*q[k]/a);

/*
		// process any fixed boundary conditions;
		for(j=0;j<3;j++)
		{
			if(L.V[n[j]]>=0)
			{
				for(k=0;k<3;k++)
				{
					if(j!=k){
						be[k]-=Me[k][j]*L.V[n[j]];
						Me[k][j]=0;
						Me[j][k]=0;
					}
				}
			}
		}
*/
		// combine block matrices into global matrices;
		for (j=0;j<3;j++)
		{
			for (k=j;k<3;k++)
				L.Put(L.Get(n[j],n[k])-Me[j][k],n[j],n[k]);
//			L.b[n[j]]-=be[j];
		}
	}

	// Apply fixed nodal values;
	for(i=0;i<NumNodes;i++) 
		if (L.V[i]>=0) L.SetValue(i,L.V[i]);

//	for(i=0;i<NumNodes;i++)  if (L.V[i]>=0) L.b[i]=L.V[i]/L.Get(i,i);

	// solve the problem;
	if (L.PCGSolve(FALSE)==FALSE) return FALSE;
	for(i=0;i<NumNodes;i++) meshnode[i].msk = L.V[i];

	free(matflag);
	free(lblflag);
	return TRUE;
}

#endif