// hsolvdoc.cpp : implementation of the Chsolvdoc class
//

 
#include <stdafx.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include "hsolv.h"
#include "hsolvDlg.h"
#include "mesh.h"
#include "spars.h"
#include "hsolvdoc.h"

/////////////////////////////////////////////////////////////////////////////
// Chsolvdoc construction/destruction

Chsolvdoc::Chsolvdoc()
{
	TheView=NULL;
	Precision=NULL;
	LengthUnits=NULL;
	ProblemType=NULL;
	Coords=NULL;
	BandWidth=NULL;
	NumNodes=NULL;
	NumEls=NULL;
	NumBlockProps=NULL;
	NumPBCs=NULL;
	NumLineProps=NULL;
	NumPointProps=NULL;
	NumCircProps=NULL;
	NumBlockLabels=NULL;

	meshnode=NULL;
	meshele=NULL;
	Tprev=NULL;
	blockproplist=NULL;
	lineproplist=NULL;
	nodeproplist=NULL;
	circproplist=NULL;
	labellist=NULL;
	pbclist=NULL;
	PathName=NULL;

	extRo=extRi=extZo=NULL;
}

Chsolvdoc::~Chsolvdoc()
{
	// This space for rent.
}

void Chsolvdoc::CleanUp()
{
	if (meshnode!=NULL)		free(meshnode);
	if (meshele!=NULL)		free(meshele);
	if (blockproplist!=NULL) free(blockproplist);
	if (lineproplist!=NULL)	 free(lineproplist);
	if (nodeproplist!=NULL)	 free(nodeproplist);
	if (circproplist!=NULL)	 free(circproplist);
	if (labellist!=NULL)	 free(labellist);
	if (pbclist!=NULL)		 free(pbclist);
	if (Tprev!=NULL)		 free(Tprev);
}

/////////////////////////////////////////////////////////////////////////////
// Chsolvdoc commands

char* StripKey(char *c)
{
	char *d;
	int i,k;

	k=(int) strlen(c);

	for(i=0;i<k;i++){
		if (c[i] == '='){
			d=c+i+1;
			return d;
		}
	} 

	return c+k;
}

BOOL Chsolvdoc::OnOpenDocument() 
{
	FILE *fp;
	int j,k;
	char s[1024],q[1024];
	char *v;
	CPointProp	  PProp;
	CBoundaryProp BProp;
	CMaterialProp MProp;
	CCircuit	  CProp;
	CBlockLabel   blk;

	sprintf(s,"%s.feh",PathName);
	if ((fp=fopen(s,"rt"))==NULL){
		fprintf(stderr,"Couldn't read from specified .feh file\n");
		return FALSE;
	}

	// define some defaults
	Precision=1.e-08;
	Depth=-1;
	ProblemType=0;
	Coords=0;
	NumPointProps=0;
	NumLineProps=0;
	NumBlockProps=0;
	NumCircProps=0;

	// parse the file

	while (fgets(s,1024,fp)!=NULL)
	{
		sscanf(s,"%s",q);

		// Precision
		if( _strnicmp(q,"[precision]",11)==0){
			v=StripKey(s);
			sscanf(v,"%lf",&Precision);
			q[0]=NULL;
		}

		// Units of length used by the problem
		if( _strnicmp(q,"[lengthunits]",13)==0){
			v=StripKey(s);
			sscanf(v,"%s",q);
			if( _strnicmp(q,"inches",6)==0) LengthUnits=0;
			else if( _strnicmp(q,"millimeters",11)==0) LengthUnits=1;
			else if( _strnicmp(q,"centimeters",1)==0) LengthUnits=2;
			else if( _strnicmp(q,"mils",4)==0) LengthUnits=4;
			else if( _strnicmp(q,"microns",6)==0) LengthUnits=5;
			else if( _strnicmp(q,"meters",6)==0) LengthUnits=3;
			q[0]=NULL;
		}

		// Depth for 2D planar problems;
		if( _strnicmp(q,"[depth]",7)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&Depth);
		   q[0]=NULL;
		}	

		// Problem Type (planar or axisymmetric)
		if( _strnicmp(q,"[problemtype]",13)==0){
			v=StripKey(s);
			sscanf(v,"%s",q);
			if( _strnicmp(q,"planar",6)==0) ProblemType=0;
			if( _strnicmp(q,"axisymmetric",3)==0) ProblemType=1;
			q[0]=NULL;
		}

		// Coordinates (cartesian or polar)
		if( _strnicmp(q,"[coordinates]",13)==0){
			v=StripKey(s);
			sscanf(v,"%s",q);
			if ( _strnicmp(q,"cartesian",4)==0) Coords=0;
			if ( _strnicmp(q,"polar",5)==0) Coords=1;
			q[0]=NULL;
		}
				
		// properties for axisymmetric external region
		if( _strnicmp(q,"[extzo]",7)==0){
			v=StripKey(s);
			sscanf(v,"%lf",&extZo);
			q[0]=NULL;
		}

		if( _strnicmp(q,"[extro]",7)==0){
			v=StripKey(s);
			sscanf(v,"%lf",&extRo);
			q[0]=NULL;
		}

		if( _strnicmp(q,"[extri]",7)==0){
			v=StripKey(s);
			sscanf(v,"%lf",&extRi);
			q[0]=NULL;
		}

		// Point Properties
		if( _strnicmp(q,"[pointprops]",12)==0){	
			v=StripKey(s);
			sscanf(v,"%i",&k);
			if (k>0) nodeproplist=(CPointProp *)calloc(k,sizeof(CPointProp));
			q[0]=NULL;
		}

		if( _strnicmp(q,"<beginpoint>",11)==0){	
			PProp.V=0;
			PProp.qp=0;
			q[0]=NULL;
		}

		if( _strnicmp(q,"<tp>",4)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&PProp.V);
		   q[0]=NULL;
		}	

		if( _strnicmp(q,"<qp>",4)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&PProp.qp);
		   q[0]=NULL;
		}	

		if( _strnicmp(q,"<endpoint>",9)==0){
			nodeproplist[NumPointProps]=PProp;
			NumPointProps++;
			q[0]=NULL;
		}

		// Boundary Properties;
		if( _strnicmp(q,"[bdryprops]",11)==0){	
			v=StripKey(s);
			sscanf(v,"%i",&k);
			if (k>0) lineproplist=(CBoundaryProp *)calloc(k,sizeof(CBoundaryProp));
			q[0]=NULL;
		}
 		
		if( _strnicmp(q,"<beginbdry>",11)==0){	
			BProp.BdryFormat=0;
			BProp.Tset=0;
			BProp.Tinf=0;
			BProp.qs=0;
			BProp.beta=0;
			BProp.h=0;
			q[0]=NULL;
		}

		if( _strnicmp(q,"<bdrytype>",10)==0){
		   v=StripKey(s);
		   sscanf(v,"%i",&BProp.BdryFormat);
		   q[0]=NULL;
		}
				if( _strnicmp(q,"<Tset>",6)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&BProp.Tset);
		   q[0]=NULL;
		}	

		if( _strnicmp(q,"<qs>",4)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&BProp.qs);
		   q[0]=NULL;
		}	
		
		if( _strnicmp(q,"<beta>",6)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&BProp.beta);
		   q[0]=NULL;
		}	
		
		if( _strnicmp(q,"<h>",3)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&BProp.h);
		   q[0]=NULL;
		}	

		if( _strnicmp(q,"<Tinf>",6)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&BProp.Tinf);
		   q[0]=NULL;
		}
	
		if( _strnicmp(q,"<endbdry>",9)==0){
			lineproplist[NumLineProps]=BProp;
			NumLineProps++;
			q[0]=NULL;
		}

		// Block Properties;
		if( _strnicmp(q,"[blockprops]",12)==0){	
			v=StripKey(s);
			sscanf(v,"%i",&k);
			if (k>0) blockproplist=(CMaterialProp *)calloc(k,sizeof(CMaterialProp));
			q[0]=NULL;
		}

        // timestep
        if( _strnicmp(q,"[dt]",4)==0){
            v=StripKey(s);
            sscanf(v,"%lf",&dT);
            q[0]=NULL;
        }

        // Previous Solution File
        if( _strnicmp(q,"[prevsoln]",10)==0){
			int i;
            v=StripKey(s);

            // have to do this carefully to accept a filename with spaces
            k=(int) strlen(v);
            for(i=0;i<k;i++)
                if(v[i]=='\"'){
                    v=v+i+1;
                    i=k;
                }

            k=(int) strlen(v);
            if(k>0) for(i=k-1;i>=0;i--){
                if(v[i]=='\"'){
                    v[i]=0;
                    i=-1;
                }
            }

            PrevSoln=v;
            q[0]=NULL;
        }

		if( _strnicmp(q,"<beginblock>",12)==0){	
			MProp.Kx=1;
			MProp.Ky=1;			// permittivity, relative
			MProp.Kt=0;
			MProp.qv=0;			// charge density, C/m^3
			MProp.npts=0;
			q[0]=NULL;
		}

		if( _strnicmp(q,"<Kx>",6)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&MProp.Kx);
		   q[0]=NULL;
		}
		
		if( _strnicmp(q,"<Ky>",6)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&MProp.Ky);
		   q[0]=NULL;
		}	

		if( _strnicmp(q,"<Kt>",6)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&MProp.Kt);
		   MProp.Kt*=1.e6;
		   q[0]=NULL;
		}	

		if( _strnicmp(q,"<qv>",5)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&MProp.qv);
		   q[0]=NULL;
		}	

		if( _strnicmp(q,"<TKPoints>",10)==0){
			v=StripKey(s);
			sscanf(v,"%i",&MProp.npts);
			if (MProp.npts>0)
			{
				for(j=0;j<MProp.npts;j++){
					fgets(s,1024,fp);
					sscanf(s,"%lf	%lf",&MProp.Kn[j].re,&MProp.Kn[j].im);
				}
			}
			q[0]=NULL;
		}

		if( _strnicmp(q,"<endblock>",9)==0){
			blockproplist[NumBlockProps]=MProp;
			NumBlockProps++;
			q[0]=NULL;
		}

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

		if( _strnicmp(q,"<beginconductor>",16)==0){	
			CProp.V=0;
			CProp.q=0;
			CProp.CircType=0;
			q[0]=NULL;
		}
		
		if( _strnicmp(q,"<tc>",4)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&CProp.V);
		   q[0]=NULL;
		}

		if( _strnicmp(q,"<qc>",4)==0){
		   v=StripKey(s);
		   sscanf(v,"%lf",&CProp.q);
		   q[0]=NULL;
		}
		
		if( _strnicmp(q,"<conductortype>",15)==0){
		   v=StripKey(s);
		   sscanf(v,"%i",&CProp.CircType);
		   q[0]=NULL;
		}
		
		if( _strnicmp(q,"<endconductor>",14)==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	%i",&blk.x,&blk.y,&blk.BlockType,
					&blk.MaxArea,&blk.InGroup,&blk.IsExternal);
				blk.IsDefault  = blk.IsExternal & 2;
				blk.IsExternal = blk.IsExternal & 1;
				blk.BlockType--;
				labellist[i]=blk;
			}
			q[0]=NULL;
		}
	}	

	fclose(fp);
	
	return TRUE;
}


BOOL Chsolvdoc::LoadPrev()
{
	if (PrevSoln.GetLength()==0) return TRUE;

	FILE *fp;
	double x,y;
	int k;
	char s[1024],q[256];

	if ((fp=fopen(PrevSoln,"rt"))==NULL){
		MsgBox("Couldn't read from specified previous solution\n");
		return FALSE;
	}

	// parse the file
	k=0;
	while (fgets(s,1024,fp)!=NULL)
	{
		sscanf(s,"%s",q);
		if( _strnicmp(q,"[solution]",11)==0){
			k=1;
			break;
		}
	}

	// case where the solution is never found.
	if (k==0)
	{
		fclose(fp);
		MsgBox("Couldn't read from specified previous solution\n");
		return FALSE;
	}

	// read in the solution
	fgets(s,1024,fp);
	sscanf(s,"%i",&k);
	if(k!=NumNodes)
	{
		fclose(fp);
		MsgBox("Previous solution mesh doesn't match current mesh\n");
		return FALSE;
	}

	Tprev=(double *)calloc(NumNodes,sizeof(double));

	for(k=0;k<NumNodes;k++)
	{
		fgets(s,1024,fp);
		sscanf(s,"%lf	%lf	%lf",&x,&y,&Tprev[k]);
	}

	return TRUE;
}

BOOL Chsolvdoc::LoadMesh()
{
	int i,j,k,q,n0,n1,n;
	char infile[256];
	FILE *fp;
	char s[1024];
//	double c[]={25.4,1.,10.,1000.,0.0254,0.001};
    double c[]={0.0254,0.001,0.01,1,2.54e-5,1.e-6};


	//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;

	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",&n);

		if(n>1){
			// strip off point BC number;
			j=n & 0xffff;
			j=j-2; 
			if (j<0) j=-1;
			
			// strip off Conductor number
			n= (n - (n & 0xffff))/0x10000 - 1;
		}
		else{
			j=-1;
			n=-1;
		}
		node.bc=j;
		node.InConductor=n;

		// convert all lengths to internal working units of mm
		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;

	int defaultLabel;
	for(i=0,defaultLabel=-1;i<NumBlockLabels;i++)
		if (labellist[i].IsDefault) defaultLabel=i;

	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) elm.lbl=defaultLabel;
		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.";
			MsgBox(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;		

	// 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",&n);
		
		// BC number;
		if (n<0)
		{
			n=(-n);
			j = (n & 0xffff) - 2;
			if (j<0) j = -1;

			// Conductor number;
			n= (n - (n & 0xffff))/0x10000 - 1;
			if (n>=0)
			{
				meshnode[n0].InConductor=n;
				meshnode[n1].InConductor=n;
			}
		}
		else j=-1;
		
		if (j>=0)
		{
			// search through elements to find one containing the line;
			// set corresponding edge equal to the bc number
			for(q=0,n=FALSE;q<nmbr[n0];q++)
			{
				elm=meshele[mbr[n0][q]];

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

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

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

				meshele[mbr[n0][q]]=elm;

				//this is a little hack: line charge distributions should be applied
				//to at most one element;
				if((lineproplist[j].BdryFormat==2) && (n)) q=nmbr[n0];
			}
		}

	}
	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;
}

CComplex CMaterialProp::GetK(double t)
{
	int i,j;
	
	// Kx returned as real part;
	// Ky returned as imag part

	if (npts==0) return (Kx+I*Ky);
	if (npts==1) return (Im(Kn[0])*(1+I));
	if (t<=Re(Kn[0])) return (Im(Kn[0])*(1+I));
	if (t>=Re(Kn[npts-1])) return (Im(Kn[npts-1])*(1+I));

	for(i=0,j=1;j<npts;i++,j++)
	{
		if((t>=Re(Kn[i])) && (t<=Re(Kn[j])))
		{
			return (1+I)*(Im(Kn[i])+Im(Kn[j]-Kn[i])*Re(t-Kn[i])/Re(Kn[j]-Kn[i]));
		}
	}

	return (Kx+I*Ky);
}