/*! \file Model.cpp
 * \brief  Model structure that mirrors the matlab workspace structure. Servers for the serial 
 * and parallel runs.
 */

#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../shared/shared.h"
#include "../io/io.h"
#include "../include/globals.h"
#include "../include/macros.h"

#include <string.h>
#include "stdio.h"

#include "./Model.h"


/*!--------------------------------------------------
	NewModel
  --------------------------------------------------*/

Model* NewModel(void) {
	/*! create a new Model object */
	Model* model=NULL;

	model=(Model*)xmalloc(sizeof(Model));

	/*!initialize all pointers to 0: */


	model->repository=NULL;
	model->meshtype=NULL;
	model->analysis_type=NULL;
	model->sub_analysis_type=NULL;
	model->numberofelements=0;
	model->numberofnodes=0;
	model->x=NULL; 
	model->y=NULL;
	model->z=NULL;
	model->elements=NULL;
	model->elements_type=NULL;
	model->numberofnodes2d=0;
	model->elements2d=NULL;
	model->deadgrids=NULL;
	model->numlayers=0;
	model->uppernodes=NULL;
	model->gridonhutter=NULL;
	
	model->vx_obs=NULL;
	model->vy_obs=NULL;
	model->vx=NULL;
	model->vy=NULL;
	model->vz=NULL;
	model->pressure=NULL;
	model->elementonbed=NULL;
	model->elementonsurface=NULL;
	model->gridonbed=NULL;
	model->gridonsurface=NULL;
	model->gridonstokes=NULL;
	model->borderstokes=NULL;
	model->thickness=NULL;
	model->surface=NULL;
	model->bed=NULL;
	model->elementoniceshelf=NULL;
	model->gridonicesheet=NULL;

	model->drag_type=0;
	model->drag=NULL;
	model->p=NULL;
	
	
	model->numberofsegs_diag=0;
	model->numberofsegs_diag_stokes=0;
	model->segmentonneumann_diag=NULL;
	model->segmentonneumann_diag_stokes=NULL;
	model-> neumannvalues_diag=NULL;
	model-> gridondirichlet_diag=NULL;
	model-> dirichletvalues_diag=NULL;
	//prognostic
	model-> segmentonneumann_prog=NULL;
	model-> neumannvalues_prog=NULL;
	model-> gridondirichlet_prog=NULL;
	model-> dirichletvalues_prog=NULL;
	//prognostic2
	model-> segmentonneumann_prog2=NULL;
	model-> neumannvalues_prog2=NULL;
	//thermal
	model-> gridondirichlet_thermal=NULL;
	model-> dirichletvalues_thermal=NULL;
	model-> geothermalflux=NULL;
	
	/*!materials: */
	model->rho_water=0;
	model->rho_ice=0;
		model->g=0;
	model->n=0;
	model->B=NULL;

	/*!control methods: */
	model->control_type=NULL;

	/*!solution parameters: */
	model->fit=NULL;
	model->meanvel=0;
	model->epsvel=0;
	model->artificial_diffusivity=0;
	model->nsteps=0;
	model->tolx=0;
	model->maxiter=NULL;
	model->mincontrolconstraint=0;
	model->maxcontrolconstraint=0;
	model->debug=0;
	model->plot=0;
	model->eps_rel=0;
	model->eps_abs=0;
	model->dt=0;
	model->ndt=0;
	model->penalty_offset=0;
	model->penalty_melting=0;
	model->penalty_lock=0;
	model->sparsity=0;
	model->connectivity=0;
	model->lowmem=0;
	model->optscal=NULL;
	model->yts=0;
	model->viscosity_overshoot=0;
	model->stokesreconditioning=0;
	model->waitonlock=0;

	/*!thermal parameters: */
	model->beta=0;
	model->meltingpoint=0;
	model->latentheat=0;
	model->heatcapacity=0;
	model->thermalconductivity=0;
	model->min_thermal_constraints=0;
	model->mixed_layer_capacity=0;
	model->thermal_exchange_velocity=0;

	
	model->numrifts=0;
	model->riftsnumpenaltypairs=NULL;
	model->riftspenaltypairs=NULL;
	model->riftsfill=NULL;
	model->riftsfriction=NULL;

	/*!penalties: */
	model->numpenalties=0;
	model->penalties=NULL;
	model->penaltypartitioning=NULL;

	/*!basal: */
	model->melting=NULL;
	model->accumulation=NULL;

	/*elements type: */
	model->ishutter=0;
	model->ismacayealpattyn=0;
	model->isstokes=0;

	return model;
}


/*!--------------------------------------------------
	DeleteModel
  --------------------------------------------------*/

void DeleteModel(Model** pmodel){

	/*!Recover structure: */
	Model* model = *pmodel;
	
	int i;

	/*!Two cases here: 
	 * - serial mode: matlab's memory manager will take care of delete model when returning from Imp. Do nothing here, so as not to confuse 
	 *                the memory manager.
     * - in parallel, anything the io layer does (FetchData) did needs to be erased explicitely in the model.
	 */

	#ifdef _PARALLEL_
	xfree((void**)&model->x);
	xfree((void**)&model->y);
	xfree((void**)&model->z);
	xfree((void**)&model->elements);
	xfree((void**)&model->elements_type);
	xfree((void**)&model->gridonhutter);
	if (strcmp(model->meshtype,"3d")==0){
		xfree((void**)&model->elements2d);
		xfree((void**)&model->deadgrids);
		xfree((void**)&model->uppernodes);
	}
	xfree((void**)&model->solverstring);
	xfree((void**)&model->elementonbed);
	xfree((void**)&model->elementonsurface);
	xfree((void**)&model->gridonbed);
	xfree((void**)&model->gridonsurface);
	xfree((void**)&model->gridonstokes);
	xfree((void**)&model->borderstokes);
	xfree((void**)&model->thickness);
	xfree((void**)&model->surface);
	xfree((void**)&model->bed);
	xfree((void**)&model->vx_obs);
	xfree((void**)&model->vy_obs);
	xfree((void**)&model->vx);
	xfree((void**)&model->vy);
	xfree((void**)&model->vz);
	xfree((void**)&model->pressure);
	xfree((void**)&model->drag);
	xfree((void**)&model->p);
	xfree((void**)&model->q);
	xfree((void**)&model->elementoniceshelf);
	xfree((void**)&model->gridonicesheet);
	xfree((void**)&model->segmentonneumann_diag);
	xfree((void**)&model->segmentonneumann_diag_stokes);
	xfree((void**)&model->neumannvalues_diag);
	xfree((void**)&model->gridondirichlet_diag);
	xfree((void**)&model->dirichletvalues_diag);
	xfree((void**)&model->segmentonneumann_prog);
	xfree((void**)&model->neumannvalues_prog);
	xfree((void**)&model->gridondirichlet_prog);
	xfree((void**)&model->dirichletvalues_prog);
	xfree((void**)&model->segmentonneumann_prog2);
	xfree((void**)&model->neumannvalues_prog2);
	xfree((void**)&model->gridondirichlet_thermal);
	xfree((void**)&model->dirichletvalues_thermal);
	xfree((void**)&model->geothermalflux);
	xfree((void**)&model->melting);
	xfree((void**)&model->accumulation);
	xfree((void**)&model->B);
	xfree((void**)&model->n);
	xfree((void**)&model->fit);
	xfree((void**)&model->optscal);
	xfree((void**)&model->maxiter);


	/*!Delete structure fields: */
	xfree((void**)&model->repository);
	xfree((void**)&model->meshtype);
	xfree((void**)&model->analysis_type);
	xfree((void**)&model->sub_analysis_type);
	
	if(model->numrifts){
		for(i=0;i<model->numrifts;i++){
			double* riftpenaltypairs=model->riftspenaltypairs[i];
			xfree((void**)&riftpenaltypairs);
		}
		xfree((void**)&model->riftspenaltypairs);
		xfree((void**)&model->riftsnumpenaltypairs);
		xfree((void**)&model->riftsfill);
		xfree((void**)&model->riftsfriction);
	}
	
	xfree((void**)&model->penalties);
	xfree((void**)&model->penaltypartitioning);
	
	xfree((void**)&model->control_type);
	
	xfree((void**)&model->epart);
	xfree((void**)&model->npart);
	xfree((void**)&model->my_grids);
	xfree((void**)&model->my_bordergrids);
	
	/*!Delete entire structure: */
	xfree((void**)pmodel);
	#endif
}

/*!--------------------------------------------------
	ModelInit
  --------------------------------------------------*/

#undef __FUNCT__ 
#define __FUNCT__ "ModelInit"

int	ModelInit(Model** pmodel,ConstDataHandle model_handle){

	int i,j;
	
	/*output: */
	Model* model=NULL;

	/*Allocate model: */
	model=NewModel();

	/*In ModelInit, we get all the data that is not difficult to get, and that is small: */
	ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL); 
	ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL); 

	ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL);
	/*!Get numberofelements and numberofnodes: */
	ModelFetchData((void**)&model->numberofnodes,NULL,NULL,model_handle,"numberofgrids","Integer",NULL);
	ModelFetchData((void**)&model->numberofelements,NULL,NULL,model_handle,"numberofelements","Integer",NULL);

	/*!In case we are running 3d, we are going to need the collapsed and non-collapsed 2d meshes, from which the 3d mesh was extruded: */
	if (strcmp(model->meshtype,"3d")==0){
	
		/*!Deal with 2d mesh: */
		ModelFetchData((void**)&model->numberofelements2d,NULL,NULL,model_handle,"numberofelements2d","Integer",NULL);
		ModelFetchData((void**)&model->numberofnodes2d,NULL,NULL,model_handle,"numberofgrids2d","Integer",NULL);
		ModelFetchData((void**)&model->numlayers,NULL,NULL,model_handle,"numlayers","Integer",NULL);
	}

	/*elements type: */
	ModelFetchData((void**)&model->ishutter,NULL,NULL,model_handle,"ishutter","Integer",NULL);
	ModelFetchData((void**)&model->ismacayealpattyn,NULL,NULL,model_handle,"ismacayealpattyn","Integer",NULL);
	ModelFetchData((void**)&model->isstokes,NULL,NULL,model_handle,"isstokes","Integer",NULL);

	/*!Get drag_type, drag and p,q: */
	ModelFetchData((void**)&model->drag_type,NULL,NULL,model_handle,"drag_type","Integer",NULL);

	/*!Get materials: */
	ModelFetchData((void**)&model->rho_water,NULL,NULL,model_handle,"rho_water","Scalar",NULL);
	ModelFetchData((void**)&model->rho_ice,NULL,NULL,model_handle,"rho_ice","Scalar",NULL);
	ModelFetchData((void**)&model->g,NULL,NULL,model_handle,"g","Scalar",NULL);

	/*Get control parameters: */
	ModelFetchData((void**)&model->control_type,NULL,NULL,model_handle,"control_type","String",NULL); 

	/*!Get solution parameters: */
	ModelFetchData((void**)&model->yts,NULL,NULL,model_handle,"yts","Scalar",NULL);
	ModelFetchData((void**)&model->meanvel,NULL,NULL,model_handle,"meanvel","Scalar",NULL);
	ModelFetchData((void**)&model->epsvel,NULL,NULL,model_handle,"epsvel","Scalar",NULL);
	ModelFetchData((void**)&model->debug,NULL,NULL,model_handle,"debug","Integer",NULL);
	ModelFetchData((void**)&model->plot,NULL,NULL,model_handle,"plot","Integer",NULL);
	ModelFetchData((void**)&model->artificial_diffusivity,NULL,NULL,model_handle,"artificial_diffusivity","Integer",NULL);
	ModelFetchData((void**)&model->nsteps,NULL,NULL,model_handle,"nsteps","Integer",NULL);
	ModelFetchData((void**)&model->tolx,NULL,NULL,model_handle,"tolx","Scalar",NULL);
	ModelFetchData((void**)&model->mincontrolconstraint,NULL,NULL,model_handle,"mincontrolconstraint","Scalar",NULL);
	ModelFetchData((void**)&model->maxcontrolconstraint,NULL,NULL,model_handle,"maxcontrolconstraint","Scalar",NULL);
	ModelFetchData((void**)&model->eps_rel,NULL,NULL,model_handle,"eps_rel","Scalar",NULL);
	ModelFetchData((void**)&model->eps_abs,NULL,NULL,model_handle,"eps_abs","Scalar",NULL);
	ModelFetchData((void**)&model->dt,NULL,NULL,model_handle,"dt","Scalar",NULL);
	ModelFetchData((void**)&model->ndt,NULL,NULL,model_handle,"ndt","Scalar",NULL);
	ModelFetchData((void**)&model->penalty_offset,NULL,NULL,model_handle,"penalty_offset","Scalar",NULL);
	ModelFetchData((void**)&model->penalty_melting,NULL,NULL,model_handle,"penalty_melting","Scalar",NULL);
	ModelFetchData((void**)&model->penalty_lock,NULL,NULL,model_handle,"penalty_lock","Integer",NULL);
	ModelFetchData((void**)&model->sparsity,NULL,NULL,model_handle,"sparsity","Scalar",NULL);
	ModelFetchData((void**)&model->connectivity,NULL,NULL,model_handle,"connectivity","Integer",NULL);
	ModelFetchData((void**)&model->lowmem,NULL,NULL,model_handle,"lowmem","Integer",NULL);
	ModelFetchData((void**)&model->solverstring,NULL,NULL,model_handle,"solverstring","String",NULL);
	ModelFetchData((void**)&model->viscosity_overshoot,NULL,NULL,model_handle,"viscosity_overshoot","Scalar",NULL);
	ModelFetchData((void**)&model->stokesreconditioning,NULL,NULL,model_handle,"stokesreconditioning","Scalar",NULL);
	ModelFetchData((void**)&model->waitonlock,NULL,NULL,model_handle,"waitonlock","Integer",NULL);

	/*!Get thermal parameters: */
	ModelFetchData((void**)&model->beta,NULL,NULL,model_handle,"beta","Scalar",NULL);
	ModelFetchData((void**)&model->meltingpoint,NULL,NULL,model_handle,"meltingpoint","Scalar",NULL);
	ModelFetchData((void**)&model->latentheat,NULL,NULL,model_handle,"latentheat","Scalar",NULL);
	ModelFetchData((void**)&model->heatcapacity,NULL,NULL,model_handle,"heatcapacity","Scalar",NULL);
	ModelFetchData((void**)&model->thermalconductivity,NULL,NULL,model_handle,"thermalconductivity","Scalar",NULL);
	ModelFetchData((void**)&model->min_thermal_constraints,NULL,NULL,model_handle,"min_thermal_constraints","Integer",NULL);
	ModelFetchData((void**)&model->mixed_layer_capacity,NULL,NULL,model_handle,"mixed_layer_capacity","Scalar",NULL);
	ModelFetchData((void**)&model->thermal_exchange_velocity,NULL,NULL,model_handle,"thermal_exchange_velocity","Scalar",NULL);

	/*rifts: */
	ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL);

	/*Assign output pointers: */
	*pmodel=model;

	return 1;
}

/*!--------------------------------------------------
	ModelEcho
  --------------------------------------------------*/
void ModelEcho(Model* model,int which_part,int rank) {

	//which_part  determines what gets echoed, otherwise, we'll get too much output.
	//1-> penalties

	int i,j;

	if(which_part==1 && my_rank==rank && (strcmp(model->meshtype,"3d")==0)){
		printf("Model penalties: \n");
		printf("   number of penalties: %i\n",model->numpenalties);
		printf("   grids: \n");

		for(i=0;i<model->numpenalties;i++){
			for(j=0;j<model->numlayers;j++){
				printf("%i ",(int)*(model->penalties+model->numlayers*i+j));
			}
			printf("\n");
		}
	}
	
	if(which_part==2 && my_rank==rank){
		printf("Model rifts: \n");
		printf("   number of rifts: %i\n",model->numrifts);
		for(i=0;i<model->numrifts;i++){
			double* penaltypairs=model->riftspenaltypairs[i];
			printf("   rift #%i\n",i);
			for (j=0;j<model->riftsnumpenaltypairs[i];j++){
				printf("      grids %g %g  elements %g %g normal [%g,%g] length %g\n",*(penaltypairs+7*j+0),*(penaltypairs+7*j+1),*(penaltypairs+7*j+2),*(penaltypairs+7*j+3),
						*(penaltypairs+7*j+4),*(penaltypairs+7*j+5),*(penaltypairs+7*j+6));
			}
			printf("		  friction %g fill %i\n",model->riftsfriction[i],model->riftsfill[i]);
		}
	}
	cleanup_and_return: 
	return;
}