#include "stdafx.h"
#include <afx.h>
#include <afxtempl.h>
#include "problem.h"
#include "femm.h"
#include "xyplot.h"
#include "femmviewDoc.h"
#include "femmviewView.h"
#include "MainFrm.h"
#include "maskprogress.h"
#include "lua.h"

#include "spars.h"

extern BOOL bLinehook;
extern CLuaConsoleDlg *LuaConsole;

// #define SIMPLE

#ifdef SIMPLE

BOOL CFemmviewDoc::MakeMask()
{
	// Good placeholder mask generator
	// This gives a valid mask that butts right up against the 
	// selected region. This mask yields the same results as the
	// Coulomb virtual work method.

	int i,j;

	for(i=0;i<meshnode.GetSize();i++) meshnode[i].msk=0;

	for(i=0;i<meshelem.GetSize();i++)
	{
		if(blocklist[meshelem[i].lbl].IsSelected)
		{
			for(j=0;j<3;j++) meshnode[meshelem[i].p[j]].msk=1;
		}
	}

	return TRUE;
}

#else

BOOL CFemmviewDoc::MakeMask()
{
	if(bHasMask) return TRUE;

	int i,j,k,d;
	CBigLinProb L;
	CMaskProgress dlg;
	double bsq,dbsq,v;
	double Me[3][3],be[3];		// element matrix;
	double p[3],q[3];			// element shape parameters;
	double a;					// element area;
	int n[3];					// numbers of nodes for a particular element;
	int *matflag;
	int *lblflag;
	
	int NumNodes=(int) meshnode.GetSize();
	int NumEls=(int) meshelem.GetSize();
	BOOL bOnAxis=FALSE;

	static int plus1mod3[3] = {1, 2, 0};
	static int minus1mod3[3] = {2, 0, 1};

	//Display progress dialog
	if (bLinehook==FALSE){
		dlg.Create(IDD_MASKPROGRESS);
		dlg.ShowWindow(SW_SHOW);
		L.bar=&dlg.m_mask_progress_status;
		dlg.m_mask_progress_status.SetPos(0);
		L.Pump();
	}

	// figure out bandwidth--helps speed somethings up;
	int bw=0;
	for(i=0;i<NumEls;i++)
	{
		for(j=0;j<3;j++)
		{
			k=j+1; if (k==3) k=0;
			d=abs(meshelem[i].p[j]-meshelem[i].p[k]);
			if (d>bw) bw=d;
		}
	}
	bw++;

	L.Create(NumNodes,bw);

	 // if the problem is axisymmetric, does the selection lie along r=0?
	if(ProblemType==1)
		for(i=0;i<NumEls;i++)       
			if(blocklist[meshelem[i].lbl].IsSelected)
			{
				for(j=0;j<3;j++)
					if(meshnode[meshelem[i].p[j]].x<1.e-6)
					{
						   bOnAxis=TRUE;
						   break;
					}
				if (bOnAxis) break;
			}
	
	// Sort through materials to see if they denote air;
	matflag=(int*)calloc(blockproplist.GetSize(),sizeof(int));
	lblflag=(int*)calloc(blocklist.GetSize(),sizeof(int));
	for(i=0;i<blockproplist.GetSize();i++)
	{
		// k==0 for air, k==1 for other than air
		k=0;
		if((blockproplist[i].mu_x!=1) || (blockproplist[i].mu_y!=1)) k=1;
		if(blockproplist[i].BHpoints!=0) k=1;
		if(blockproplist[i].LamType!=0) k=1;		
		if(blockproplist[i].H_c!=0) k=1;
		if((blockproplist[i].Jr!=0) || (blockproplist[i].Ji!=0)) k=1;
		if(blockproplist[i].Cduct!=0) k=1;
		if((blockproplist[i].Theta_hn!=0) || 
		   (blockproplist[i].Theta_hx!=0) || 
		   (blockproplist[i].Theta_hy!=0)) k=1;
		matflag[i]=k;
	}

	// Now, sort through the labels to see which ones correspond to air blocks.
	for(i=0;i<blocklist.GetSize();i++)
	{
		lblflag[i]=matflag[blocklist[i].BlockType];
		if(blocklist[i].InCircuit>=0) lblflag[i]=1;
	}

	// Determine which nodal values should be fixed
	// and what values they should be fixed at;
	for(i=0;i<NumNodes;i++) L.V[i]=-1;

	// Figure out which nodes are exterior edges and set them to zero;
	for(i=0;i<NumEls;i++)
	{
		for(j=0;j<3;j++)
		{
			if (meshelem[i].n[j]==1)
			{
				k=meshelem[i].p[plus1mod3[j]];
				if((!bOnAxis) || (IsKosher(k))) L.V[k]=0;
				k=meshelem[i].p[minus1mod3[j]];
				if((!bOnAxis) || (IsKosher(k))) L.V[k]=0;
			}
		}
	}

	// Set all nodes in a selected block equal to 1;
	for(i=0;i<NumEls;i++)
	{
		if(blocklist[meshelem[i].lbl].IsSelected)
		{
			for(j=0;j<3;j++){
				L.V[meshelem[i].p[j]]=1;
			}
		}
		else if(lblflag[meshelem[i].lbl]!=0)
		{
			for(j=0;j<3;j++) L.V[meshelem[i].p[j]]=0;
		}
	}

	// Any nodes that have point currents applied to them but are not in the
	// selected region should also be set to zero so that they don't mess up
	// the force calculation
	if(nodeproplist.GetSize()>0)
	{
		CComplex *p;
		int npts;
		
		p=(CComplex *)calloc(nodelist.GetSize(),sizeof(CComplex));
		for(i=0,npts=0;i<nodelist.GetSize();i++)
			if(nodelist[i].BoundaryMarker>=0)
			{
				p[npts]=nodelist[i].CC();
				npts++;
			}
		
		if(npts>0) 
			for(i=0;i<NumNodes;i++) 
				for(j=0;j<npts;j++)
					if(abs(p[j]-meshnode[i].CC())<1.e-8)
					{
						if (L.V[i]<0) L.V[i]=0.;
						npts--;
						if(npts>0){
							p[j]=p[npts];
							j=npts;
						}
						else{
							j=npts;
							i=NumNodes;
						}
					}
		free(p);
	}

	// build up element matrices;
	for(i=0;i<NumEls;i++){

		// zero out Me;
		for(j=0;j<3;j++)
		{
			for(k=0;k<3;k++) Me[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	
	
		for(k=0;k<3;k++) n[k] = meshelem[i].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.;

		// quick check for consistency--
		// if the block is not air and is not selected,
		// all of the nodes in the block better be defined
		// to be zero;  Otherwise, the region for force
		// integration has been selected in an invalid way;
		if ((!blocklist[meshelem[i].lbl].IsSelected) && (lblflag[meshelem[i].lbl]))
		{
			for(j=0,k=0;j<3;j++) if (L.V[n[j]]==0) k++;
			if(k<3){
				CString outmsg;
				
				if (bLinehook!=FALSE){
					outmsg+="*** The selected region is invalid. A valid selection\r\n";
					outmsg+="*** cannot abut a region which is not free space.\r\n";
					MessageBeep(MB_ICONEXCLAMATION);
					LuaConsole->ToOutput(outmsg);
				}
				else{
					outmsg ="The selected region is invalid. A valid selection\n";
					outmsg+="cannot abut a region which is not free space.";
					MsgBox(outmsg);
				}
				free(matflag);
				free(lblflag);
				return FALSE;
			}
		}

		switch(WeightingScheme)
		{	
			case 1:
				// all elements evenly weighted
				v=1;
				break;

			case 2:
				// weights each element with the sqrt of its own area;
				v=sqrt(a);
				break;

			case 3:
				// determine a weighting for the element 
				// based on an error measure;
				for(j=0,bsq=0,dbsq=0;j<3;j++)
				{
					dbsq+=Re((meshelem[i].B1-meshelem[i].b1[j])*
						 conj(meshelem[i].B1-meshelem[i].b1[j]) +
							 (meshelem[i].B2-meshelem[i].b2[j])*
						 conj(meshelem[i].B2-meshelem[i].b2[j]));
					bsq +=Re(meshelem[i].B1*conj(meshelem[i].B1) +
							 meshelem[i].B2*conj(meshelem[i].B2));
				}
				if(bsq!=0) v=dbsq/bsq;
				else(v=1);
				break;

			case 4:
				// determine a weighting for the element 
				// based on an error measure;
				for(k=0,dbsq=0,bsq=0;k<3;k++)
					for(j=0;j<NumList[n[k]];j++)
					{
						dbsq+=Re((meshelem[i].B1-meshelem[ConList[n[k]][j]].B1)*
							 conj(meshelem[i].B1-meshelem[ConList[n[k]][j]].B1) +
							 (meshelem[i].B2-meshelem[ConList[n[k]][j]].B2)*
							 conj(meshelem[i].B2-meshelem[ConList[n[k]][j]].B2));
						bsq +=Re(meshelem[i].B1*conj(meshelem[i].B1) +
							 meshelem[i].B2*conj(meshelem[i].B2));
					}
				if(bsq!=0) v=dbsq/bsq;
				else(v=1);
				break;

			default:
				// Each element weighted by its region's
				// mesh size specification;
				v=blocklist[meshelem[i].lbl].MaxArea;
				if (v<=0) v=sqrt(a); else v=sqrt(v);
				break;			
		}
	
		// build element matrix;
		for(j=0;j<3;j++)
			for(k=0;k<3;k++)
				Me[j][k]+=v*(p[j]*p[k]+q[j]*q[k])/a; 
		
		// process any prescribed nodal values;
		// doing it here saves a lot of time.
		for(j=0;j<3;j++)
		{
			if(L.V[n[j]]>=0)
			{
				for(k=0;k<3;k++)
				{
					if(j!=k){
						be[k]-=Me[k][j]*L.V[n[j]];
						Me[k][j]=0;
						Me[j][k]=0;
					}
				}
				be[j]=L.V[n[j]]*Me[j][j];
			}
		}

		// combine block matrices into global matrices;
		for (j=0;j<3;j++)
		{
			for (k=j;k<3;k++)
				if(Me[j][k]!=0)	L.Put(L.Get(n[j],n[k])-Me[j][k],n[j],n[k]);
			L.b[n[j]]-=be[j];
		}
	}

	// solve the problem;
	bLinehook=BuildMask;
	if (L.PCGSolve(FALSE)==FALSE) return FALSE;
	bLinehook=FALSE;
	for(i=0;i<NumNodes;i++) 
	{
		switch(WeightingScheme)
		{	
			case 0:
				if(L.V[i]>0.5) meshnode[i].msk=1;
				else meshnode[i].msk=0;
				break;

			default:
				meshnode[i].msk = L.V[i];
				break;
		}
	}
	free(matflag);
	free(lblflag);
	bHasMask=TRUE;

	return TRUE;
}

#endif

BOOL CFemmviewDoc::IsKosher(int k)
{
    // If:
    //    1) this is an axisymmetric problem;
    //    2) the selected geometry lies along the r=0 axis, and
    //    3) we have a node on the r=0 axis that we are trying to determine
    //     if we should set to zero.
    // This routine determines whether the node is at the extents of
    // the r=0 domain (or lies at a break in some sub-interval).
	//
	// Returns TRUE if it is OK to define the node as zero;

	if((ProblemType==0) || (meshnode[k].x>1.e-6)) return TRUE;

    int i,j,n;
    int score=0;

    for(i=0;i<NumList[k];i++)
    {
        for(j=0;j<3;j++)
        {
            n=meshelem[ConList[k][i]].p[j];
            if((n!=k) && (meshnode[n].x<1.e-6))
            {
                score++;
                if(score>1) return FALSE;
            }
        }
    }

    return TRUE;
}