prob1big.cpp

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

#include<stdafx.h>
#include<stdio.h>
#include<math.h>
#include "fkn.h"
#include "fknDlg.h"
#include "complex.h"
#include "mesh.h"
#include "spars.h"
#include "FemmeDocCore.h"

BOOL CFemmeDocCore::Static2D(CBigLinProb &L)
{
	int i,j,k,w,s,sdi_iter,sdin;
	double Me[3][3],be[3];		// element matrices;
	double Mx[3][3],My[3][3],Mn[3][3];
	double l[3],p[3],q[3];		// element shape parameters;
	int n[3];					// numbers of nodes for a particular element;
	double a,K,r,t,x,y,B,B1,B2,mu,v[3],u[3],dv,res,lastres,Cduct;
	double *V_old,*V_sdi,*CircInt1,*CircInt2,*CircInt3;
	double c=PI*4.e-05;
	double units[]={2.54,0.1,1.,100.,0.00254,1.e-04};
	int Iter=0,pctr;
	BOOL LinearFlag=TRUE;
	BOOL SDIflag=FALSE;
	res=0;

	CElement *El;
	V_old=(double *) calloc(NumNodes,sizeof(double));

	// check to see if there are any SDI boundaries...
	// lineproplist[ meshele[i].e[j] ].BdryFormat==0
	for(i=0;i<NumLineProps;i++)
		if(lineproplist[i].BdryFormat==3) SDIflag=TRUE;

	if(SDIflag==TRUE){
		// there is an SDI boundary defined; check to see if it is in use
		SDIflag=FALSE;
		for(i=0;i<NumEls;i++)
			for(j=0;j<3;j++)
				if (lineproplist[meshele[i].e[j]].BdryFormat==3){
					SDIflag=TRUE;
					printf("Problem has SDI boundaries\n");
					i=NumEls;
					j=3;
				}
	}

	if (SDIflag==TRUE)
	{
		V_sdi=(double *) calloc(NumNodes,sizeof(double));
		sdin=2;
	}
	else sdin=1;

	// check to see if any circuits have been defined and process them;
	if (NumCircProps>0)
	{
		CircInt1=(double *)calloc(NumCircProps,sizeof(double));
		CircInt2=(double *)calloc(NumCircProps,sizeof(double));
		CircInt3=(double *)calloc(NumCircProps,sizeof(double));
		for(i=0;i<NumEls;i++){
			if(meshele[i].lbl>=0)
			if(labellist[meshele[i].lbl].InCircuit!=-1){
				El=&meshele[i];
				
				// get element area;
				for(k=0;k<3;k++) n[k]=El->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.;
			//	r=(meshnode[n[0]].x+meshnode[n[1]].x+meshnode[n[2]].x)/3.;

				// if coils are wound, they act like they have
				// a zero "bulk" conductivity...
				Cduct=blockproplist[El->blk].Cduct;
				if (abs(labellist[El->lbl].Turns)!=1) Cduct=0;
				
				// evaluate integrals;
				CircInt1[labellist[El->lbl].InCircuit]+=a;
				CircInt2[labellist[El->lbl].InCircuit]+=
					a*Cduct;
				CircInt3[labellist[El->lbl].InCircuit]+=
					blockproplist[El->blk].Jr*a*100.;
			}
		}

		for(i=0;i<NumCircProps;i++)
		{
			// Case 0 :: voltage gradient is applied to the region;
			// Case 1 :: flat current density is applied to the region;
			if (circproplist[i].CircType==0)
			{
				if(CircInt2[i]==0){
					circproplist[i].Case=1;
					if(CircInt1[i]==0.) circproplist[i].J=0.;
					else circproplist[i].J=0.01*(circproplist[i].Amps_re -
						CircInt3[i])/CircInt1[i];
				}
				else{
					circproplist[i].Case=0;
					circproplist[i].dV=-0.01*(circproplist[i].Amps_re - 
						CircInt3[i])/CircInt2[i];
				}
			}
			else{
				circproplist[i].Case=0;
				circproplist[i].dV=circproplist[i].dVolts_re;
			}
		}
	}





	// first, need to define permeability in each block.  In nonlinear
	// case, this is sort of a hassle.  Linear permeability is simply
	// copied from the associated block definition, but nonlinear
	// permeability must be updated from iteration to iteration...

	// build element matrices using the matrices derived in Allaire's book.

do{
  for(sdi_iter=0;sdi_iter<sdin;sdi_iter++)
  {

	TheView->SetDlgItemText(IDC_FRAME1,"Matrix Construction");
	TheView->m_prog1.SetPos(0);
	pctr=0;

	if(Iter>0) L.Wipe();

	for(i=0;i<NumEls;i++){
		
		// update ``building matrix'' progress bar...
		j=(i*20)/NumEls+1;
		if(j>pctr){
			j=pctr*5; if (j>100) j=100;
			TheView->m_prog1.SetPos(j);
			pctr++;
		}

		// zero out Me, be;
		for(j=0;j<3;j++){
			for(k=0;k<3;k++){
				Me[j][k]=0.;
				Mx[j][k]=0.;
				My[j][k]=0.;
				Mn[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
		El=&meshele[i];		
	
		for(k=0;k<3;k++) n[k]=El->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;
		for(j=0,k=1;j<3;k++,j++){
			if (k==3) k=0;
			l[j]=sqrt( pow(meshnode[n[k]].x-meshnode[n[j]].x,2.) +
					   pow(meshnode[n[k]].y-meshnode[n[j]].y,2.) );
		}
		a=(p[0]*q[1]-p[1]*q[0])/2.;
		r=(meshnode[n[0]].x+meshnode[n[1]].x+meshnode[n[2]].x)/3.;

		// x-contribution; only need to do main diagonal and above;
		K = (-1./(4.*a));
		for(j=0;j<3;j++)
			for(k=j;k<3;k++)
			{	
				Mx[j][k] += K*p[j]*p[k];
				if (j!=k) Mx[k][j]+=K*p[j]*p[k];
			}
		// y-contribution; only need to do main diagonal and above;
		K = (-1./(4.*a));
		for(j=0;j<3;j++)
			for(k=j;k<3;k++)
			{
				My[j][k] +=K*q[j]*q[k];
				if (j!=k) My[k][j]+=K*q[j]*q[k];
			}

		// contributions to Me, be from derivative boundary conditions;
		for(j=0;j<3;j++){
			if (El->e[j] >= 0)
			if (lineproplist[El->e[j]].BdryFormat==2)
			{
				// conversion factor is 10^(-4) (I think...)
				K=-0.0001*c*lineproplist[ El->e[j] ].c0*l[j]/6.;
				k=j+1; if(k==3) k=0;
				Me[j][j]+=K*2.;
				Me[k][k]+=K*2.;
				Me[j][k]+=K;
				Me[k][j]+=K;

				K=(lineproplist[ El->e[j] ].c1*l[j]/2.)*0.0001;
				be[j]+=K;
				be[k]+=K;
			}
		}
		
		// contribution to be from current density in the block
		for(j=0;j<3;j++){
			if(labellist[El->lbl].InCircuit>=0)
			{
				k=labellist[El->lbl].InCircuit;
				if(circproplist[k].Case==1) t=circproplist[k].J.Re();
				if(circproplist[k].Case==0)
					t=-circproplist[k].dV.Re()*blockproplist[El->blk].Cduct;
			}
			else t=0;
			K=-(blockproplist[El->blk].Jr+t)*a/3.;
			be[j]+=K;
		}

		// contribution to be from magnetization in the block;
		for(j=0;j<3;j++){
			t=labellist[El->lbl].MagDir;
			k=j+1; if(k==3) k=0;
			// need to scale so that everything is in proper units...
			// conversion is 0.0001
			K=0.0001*blockproplist[El->blk].H_c*(
			  cos(t*PI/180.)*(meshnode[n[k]].x-meshnode[n[j]].x) +
			  sin(t*PI/180.)*(meshnode[n[k]].y-meshnode[n[j]].y) )/2.;
			be[j]+=K;
			be[k]+=K;
		}	

//////// Nonlinear Part

		// update permeability for the element;
		if (Iter==0){ 
			k=meshele[i].blk;
		
			if (blockproplist[k].BHpoints != 0) LinearFlag=FALSE;

			if (blockproplist[k].LamType==0){
				t=blockproplist[k].LamFill;
				meshele[i].mu1=blockproplist[k].mu_x*t + (1.-t);
				meshele[i].mu2=blockproplist[k].mu_y*t + (1.-t);
			}
			if (blockproplist[k].LamType==1){
				t=blockproplist[k].LamFill;
				mu=blockproplist[k].mu_x;
				meshele[i].mu1=mu*t + (1.-t);
				meshele[i].mu2=mu/(t + mu*(1.-t));
			}
			if (blockproplist[k].LamType==2){
				t=blockproplist[k].LamFill;
				mu=blockproplist[k].mu_y;
				meshele[i].mu2=mu*t + (1.-t);
				meshele[i].mu1=mu/(t + mu*(1.-t));
			}
			if (blockproplist[k].LamType>2)
			{
				meshele[i].mu1=1;
				meshele[i].mu2=1;
			}
		}
		else{ 
			k=meshele[i].blk;
		
			if ((blockproplist[k].LamType==0) &&
			    (meshele[i].mu1==meshele[i].mu2)
				&&(blockproplist[k].BHpoints>0))
			{
				for(j=0,B1=0.,B2=0.;j<3;j++){
					B1+=L.V[n[j]]*q[j];
					B2+=L.V[n[j]]*p[j];
				}
				B=c*sqrt(B1*B1+B2*B2)/(0.02*a); 
				// correction for lengths in cm of 1/0.02
	
				// find out new mu from saturation curve;
				blockproplist[k].GetBHProps(B,mu,dv);
				mu=1./(muo*mu);
				meshele[i].mu1=mu;
				meshele[i].mu2=mu;
				for(j=0;j<3;j++){
					for(w=0,v[j]=0;w<3;w++)
						v[j]+=(Mx[j][w]+My[j][w])*L.V[n[w]];
				}
				K=-200.*c*c*c*dv/a;
				for(j=0;j<3;j++)
					for(w=0;w<3;w++)
						Mn[j][w]=K*v[j]*v[w];
			}

			if ((blockproplist[k].LamType==1) && (blockproplist[k].BHpoints>0)){
				t=blockproplist[k].LamFill;