/*!\file:  TransformStiffnessMatrixCoord.cpp
 * \brief transform stiffness matrix coordinate system
 */ 
#include "./elements.h"

void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum){

	int* cs_array=NULL;

	/*All nodes have the same Coordinate System*/
	cs_array=xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=cs_enum;

	/*Call core*/
	TransformStiffnessMatrixCoord(Ke,nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}

void TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array){

	int     i,j;
	int     numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case XYEnum:   numdofs+=2; break;
			case XYZPEnum: numdofs+=4; break;
			default: _error2_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current stiffness matrix*/
	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes,numnodes,cs_array);

	/*Transform matrix: R^T*Ke*R */
	TripleMultiply(transform,numdofs,numdofs,1,
				values,Ke->nrows,Ke->ncols,0,
				transform,numdofs,numdofs,0,
				&Ke->values[0],0);

	/*Free Matrix*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}
