/*!\file Solverx
 * \brief solver
 */

#include "./Solverx.h"
#include "../../shared/shared.h"
#include "../../include/include.h"

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
        
void DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum){

	/*output: */
	IS          isv=NULL;
	IS          isp=NULL;

	int         start,end;
	IssmDouble*     df_local=NULL;
	int         df_local_size;
	int         i;

	int*     pressure_indices=NULL;
	int*     velocity_indices=NULL;
	int      pressure_num=0;
	int      velocity_num=0;
	int      pressure_count=0;
	int      velocity_count=0;

	if(typeenum==StokesSolverEnum){

		/*Ok, recover doftypes vector values and indices: */
		VecGetOwnershipRange(df,&start,&end);
		VecGetLocalSize(df,&df_local_size);
		VecGetArray(df,&df_local);

		pressure_num=0;
		velocity_num=0;
		for(i=0;i<df_local_size;i++){
			if (df_local[i]==PressureEnum)pressure_num++;
			else velocity_num++;
		}

		/*Allocate indices: */
		if(pressure_num)pressure_indices=xNew<int>(pressure_num);
		if(velocity_num)velocity_indices=xNew<int>(velocity_num);

		pressure_count=0;
		velocity_count=0;
		for(i=0;i<df_local_size;i++){
			if (df_local[i]==PressureEnum){
				pressure_indices[pressure_count]=start+i;
				pressure_count++;
			}
			if (df_local[i]==VelocityEnum){
				velocity_indices[velocity_count]=start+i;
				velocity_count++;
			}
		}
		VecRestoreArray(df,&df_local);
		
		/*Create indices sets: */
		#if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
		ISCreateGeneral(PETSC_COMM_WORLD,pressure_num,pressure_indices,&isp);
		ISCreateGeneral(PETSC_COMM_WORLD,velocity_num,velocity_indices,&isv);
		#else
		ISCreateGeneral(PETSC_COMM_WORLD,pressure_num,pressure_indices,PETSC_COPY_VALUES,&isp);
		ISCreateGeneral(PETSC_COMM_WORLD,velocity_num,velocity_indices,PETSC_COPY_VALUES,&isv);
		#endif
	}

	/*Free ressources:*/
	xDelete<int>(pressure_indices);
	xDelete<int>(velocity_indices);

	/*Assign output pointers:*/
	*pisv=isv;
	*pisp=isp;
}
