problem.cpp

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

				b[k]    += (Md*x[k] + Mo*x[k+1]);
				b[k+1]  += (Mo*x[k] + Md*x[k+1]);
			}

			// set boundary conditions
			m1[0] = 0;
			b[0]  = 0;
			b[n] += Hdata[i];

			// solve tridiagonal equation
			// solution ends up in b;
			for(k = 0; k < n; k++)
			{
				c = m1[k]/m0[k];
				m0[k+1] -= m1[k]*c;
				b[k+1]  -= b[k]*c;
			}
			b[n] = b[n]/m0[n];
			for(k = n-1; k >= 0; k--)
				b[k] = (b[k] - m1[k]*b[k + 1])/m0[k];
	
			iter++;

			lastres=res;
			res=abs(b[n]-x[n])/d;
			
			if (res<1.e-8) Converged=TRUE;

			// Do the same relaxation scheme as is implemented
			// in the solver to make sure that this effective
			// lamination permeability calculation converges
			if(iter>5)
			{
				if ((res>lastres) && (Relax>0.1)) Relax/=2.;
				else Relax+= 0.1 * (1. - Relax);
			}

			for(k=0;k<=n;k++) x[k]=Relax*b[k]+(1.0-Relax)*x[k];

		}while(Converged!=TRUE);


		mu = x[n]/(Hdata[i]*d);

		free(x );
		free(b );
		free(m0);
		free(m1);

		return mu;
}

CComplex CMaterialProp::GetdHdB(double B)
{
	double b,z,l;
	CComplex h;
	int i;

	b=fabs(B);
	
	if(BHpoints==0)	return CComplex(b/(mu_x*muo));

	if(b>Bdata[BHpoints-1])
		return slope[BHpoints-1];

	for(i=0;i<BHpoints-1;i++)
		if((b>=Bdata[i]) && (b<=Bdata[i+1])){
			l=(Bdata[i+1]-Bdata[i]);
			z=(b-Bdata[i])/l;
			h=6.*z*(z-1.)*Hdata[i]/l +
			  (1.-4.*z+3.*z*z)*slope[i] +
			  6.*z*(1.-z)*Hdata[i+1]/l +
			  z*(3.*z-2.)*slope[i+1];
			return h;
		}

	return CComplex(0);
}

double CMaterialProp::GetH(double x)
{
	return Re(GetH(CComplex(x)));
}

CComplex CMaterialProp::GetH(CComplex x)
{
	double b,z,z2,l;
	CComplex p,h;
	int i;

	b=abs(x);
	if((BHpoints==0) || (b==0))	return 0;
	p=x/b;

	if(b>Bdata[BHpoints-1])
		return p*(Hdata[BHpoints-1] + slope[BHpoints-1]*(b-Bdata[BHpoints-1]));

	for(i=0;i<BHpoints-1;i++)
		if((b>=Bdata[i]) && (b<=Bdata[i+1])){
			l=Bdata[i+1]-Bdata[i];
			z=(b-Bdata[i])/l;
			z2=z*z;
			h=(1.-3.*z2+2.*z2*z)*Hdata[i] +
			  z*(1.-2.*z+z2)*l*slope[i] +
			  z2*(3.-2.*z)*Hdata[i+1] +
			  z2*(z-1.)*l*slope[i+1];
			return p*h;
		}

	return 0;
}

// GetEnergy for the magnetostatic case
double CMaterialProp::GetEnergy(double x)
{
	double b,z,z2,l,nrg;
	double b0,b1,h0,h1,dh0,dh1;
	int i;

	b=fabs(x);
	nrg=0.;
	
	if(BHpoints==0)	return 0;

	for(i=0;i<BHpoints-1;i++){

		b0=Bdata[i];	h0=Re(Hdata[i]);
		b1=Bdata[i+1];	h1=Re(Hdata[i+1]);
		dh0=Re(slope[i]);
		dh1=Re(slope[i+1]);

		if((b>=b0) && (b<=b1)){
			l=b1-b0;
			z=(b-b0)/l;
			z2=z*z;

			nrg += (dh0*l*l*(6. + z*(-8. + 3.*z))*z2)/12. + 
				(h0*l*z*(2. + (-2. + z)*z2))/2. - 
				(h1*l*(-2. + z)*z2*z)/2. +			
				(dh1*l*l*(-4. + 3.*z)*z2*z)/12;

			return nrg;
		}
		else{
			// point isn't in this section, but add in the
			// energy required to pass through it.
			b0=Bdata[i];	h0=Re(Hdata[i]);
			b1=Bdata[i+1];	h1=Re(Hdata[i+1]);
			dh0=Re(slope[i]);
			dh1=Re(slope[i+1]);
			nrg += ((b0 - b1)*((b0 - b1)*(dh0 - dh1) -
				6.*(h0 + h1)))/12.;
		}
	}		

	// if we've gotten to this point, the point is off the scale,
	// so we have to extrapolate the rest of the way...

	h0=Re(Hdata[BHpoints-1]);
	dh0=Re(slope[BHpoints-1]);
	b0=Bdata[BHpoints-1];

	nrg += ((b - b0)*(b*dh0 - b0*dh0 + 2*h0))/2.;
	
	return nrg;
}

double CMaterialProp::GetCoEnergy(double b)
{
	return (fabs(b)*GetH(b) - GetEnergy(b));
}

double CMaterialProp::DoEnergy(double b1, double b2)
{
	// calls the raw routine to get point energy, 
	// but deals with the load of special cases that
	// arise because I insisted on trying to deal with
	// laminations on a continuum basis.

	double nrg,biron,bair;

	// easiest case: the material is linear!
	if (BHpoints==0)
	{
		double h1,h2;
		
		if(LamType==0){		// laminated in-plane
			h1=b1/((1.+LamFill*(mu_x-1.))*muo);
			h2=b2/((1.+LamFill*(mu_y-1.))*muo);
		}

		if(LamType==1){		// laminated parallel to x;
			h1=b1/((1.+LamFill*(mu_x-1.))*muo);
			h2=b1*(LamFill/(mu_y*muo) + (1. - LamFill)/muo);
		}

		if(LamType==2){		// laminated parallel to x;
			h2=b1/((1.+LamFill*(mu_y-1.))*muo);
			h1=b1*(LamFill/(mu_x*muo) + (1. - LamFill)/muo);
		}

		if(LamType>2){
			h1=b1/muo;
			h2=b2/muo;
		}

		return ((h1*b1+h2*b2)/2.);
	}

	// Rats! The material is nonlinear.  Now, we have to do
	// a bit of work to get the energy.
	if(LamType==0) nrg = GetEnergy(sqrt(b1*b1+b2*b2));

	if(LamType==1){
		biron=sqrt((b1/LamFill)*(b1/LamFill) + b2*b2);
		bair=b2;
		nrg = LamFill*GetEnergy(biron)+(1-LamFill)*bair*bair/(2.*muo);
	}

	if(LamType==2){
		biron=sqrt((b2/LamFill)*(b2/LamFill) + b1*b1);
		bair=b1;
		nrg = LamFill*GetEnergy(biron)+(1-LamFill)*bair*bair/(2.*muo);
	}

	return nrg;
}

double CMaterialProp::DoCoEnergy(double b1, double b2)
{
	double nrg,biron,bair;

	// easiest case: the material is linear!
	// in this case, energy and coenergy are the same!
	if (BHpoints==0) return DoEnergy(b1,b2);

	if(LamType==0) nrg = GetCoEnergy(sqrt(b1*b1+b2*b2));

	if(LamType==1){
		biron=sqrt((b1/LamFill)*(b1/LamFill) + b2*b2);
		bair=b2;
		nrg = LamFill*GetCoEnergy(biron)+(1-LamFill)*bair*bair/(2.*muo);
	}

	if(LamType==2){
		biron=sqrt((b2/LamFill)*(b2/LamFill) + b1*b1);
		bair=b1;
		nrg = LamFill*GetCoEnergy(biron)+(1-LamFill)*bair*bair/(2.*muo);
	}

	return nrg;
}


double CMaterialProp::DoEnergy(CComplex b1, CComplex b2)
{
	// This one is meant for the frequency!=0 case.
	// Fortunately, there's not so much effort in this case.
	CComplex mu1,mu2,h1,h2;
	
	GetMu(b1,b2,mu1,mu2);
	h1=b1/(mu1*muo);
	h2=b2/(mu2*muo);
	return (Re(h1*conj(b1)+h2*conj(b2))/4.);

}

double CMaterialProp::DoCoEnergy(CComplex b1, CComplex b2)
{
	return DoEnergy(b1,b2);
}

void CMaterialProp::GetMu(CComplex b1, CComplex b2,
							CComplex &mu1, CComplex &mu2)
{
	// gets the permeability, given a flux density
	// version for frequency!=0

	CComplex biron;

	// easiest case: the material is linear!
	if (BHpoints==0)
	{		
		mu1=mu_fdx;
		mu2=mu_fdy;
		return;
	}

	// Rats! The material is nonlinear. 
	else{
		CComplex muiron;
	
		if(LamType==0){
			biron=sqrt(b1*conj(b1)+b2*conj(b2));
			if(abs(biron)<1.e-08) mu1=slope[0]; //catch degenerate case
			else mu1=biron/GetH(biron);
			mu2=mu1;
		}

		if(LamType==1){
			biron=sqrt((b1/LamFill)*(b1/LamFill) + b2*b2);
			if(abs(biron)<1.e-08) muiron=slope[0];
			else muiron=biron/GetH(biron);
			mu1=muiron*LamFill;
			mu2=1./(LamFill/muiron + (1. - LamFill)/muo);
		
		}

		if(LamType==2){
			biron=sqrt((b2/LamFill)*(b2/LamFill) + b1*b1);
			if(abs(biron)<1.e-08) muiron=slope[0];
			else muiron=biron/GetH(biron);
			mu2=muiron*LamFill;
			mu1=1./(LamFill/muiron + (1. - LamFill)/muo);
		}
	}	

	// convert to relative permeability
	mu1/=muo;
	mu2/=muo;

	return;
}


void CMaterialProp::GetMu(double b1, double b2,
						double &mu1, double &mu2)
{
	// gets the permeability, given a flux density
	//

	double biron;
	
	mu1=mu2=muo;			// default

	// easiest case: the material is linear!
	if (BHpoints==0)
	{		
		if(LamType==0){		// laminated in-plane
			mu1=((1.+LamFill*(mu_x-1.))*muo);
			mu2=((1.+LamFill*(mu_y-1.))*muo);
		}

		if(LamType==1){		// laminated parallel to x;
			mu1=((1.+LamFill*(mu_x-1.))*muo);
			mu2=1./(LamFill/(mu_y*muo) + (1. - LamFill)/muo);
		}

		if(LamType==2){		// laminated parallel to x;
			mu2=((1.+LamFill*(mu_y-1.))*muo);
			mu1=1./(LamFill/(mu_x*muo) + (1. - LamFill)/muo);
		}
	}

	// Rats! The material is nonlinear.
	else{
		double muiron;
	
		if(LamType==0){
			biron=sqrt(b1*b1+b2*b2);
			if(biron<1.e-08) mu1=Re(slope[0]); //catch degenerate case
			else mu1=biron/GetH(biron);
			mu2=mu1;
		}

		if(LamType==1){
			biron=sqrt((b1/LamFill)*(b1/LamFill) + b2*b2);
			if(biron<1.e-08) muiron=Re(slope[0]);
			else muiron=biron/GetH(biron);
			mu1=muiron*LamFill;
			mu2=1./(LamFill/muiron + (1. - LamFill)/muo);	
		}

		if(LamType==2){
			biron=sqrt((b2/LamFill)*(b2/LamFill) + b1*b1);
			if(biron<1.e-08) muiron=Re(slope[0]);
			else muiron=biron/GetH(biron);
			mu2=muiron*LamFill;
			mu1=1./(LamFill/muiron + (1. - LamFill)/muo);
		}

	}
	
	// convert to relative permeability
	mu1/=muo;
	mu2/=muo;

	return;
}









CBoundaryProp::CBoundaryProp()
{
		BdryName="New Boundary";
		BdryFormat=0;				// type of boundary condition we are applying
									// 0 = constant value of A
									// 1 = Small skin depth eddy current BC
									// 2 = Mixed BC

		A0=0.; A1=0.;
		A2=0.; phi=0.;			// set value of A for BdryFormat=0;
	
		Mu=0.; Sig=0.;			// material properties necessary to apply
								// eddy current BC
		
		c0=0.; c1=0.;			// coefficients for mixed BC

}

CPointProp::CPointProp()
{
		PointName="New Point Property";
		Jr=0.; Ji=0.;					// applied point current, A 
		Ar=0.; Ai=0.;					// prescribed nodal value;
}

CCircuit::CCircuit()
{
		CircName="New Circuit";
		CircType=0;
		Amps=0.;
}