femmedoccore.cpp

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

		}	

		if( _strnicmp(q,"<BHPoints>",10)==0){
			v=StripKey(s);
			sscanf(v,"%i",&MProp.BHpoints);
			if (MProp.BHpoints>0)
			{
				MProp.Hdata=(CComplex *)calloc(MProp.BHpoints,sizeof(CComplex));
				MProp.Bdata=(double *)  calloc(MProp.BHpoints,sizeof(double));
				for(j=0;j<MProp.BHpoints;j++){
					fgets(s,1024,fp);
					sscanf(s,"%lf	%lf",&MProp.Bdata[j],&MProp.Hdata[j].re);
					MProp.Hdata[j].im=0;
				}
			}
			q[0]=NULL;
		}

		if( _strnicmp(q,"<endblock>",9)==0){
			blockproplist[NumBlockProps]=MProp;
			blockproplist[NumBlockProps].GetSlopes(Frequency*2.*PI);
			NumBlockProps++;
			MProp.BHpoints=0;
			MProp.Bdata=NULL;
			MProp.Hdata=NULL;
			q[0]=NULL;
		}

		// Circuit Properties
		if( _strnicmp(q,"[circuitprops]",15)==0){	
			v=StripKey(s);
			sscanf(v,"%i",&k);
			if(k>0) circproplist=(CCircuit *)calloc(k,sizeof(CCircuit));
			q[0]=NULL;
		}

		if( _strnicmp(q,"<begincircuit>",14)==0){	
			CProp.dVolts_re=0.;
			CProp.dVolts_im=0.;
			CProp.Amps_re=0.;
			CProp.Amps_im=0.;
			CProp.CircType=0;
			q[0]=NULL;
		}
		
		if( _strnicmp(q,"<voltgradient_re>",17)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&CProp.dVolts_re);
		   q[0]=NULL;
		}

		if( _strnicmp(q,"<voltgradient_im>",17)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&CProp.dVolts_im);
		   q[0]=NULL;
		}

		if( _strnicmp(q,"<totalamps_re>",14)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&CProp.Amps_re);
		   q[0]=NULL;
		}

		if( _strnicmp(q,"<totalamps_im>",14)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&CProp.Amps_im);
		   q[0]=NULL;
		}
		
		if( _strnicmp(q,"<circuittype>",13)==0){
		   v=StripKey(s);
		   sscanf(v,"%i",&CProp.CircType);
		   q[0]=NULL;
		}
		
		if( _strnicmp(q,"<endcircuit>",12)==0){
			circproplist[NumCircProps]=CProp;
			NumCircProps++;
			q[0]=NULL;
		}


		// read in regional attributes
		if(_strnicmp(q,"[numblocklabels]",13)==0){
			int i;
			v=StripKey(s);
			sscanf(v,"%i",&k);
			if (k>0) labellist=(CBlockLabel *)calloc(k, sizeof(CBlockLabel));
			NumBlockLabels=k;
			for(i=0;i<k;i++)
			{
				fgets(s,1024,fp);
				sscanf(s,"%lf	%lf	%i	%lf	%i	%lf	%i	%i	%i",&blk.x,&blk.y,&blk.BlockType,&blk.MaxArea,
					&blk.InCircuit,&blk.MagDir,&blk.InGroup,&blk.Turns,&blk.IsExternal);
				blk.BlockType--;
				blk.InCircuit--;
				labellist[i]=blk;
			}
			q[0]=NULL;
		}
	}	

	// need to set these so that valid BH data doesn't get wiped
	// by the destructor of MProp
	MProp.BHpoints=0;
	MProp.Bdata=NULL;
	MProp.Hdata=NULL;

	fclose(fp);
	
	if (NumCircProps==0) return TRUE;

	// Process circuits for serial connections.
	// The program deals with serial "circuits" by making a separate
	// circuit property for each block in the serial circuit.  Then,
	// each of this larger number of circuits can be processed using
	// the previous approach which considered all circuits to be
	// parallel connected.

	// first, make enough space for all possible circuits;
	CCircuit *temp=(CCircuit *)calloc(NumCircProps+NumBlockLabels,sizeof(CCircuit));
	for(k=0;k<NumCircProps;k++){
		temp[k]=circproplist[k];
		temp[k].OrigCirc=-1;
	}
	free(circproplist);
	circproplist=temp;
	NumCircPropsOrig=NumCircProps;

	// now, go through the block label list and make a new "circuit" 
	// for every block label that is an element of a "serial" circuit.
	CCircuit ncirc;
	for(k=0;k<NumBlockLabels;k++)
	if(labellist[k].InCircuit>=0){
		ic=labellist[k].InCircuit;
		if(circproplist[ic].CircType==1)
		{
			ncirc=circproplist[ic];
			ncirc.OrigCirc=ic;
			ncirc.Amps_im*=labellist[k].Turns;
			ncirc.Amps_re*=labellist[k].Turns;
			circproplist[NumCircProps]=ncirc;
			labellist[k].InCircuit=NumCircProps;
			NumCircProps++;
		}
	}

	// now, all "circuits" look like parallel circuits, so 
	for(k=0;k<NumCircProps;k++) 
		if(circproplist[k].CircType==1) circproplist[k].CircType=0;

	return TRUE;
}

BOOL CFemmeDocCore::LoadMesh()
{
	int i,j,k,q,n0,n1;
	char infile[256];
	FILE *fp;
	char s[1024];
	double c[]={2.54,0.1,1.,100.,0.00254,1.e-04};


	//read meshnodes;
	sprintf(infile,"%s.node",PathName);
	if((fp=fopen(infile,"rt"))==NULL){
		return FALSE;
	}
	fgets(s,1024,fp);
	sscanf(s,"%i",&k); NumNodes=k;

#ifdef CRIPPLEWARE
	if (NumNodes>999){
		fclose(fp);
		AfxMessageBox("This demo version only allows meshes with less than 1000 nodes");
		return FALSE;
	}
#endif

	meshnode=(CNode *)calloc(k,sizeof(CNode));
	CNode node;
	for(i=0;i<k;i++){
		fscanf(fp,"%i",&j);
		fscanf(fp,"%lf",&node.x);
		fscanf(fp,"%lf",&node.y);
		fscanf(fp,"%i",&j);
		if(j>1) j=j-2; else j=-1;
		node.bc=j;

		// convert all lengths to centimeters (better conditioning this way...)
		node.x*=c[LengthUnits];
		node.y*=c[LengthUnits];
		
		meshnode[i]=node;
	}
	fclose(fp);

	//read in periodic boundary conditions;
	sprintf(infile,"%s.pbc",PathName);
	if((fp=fopen(infile,"rt"))==NULL){
		return FALSE;
	}
	fgets(s,1024,fp);
	sscanf(s,"%i",&k); NumPBCs=k;

	if (k!=0) pbclist=(CCommonPoint *)calloc(k,sizeof(CCommonPoint));
	CCommonPoint pbc;
	for(i=0;i<k;i++){
		fscanf(fp,"%i",&j);
		fscanf(fp,"%i",&pbc.x);
		fscanf(fp,"%i",&pbc.y);
		fscanf(fp,"%i",&pbc.t);
		pbclist[i]=pbc;
	}
	fclose(fp);

	// read in elements;
	sprintf(infile,"%s.ele",PathName);
	if((fp=fopen(infile,"rt"))==NULL){
		return FALSE;
	}
	fgets(s,1024,fp);
	sscanf(s,"%i",&k); NumEls=k;

	meshele=(CElement *)calloc(k,sizeof(CElement));
	CElement elm;
	for(i=0;i<k;i++){
		fscanf(fp,"%i",&j);
		fscanf(fp,"%i",&elm.p[0]);
		fscanf(fp,"%i",&elm.p[1]);
		fscanf(fp,"%i",&elm.p[2]);
		fscanf(fp,"%i",&elm.lbl);
		elm.lbl--;
		if(elm.lbl<0){
			CString msg;
			msg ="Material properties have not been defined for\n"; 
			msg+="all regions. Press the \"Run Mesh Generator\"\n";
			msg+="button to highlight the problem regions.";
			AfxMessageBox(msg);
			fclose(fp);
			sprintf(infile,"%s.ele",PathName);	DeleteFile(infile);
			sprintf(infile,"%s.node",PathName);	DeleteFile(infile);
			sprintf(infile,"%s.pbc",PathName);	DeleteFile(infile);
			sprintf(infile,"%s.poly",PathName);	DeleteFile(infile);
			sprintf(infile,"%s.edge",PathName);	DeleteFile(infile);
			exit(1);
		}
		// look up block type out of the list of block labels
		elm.blk=labellist[elm.lbl].BlockType;

		meshele[i]=elm;
	}
	fclose(fp);

	// initialize edge bc's and element permeabilities;
	for(i=0;i<NumEls;i++)
		for(j=0;j<3;j++)
		{
			meshele[i].e[j] = -1;		
			meshele[i].mu1  = -1.;
			meshele[i].mu2  = -1.;
		}

	// read in edges to which boundary conditions are applied;

		// first, do a little bookkeeping so that element
		// associated with a give edge can be identified fast
		int *nmbr;
		int **mbr;

		nmbr=(int *)calloc(NumNodes,sizeof(int));

		// Make a list of how many elements that tells how
		// many elements to which each node belongs.
		for(i=0;i<NumEls;i++)
			for(j=0;j<3;j++)
				nmbr[meshele[i].p[j]]++;

		// mete out some memory to build a list of the
		// connectivity...
		mbr=(int **)calloc(NumNodes,sizeof(int *));
		for(i=0;i<NumNodes;i++){
			mbr[i]=(int *)calloc(nmbr[i],sizeof(int));
			nmbr[i]=0;
		}

		// fill up the connectivity information;
		for(i=0;i<NumEls;i++)
			for(j=0;j<3;j++)
			{
				k=meshele[i].p[j];
				mbr[k][nmbr[k]]=i;
				nmbr[k]++;
			}		

	sprintf(infile,"%s.edge",PathName);
	if((fp=fopen(infile,"rt"))==NULL)
	{
		return FALSE;
	}
	fscanf(fp,"%i",&k);	// read in number of lines

	fscanf(fp,"%i",&j);	// read in boundarymarker flag;
	for(i=0;i<k;i++)
	{
		fscanf(fp,"%i",&j);
		fscanf(fp,"%i",&n0);
		fscanf(fp,"%i",&n1);
		fscanf(fp,"%i",&j);
		
		if(j<0){
			j = -(j+2);
			// search through elements to find one containing the line;
			// set corresponding edge equal to the bc number.
			for(q=0;q<nmbr[n0];q++)
			{
				elm=meshele[mbr[n0][q]];

				if ((elm.p[0] == n0) && (elm.p[1] == n1)) elm.e[0]=j;
				if ((elm.p[0] == n1) && (elm.p[1] == n0)) elm.e[0]=j;

				if ((elm.p[1] == n0) && (elm.p[2] == n1)) elm.e[1]=j;
				if ((elm.p[1] == n1) && (elm.p[2] == n0)) elm.e[1]=j;

				if ((elm.p[2] == n0) && (elm.p[0] == n1)) elm.e[2]=j;
				if ((elm.p[2] == n1) && (elm.p[0] == n0)) elm.e[2]=j;

				meshele[mbr[n0][q]]=elm;
			}
		}
		
	}
	fclose(fp);

	// free up the connectivity information
	free(nmbr);
	for(i=0;i<NumNodes;i++) free(mbr[i]);
	free(mbr);

	// clear out temporary files
	sprintf(infile,"%s.ele",PathName);	DeleteFile(infile);
	sprintf(infile,"%s.node",PathName);	DeleteFile(infile);
	sprintf(infile,"%s.pbc",PathName);	DeleteFile(infile);
	sprintf(infile,"%s.poly",PathName);	DeleteFile(infile);

	return TRUE;
}

void CFemmeDocCore::GetFillFactor(int lbl)
{
	// Get the fill factor associated with a stranded and
	// current-carrying region.  For AC problems, also compute
	// the apparent conductivity and permeability for use in
	// post-processing the voltage.
	
	CMaterialProp* bp= &blockproplist[labellist[lbl].BlockType];
	CBlockLabel* bl= &labellist[lbl];
	double atot,awire,r,FillFactor;
	int i,wiretype;
	CComplex ufd;
	double W=2.*PI*Frequency;

	if ((Frequency==0) || (blockproplist[labellist[lbl].BlockType].LamType<3))
	{
		bl->ProximityMu=0;
		return;
	}

	wiretype=bp->LamType-3;
	// wiretype = 0 for magnet wire
	// wiretype = 1 for stranded but non-litz wire
	// wiretype = 2 for litz wire
	// wiretype = 3 for rectangular wire
	r=bp->WireD*0.0005;
	
	for(i=0,atot=0;i<NumEls;i++)
		  if(meshele[i].lbl==lbl) atot+=ElmArea(i);

	awire=PI*r*r;
	if (wiretype==3) awire*=(4./PI); // if rectangular wire;
	awire*=((double) bp->NStrands);
	awire*=((double) bl->Turns);
	
	if (atot==0) return;
	FillFactor=fabs(awire/atot);

	double w,d,h,o,fill,dd;
	
	// if stranded but non-litz, use an effective wire radius that
	// gives the same cross-section as total stranded area
	if (wiretype==1) r*=sqrt((double) bp->NStrands);

	if (wiretype!=3){
		d=r*sqrt(3.);
		h=PI/sqrt(3.)*r;
		w=r*sqrt(PI/(2.*sqrt(3.)*FillFactor));
		dd=sqrt(3.)*w;
	}
	else{
		d=2*r;
		h=2*r;
		w=r/sqrt(FillFactor);
		dd=w;
	}
	o=bp->Cduct*(h/w)*5.e5; // conductivity in S/m
	fill=d/dd; //fill for purposes of equivalent foil analysis

	// At this point, sanity-check the fill factor;
	if (fill>1)
	{
		CString mymsg;
		mymsg.Format("Block label at (%g,%g) has a fill factor",bl->x,bl->y);
		mymsg +=     "greater than the theoretical maximum.  Couldn't solve the problem.";
		AfxMessageBox(mymsg);
		exit(5);
	}

	// effective permeability for the equivalent foil.  Note that this is
	// the same equation as effective permeability of a lamination...
	if (o!=0) ufd=muo*tanh(sqrt(I*W*o*muo)*d/2.)/(sqrt(I*W*o*muo)*d/2.);
	else ufd=0;
	
	// relative complex permeability
	bl->ProximityMu=(fill*ufd+(1.-fill)*muo)/muo;
}

double CFemmeDocCore::ElmArea(int i)
{
	// returns element cross-section area in meter^2
	int j,n[3];
	double b0,b1,c0,c1;

	for(j=0;j<3;j++) n[j]=meshele[i].p[j];

	b0=meshnode[n[1]].y - meshnode[n[2]].y;
	b1=meshnode[n[2]].y - meshnode[n[0]].y;
	c0=meshnode[n[2]].x - meshnode[n[1]].x;
	c1=meshnode[n[0]].x - meshnode[n[2]].x;
	return 0.0001*(b0*c1-b1*c0)/2.;

}