prob2big.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"

// #define NEWTON

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

#ifndef NEWTON
	CComplex murel,muinc;
#endif

	deg45=1+I;
	w=Frequency*2.*PI;
		
	CComplex *CircInt1,*CircInt2,*CircInt3;



	// Can't handle LamType==1 or LamType==2 in AC problems.
	// Detect if this is being attempted.
	for(i=0;i<NumEls;i++)
	{
		if( (blockproplist[meshele[i].blk].LamType==1) ||
			(blockproplist[meshele[i].blk].LamType==2) ){
			AfxMessageBox("On-edge lamination not supported in AC analyses");
			return FALSE;
		}
	}

	// Go through and evaluate permeability for regions subject to prox effects
	for(i=0;i<NumBlockLabels;i++) GetFillFactor(i);
	
	V_old=(CComplex *) calloc(NumNodes+NumCircProps,sizeof(CComplex));

	// check to see if any circuits have been defined and process them;
	if (NumCircProps>0)
	{
		CircInt1=(CComplex *)calloc(NumCircProps,sizeof(CComplex));
		CircInt2=(CComplex *)calloc(NumCircProps,sizeof(CComplex));
		CircInt3=(CComplex *)calloc(NumCircProps,sizeof(CComplex));
		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;
				
				// total cross-section of circuit;
				CircInt1[labellist[El->lbl].InCircuit]+=a; 

				// integral of conductivity over the circuit;
				CircInt2[labellist[El->lbl].InCircuit]+=a*Cduct;

				// integral of applied J over current;
				CircInt3[labellist[El->lbl].InCircuit]+=
					(blockproplist[El->blk].Jr+I*blockproplist[El->blk].Ji)*a*100.;
			}
		}

		for(i=0;i<NumCircProps;i++)
		{
			// Case 0 :: a priori known voltage gradient is applied to the region;
			// Case 1 :: flat current density is applied to the region;
			// Case 2 :: voltage gradient applied to the region, but we gotta solve for it...

			if (circproplist[i].CircType==0) // specified current
			{
				if(CircInt2[i]==0){ //circuit composed of zero cond. materials
					circproplist[i].Case=1;
					if (CircInt1[i]==0.) circproplist[i].J=0.;
					else circproplist[i].J=0.01*(
						(circproplist[i].Amps_re+I*circproplist[i].Amps_im) -
						 CircInt3[i])/CircInt1[i];
				}
				else{
					circproplist[i].Case=2; // need to include an extra
										    // entry in matrix to solve for
					                        // voltage grad in the circuit
				}
			}
			else{
				// case where voltage gradient is specified a priori...
				circproplist[i].Case=0;
				circproplist[i].dV=circproplist[i].dVolts_re +
								   I*circproplist[i].dVolts_im;
			}
		}
	}



	// 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=(CComplex *) calloc(NumNodes+NumCircProps,sizeof(CComplex));
		sdin=2;
	}
	else sdin=1;

	// compute effective permeability for each block type;
	Mu=(CComplex **)calloc(NumBlockProps,sizeof(CComplex *));
	for(i=0;i<NumBlockProps;i++) Mu[i]=(CComplex *)calloc(2,sizeof(CComplex));

	for(k=0;k<NumBlockProps;k++){

		if (blockproplist[k].LamType==0){
			Mu[k][0]=blockproplist[k].mu_x*exp(-I*blockproplist[k].Theta_hx*DEG);
			Mu[k][1]=blockproplist[k].mu_y*exp(-I*blockproplist[k].Theta_hy*DEG);
			
			if(blockproplist[k].Lam_d!=0){
				if (blockproplist[k].Cduct != 0){
					halflag=exp(-I*blockproplist[k].Theta_hx*DEG/2.);
					ds=sqrt(2./(0.4*PI*w*blockproplist[k].Cduct*blockproplist[k].mu_x));
					K=halflag*deg45*blockproplist[k].Lam_d*0.001/(2.*ds);
					Mu[k][0]=((Mu[k][0]*tanh(K))/K)*blockproplist[k].LamFill +
							(1.- blockproplist[k].LamFill);

					halflag=exp(-I*blockproplist[k].Theta_hy*DEG/2.);
					ds=sqrt(2./(0.4*PI*w*blockproplist[k].Cduct*blockproplist[k].mu_y));
					K=halflag*deg45*blockproplist[k].Lam_d*0.001/(2.*ds);
					Mu[k][1]=((Mu[k][1]*tanh(K))/K)*blockproplist[k].LamFill + 
							(1. - blockproplist[k].LamFill);
				}
				else{
					Mu[k][0]=Mu[k][0]*blockproplist[k].LamFill +
							(1.- blockproplist[k].LamFill);
					Mu[k][1]=Mu[k][1]*blockproplist[k].LamFill + 
							(1. - blockproplist[k].LamFill);
				}
			}
		}
		else{
			Mu[k][0]=1;
			Mu[k][1]=1;
		}

	}

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();
  
	// build element matrices using the matrices derived in Allaire's book.
	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.;

		// x-contribution; 
		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; 
		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];
			}

		// contribution from eddy currents;	
		K=-I*a*w*blockproplist[meshele[i].blk].Cduct*c/12.;

		// in-plane laminated blocks appear to have no conductivity;
		// eddy currents are accounted for in these elements by their
		// frequency-dependent permeability.
		if((blockproplist[El->blk].LamType==0) &&
			(blockproplist[El->blk].Lam_d>0)) K=0;
		
		// if this element is part of a wound coil, 
		// it should have a zero "bulk" conductivity...
		if(abs(labellist[El->lbl].Turns)!=1) K=0;

		for(j=0;j<3;j++)
		{
			for(k=j;k<3;k++){
				Me[j][k]+=K;
				Me[k][j]+=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]+=2*K;
					Me[k][k]+=2*K;
					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;
				}
			
				if (lineproplist[El->e[j]].BdryFormat==1)
				{
					ds=sqrt(2./(0.4*PI*w*lineproplist[El->e[j]].Sig*
						lineproplist[El->e[j]].Mu));
					K=deg45/(-ds*lineproplist[El->e[j]].Mu*100.);
					K*=(l[j]/6.);
					k=j+1; if(k==3) k=0;
					Me[j][j]+=2*K;
					Me[k][k]+=2*K;
					Me[j][k]+=K;
					Me[k][j]+=K;
				}
			}
		}
		
		// contribution to be from current density in the block
		for(j=0;j<3;j++){
			Jv=0;
			if(labellist[El->lbl].InCircuit>=0)
			{
				k=labellist[El->lbl].InCircuit;
				if(circproplist[k].Case==1) Jv=circproplist[k].J;
				if(circproplist[k].Case==0)
					Jv=-circproplist[k].dV*blockproplist[El->blk].Cduct;
			}
			K=-(blockproplist[El->blk].Jr+I*blockproplist[El->blk].Ji+Jv)*a/3.;
			be[j]+=K;

			if(labellist[El->lbl].InCircuit>=0){
				k=labellist[El->lbl].InCircuit;
				if(circproplist[k].Case==2) L.b[NumNodes+k]+=K;
			}
		}

		// do Case 2 circuit stuff for element
		if(labellist[El->lbl].InCircuit>=0){
			k=labellist[El->lbl].InCircuit;
			if(circproplist[k].Case==2){
				K=-I*a*w*blockproplist[meshele[i].blk].Cduct*c;
				for(j=0;j<3;j++) L.Put(L.Get(n[j],NumNodes+k)+K/3.,n[j],NumNodes+k);
				L.Put(L.Get(NumNodes+k,NumNodes+k)+K,NumNodes+k,NumNodes+k);
			}
		}


///////////////////////////////////////////////////////////////
//
//	New Nonlinear stuff
//
///////////////////////////////////////////////////////////////

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