/*!\file:  TransformSolutionCoord.cpp
 * \brief transform solution vector coordinate system
 */ 
#include "./elements.h"

void TransformSolutionCoord(IssmDouble* solution,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*/
	TransformSolutionCoord(solution,nodes,numnodes,cs_array);

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

void TransformSolutionCoord(IssmDouble* solution,Node** nodes,int numnodes,int* cs_array){

	int     i;
	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: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current solution vector*/
	values=xNew<IssmDouble>(numdofs);
	for(i=0;i<numdofs;i++) values[i]=solution[i];

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

	/*Transform matrix: R*U */
	MatrixMultiply(transform,numdofs,numdofs,0,
				values,numdofs,1,0,
				&solution[0],0);

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