//MODEL class definition
//
//   Usage:
//      md = new model()

function model () {
	//methods
		this.disp = function() { //{{{
			console.log(sprintf("class model echo: "));
			console.log(sprintf("//19s: //-22s -- //s","mesh"            ,"[1x1 " + typeof(this.mesh) + "]","mesh properties"));
			console.log(sprintf("//19s: //-22s -- //s","mask"            ,"[1x1 " + typeof(this.mask) + "]","defines grounded and floating elements"));
			console.log(sprintf("//19s: //-22s -- //s","geometry"        ,"[1x1 " + typeof(this.geometry) + "]","surface elevation, bedrock topography, ice thickness,..."));
			console.log(sprintf("//19s: //-22s -- //s","constants"       ,"[1x1 " + typeof(this.constants) + "]","physical constants"));
			console.log(sprintf("//19s: //-22s -- //s","smb"             ,"[1x1 " + typeof(this.smb) + "]","surface mass balance"));
			console.log(sprintf("//19s: //-22s -- //s","basalforcings"   ,"[1x1 " + typeof(this.basalforcings) + "]","bed forcings"));
			console.log(sprintf("//19s: //-22s -- //s","materials"       ,"[1x1 " + typeof(this.materials) + "]","material properties"));
			console.log(sprintf("//19s: //-22s -- //s","damage"          ,"[1x1 " + typeof(this.damage) + "]","parameters for damage evolution solution"));
			console.log(sprintf("//19s: //-22s -- //s","friction"        ,"[1x1 " + typeof(this.friction) + "]","basal friction/drag properties"));
			console.log(sprintf("//19s: //-22s -- //s","flowequation"    ,"[1x1 " + typeof(this.flowequation) + "]","flow equations"));
			console.log(sprintf("//19s: //-22s -- //s","timestepping"    ,"[1x1 " + typeof(this.timestepping) + "]","time stepping for trans models"));
			console.log(sprintf("//19s: //-22s -- //s","initialization"  ,"[1x1 " + typeof(this.initialization) + "]","initial guess/state"));
			console.log(sprintf("//19s: //-22s -- //s","rifts"           ,"[1x1 " + typeof(this.rifts) + "]","rifts properties"));
			console.log(sprintf("//19s: //-22s -- //s","slr"             ,"[1x1 " + typeof(this.slr) + "]","slr forcings"));
			console.log(sprintf("//19s: //-22s -- //s","debug"           ,"[1x1 " + typeof(this.debug) + "]","debugging tools (valgrind, gprof)"));
			console.log(sprintf("//19s: //-22s -- //s","verbose"         ,"[1x1 " + typeof(this.verbose) + "]","verbosity level in solve"));
			console.log(sprintf("//19s: //-22s -- //s","settings"        ,"[1x1 " + typeof(this.settings) + "]","settings properties"));
			console.log(sprintf("//19s: //-22s -- //s","toolkits"        ,"[1x1 " + typeof(this.toolkits) + "]","PETSc options for each solution"));
			console.log(sprintf("//19s: //-22s -- //s","cluster"         ,"[1x1 " + typeof(this.cluster) + "]","cluster parameters (number of cpus...)"));
			console.log(sprintf("//19s: //-22s -- //s","balancethickness","[1x1 " + typeof(this.balancethickness) + "]","parameters for balancethickness solution"));
			console.log(sprintf("//19s: //-22s -- //s","stressbalance"   ,"[1x1 " + typeof(this.stressbalance) + "]","parameters for stressbalance solution"));
			console.log(sprintf("//19s: //-22s -- //s","groundingline"   ,"[1x1 " + typeof(this.groundingline) + "]","parameters for groundingline solution"));
			console.log(sprintf("//19s: //-22s -- //s","hydrology"       ,"[1x1 " + typeof(this.hydrology) + "]","parameters for hydrology solution"));
			console.log(sprintf("//19s: //-22s -- //s","masstransport"   ,"[1x1 " + typeof(this.masstransport) + "]","parameters for masstransport solution"));
			console.log(sprintf("//19s: //-22s -- //s","thermal"         ,"[1x1 " + typeof(this.thermal) + "]","parameters for thermal solution"));
			console.log(sprintf("//19s: //-22s -- //s","steadystate"     ,"[1x1 " + typeof(this.steadystate) + "]","parameters for steadystate solution"));
			console.log(sprintf("//19s: //-22s -- //s","trans"       ,"[1x1 " + typeof(this.trans) + "]","parameters for trans solution"));
			console.log(sprintf("//19s: //-22s -- //s","levelset"        ,"[1x1 " + typeof(this.levelset) + "]","parameters for moving boundaries (level-set method)"));
			console.log(sprintf("//19s: //-22s -- //s","calving"         ,"[1x1 " + typeof(this.calving) + "]","parameters for calving"));
			console.log(sprintf("//19s: //-22s -- //s","gia"             ,"[1x1 " + typeof(this.gia) + "]","parameters for gia solution"));
			console.log(sprintf("//19s: //-22s -- //s","autodiff"        ,"[1x1 " + typeof(this.autodiff) + "]","automatic differentiation parameters"));
			console.log(sprintf("//19s: //-22s -- //s","flaim"           ,"[1x1 " + typeof(this.flaim) + "]","flaim parameters"));
			console.log(sprintf("//19s: //-22s -- //s","inversion"       ,"[1x1 " + typeof(this.inversion) + "]","parameters for inverse methods"));
			console.log(sprintf("//19s: //-22s -- //s","qmu"             ,"[1x1 " + typeof(this.qmu) + "]","dakota properties"));
			console.log(sprintf("//19s: //-22s -- //s","outputdefinition","[1x1 " + typeof(this.outputdefinition) + "]","output definition"));
			console.log(sprintf("//19s: //-22s -- //s","results"         ,"[1x1 " + typeof(this.results) + "]","model results"));
			console.log(sprintf("//19s: //-22s -- //s","radaroverlay"    ,"[1x1 " + typeof(this.radaroverlay) + "]","radar image for plot overlay"));
			console.log(sprintf("//19s: //-22s -- //s","miscellaneous"   ,"[1x1 " + typeof(this.miscellaneous) + "]","miscellaneous fields"));
		} //}}}
		this.setdefaultparameters = function () { // {{{

			//initialize subclasses
			this.mesh             = new mesh2d();
			this.mask             = new mask();
			this.constants        = new constants();
			this.geometry         = new geometry();
			this.initialization   = new initialization();
			this.smb              = new SMBforcing();
			this.basalforcings    = new basalforcings();
			this.friction         = new friction();
			this.rifts            = new rifts();
			this.slr              = new slr();
			this.timestepping     = new timestepping();
			this.groundingline    = new groundingline();
			this.materials        = new matice();
			this.damage           = new damage();
			this.flowequation     = new flowequation();
			this.debug            = new debug();
			this.verbose          = new verbose();
			this.settings         = new settings();
			this.toolkits         = new toolkits();
			this.cluster          = new local();
			this.balancethickness = new balancethickness();
			this.stressbalance    = new stressbalance();
			this.hydrology        = new hydrologyshreve();
			this.masstransport    = new masstransport();
			this.thermal          = new thermal();
			this.steadystate      = new steadystate();
			this.trans            = new trans();
			this.levelset		  = new levelset();
			this.calving          = new calving();
			this.gia              = new gia();
			this.autodiff         = new autodiff();
			this.flaim            = new flaim();
			this.inversion        = new inversion();
			this.qmu              = new qmu();
			this.radaroverlay     = new radaroverlay();
			this.results          = {};
			this.outputdefinition = new outputdefinition();
			this.miscellaneous    = new miscellaneous();
			this.priv             = new priv();
		} //}}}
		this.checkmessage = function(string){ //{{{
			console.log('model not consistent: ' + string);
			md.private.isconsistent=false;
		} //}}}
		this.fix = function(){ //{{{

			for (var field in this){

				//Some properties do not need to be fixed
				if (field == 'results' | field =='radaroverlay' | field == 'toolkits' | field =='cluster' | field == 'flaim' | field == 'priv') continue;

				//Check that current field is a class
				if(typeof this[field] == 'function'){
					continue;
				}

				//Fix current object
				this[field].fix(this);
			}

		} //}}}
        this.extrude = function(md, ...args) { //{{{
        //EXTRUDE - vertically extrude a 2d mesh
        //
        //   vertically extrude a 2d mesh and create corresponding 3d mesh.
        //   The vertical distribution can:
        //    - follow a polynomial law
        //    - follow two polynomial laws, one for the lower part and one for the upper part of the mesh
        //    - be discribed by a list of coefficients (between 0 and 1)
        //   
        //
        //   Usage:
        //      md=extrude(md,numlayers,extrusionexponent);
        //      md=extrude(md,numlayers,lowerexponent,upperexponent);
        //      md=extrude(md,listofcoefficients);
        //
        //   Example:
        //      md=extrude(md,15,1.3);
        //      md=extrude(md,15,1.3,1.2);
        //      md=extrude(md,[0 0.2 0.5 0.7 0.9 0.95 1]);
        //
        //   See also: MODELEXTRACT, COLLAPSE

        //some checks on list of arguments
            var argc = args.length;

            if ((argc > 4) | (argc < 2)) {
                throw "extrude error message";
            }

            //Extrude the mesh
            if (argc==2) { //list of coefficients
                clist=args[0];

                if (clist.some(function(x) { return x < 0; }) | clist.some(function(x) { return x > 1; }))
                    throw 'extrusioncoefficients must be between 0 and 1';
                clistcopy.push.apply(clistcopy, [0, 1]);

                extrusionlist=clistcopy.reverse().filter(function(e, i, arr) { //only keep unique elements
                    return arr.indexOf(e, i+1) === -1;
                }).sort();

                console.log("extrusionlist is " + extrusionlist) //@REMOVE

                numlayers=extrusionlist.length;
            } else if (argc==3) { //one polynomial law
                if (args[2]<=0) {
                    //help extrude;
                    throw 'extrusionexponent must be >=0';
                }
                numlayers=args[0];

                extrusionlist = [];

                for (var i = 0; i < numlayers; i++) {
                    extrusionlist.push(Math.pow(i/(numlayers-1), args[2]));
                };

            } else if (argc==4) { //two polynomial laws
                numlayers=args[0];
                lowerexp=args[1];
                upperexp=args[2];

                if (args[1]<=0 | args[2]<=0) {
                    throw 'lower and upper extrusionexponents must be >=0';
                }

                lowerextrusionlist = [];
                for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
                    lowerextrusionlist.push(Math.pow(i, lowerexp)/2)
                };

                upperextrusionlist = [];
                for (var i = 0; i <= 1; i += (2/(numlayers-1))) {
                    upperextrusionlist.push(Math.pow(i, upperexp)/2)
                };

                //lowerextrusionlist=[(0:2/(numlayers-1):1).^lowerexp]/2;
                //upperextrusionlist=[(0:2/(numlayers-1):1).^upperexp]/2;

                var temp = lowerextrusionlist;
                temp.push.apply(temp, [1]);

                var temp2 = upperextrusionlist.map(function(x) { return -x; });
                temp.push.apply(temp, temp2);

                clistcopy.push.apply(clistcopy, [0, 1]);
                //extrusionlist=sort(unique([lowerextrusionlist 1-upperextrusionlist]));
                extrusionlist=temp.sort();
            }

            if (numlayers<2) {
                error('number of layers should be at least 2');
            }

            if (md.mesh.domaintype() !== '3D') {
                error('Cannot extrude a 3d mesh (extrude cannot be called more than once)');
            }

            //Initialize with the 2d mesh
            mesh2d = md.mesh;
            md.mesh=mesh3dprisms();
            md.mesh.x                           = mesh2d.x;
            md.mesh.y                           = mesh2d.y;
            md.mesh.elements                    = mesh2d.elements;
            md.mesh.numberofelements            = mesh2d.numberofelements;
            md.mesh.numberofvertices            = mesh2d.numberofvertices;

            md.mesh.lat                         = mesh2d.lat;
            md.mesh.long                        = mesh2d.long;
            md.mesh.epsg                        = mesh2d.epsg;

            md.mesh.vertexonboundary            = mesh2d.vertexonboundary;
            md.mesh.vertexconnectivity          = mesh2d.vertexconnectivity;
            md.mesh.elementconnectivity         = mesh2d.elementconnectivity;
            md.mesh.average_vertex_connectivity = mesh2d.average_vertex_connectivity;

            md.mesh.extractedvertices           = mesh2d.extractedvertices;
            md.mesh.extractedelements           = mesh2d.extractedelements;

            x3d=[]; 
            y3d=[];
            z3d=[];  //the lower node is on the bed
            thickness3d=md.geometry.thickness; //thickness and bed for these nodes
                bed3d=md.geometry.base;

            //Create the new layers
            for (var i = 1; i <= numlayers; i++) {
                x3d=[x3d, md.mesh.x]; 
                y3d=[y3d, md.mesh.y];
                //nodes are distributed between bed and surface accordingly to the given exponent
                z3d=[z3d, bed3d+thickness3d*extrusionlist(i)]; 
            };

            number_nodes3d=size(x3d,1); //number of 3d nodes for the non extruded part of the mesh

                //Extrude elements 
                elements3d=[];
                for (var i = 0; i < numlayers; i++) {
                    elements3d.push(md.mesh.elements+(i-1)*md.mesh.numberofvertices, md.mesh.elements+i*md.mesh.numberofvertices); //Create the elements of the 3d mesh for the non extruded part
                };

                number_el3d=size(elements3d,1); //number of 3d nodes for the non extruded part of the mesh

                function fillArray(value, len) {
                    if (len == 0) return [];
                    var a = [value];
                    while (a.length * 2 <= len) a = a.concat(a);
                    if (a.length < len) a = a.concat(a.slice(0, len - a.length));
                    return a;
                }

                //Keep a trace of lower and upper nodes
                //lowervertex=NaN*ones(number_nodes3d,1);
                //uppervertex=NaN*ones(number_nodes3d,1);

                lowervertex = fillArray(null, number_nodes3d);
                uppervertex = fillArray(null, number_nodes3d);

                //lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices;
                //uppervertex(1:(numlayers-1)*md.mesh.numberofvertices)=md.mesh.numberofvertices+1:number_nodes3d;

                for (var i = md.mesh.numberofvertices+1, k = 1; i < lowervertex.length && k <= (numlayers-1)*md.mesh.numberofvertices; i++, k++) {
                    lowervertex[i] = k;
                };

                for (var i = 1, k = md.mesh.numberofvertices+1; i <= (numlayers-1)*md.mesh.numberofvertices && k <= number_nodes3d; i++, k++) {
                    uppervertex[i] = k;
                };

                md.mesh.lowervertex=lowervertex;
                md.mesh.uppervertex=uppervertex;

                //same for lower and upper elements
                lowerelements = fillArray(null, number_el3d);
                upperelements = fillArray(null, number_el3d);
                //lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements;
                //upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements;

                for (var i = md.mesh.numberofvertices+1, k = 1; i < lowerelements.length && k <= (numlayers-2)*md.mesh.numberofelements; i++, k++) {
                    lowerelements[i] = k;
                };

                for (var i = 1, k = md.mesh.numberofelements + 1; i <= (numlayers-2)*md.mesh.numberofelements && k <= (numlayers-1)*md.mesh.numberofelements; i++, k++) {
                    upperelements[i] = k;
                };

                md.mesh.lowerelements=lowerelements;
                md.mesh.upperelements=upperelements;

                //Save old mesh 
                md.mesh.x2d=md.mesh.x;
                md.mesh.y2d=md.mesh.y;
                md.mesh.elements2d=md.mesh.elements;
                md.mesh.numberofelements2d=md.mesh.numberofelements;
                md.mesh.numberofvertices2d=md.mesh.numberofvertices;

                //Build global 3d mesh 
                md.mesh.elements=elements3d;
                md.mesh.x=x3d;
                md.mesh.y=y3d;
                md.mesh.z=z3d;
                md.mesh.numberofelements=number_el3d;
                md.mesh.numberofvertices=number_nodes3d;
                md.mesh.numberoflayers=numlayers;

                //Ok, now deal with the other fields from the 2d mesh:

                    //bedinfo and surface info
                md.mesh.vertexonbase=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',1);
                md.mesh.vertexonsurface=project3d(md,'vector',ones(md.mesh.numberofvertices2d,1),'type','node','layer',md.mesh.numberoflayers);
                md.mesh.vertexonboundary=project3d(md,'vector',md.mesh.vertexonboundary,'type','node');

                //lat long
                md.mesh.lat=project3d(md,'vector',md.mesh.lat,'type','node');
                md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');

                md.geometry=extrude(md.geometry,md);
                md.friction  = extrude(md.friction,md);
                md.inversion = extrude(md.inversion,md);
                md.smb = extrude(md.smb,md);
                md.initialization = extrude(md.initialization,md);

                md.flowequation=md.flowequation.extrude(md);
                md.stressbalance=extrude(md.stressbalance,md);
                md.thermal=md.thermal.extrude(md);
                md.masstransport=md.masstransport.extrude(md);
                md.levelset=extrude(md.levelset,md);
                md.calving=extrude(md.calving,md);
                md.hydrology = extrude(md.hydrology,md);

                //connectivity
                //if ~isnan(md.mesh.elementconnectivity)

                if (md.mesh.elementconnectivity.every(function(e, i, arr) { return e !== null && typeof e !== 'undefined'; })) {
                    var temparr = [];
                    for (var i = 0; i < numlayers; i++) {
                        temparr.push(md.mesh.elementconnectivity);
                    };

                    //md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);

                    md.mesh.elementconnectivity = temparr;

                    //md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;

                    for (var i = 0; i < md.mesh.elementconnectivity.length; i++) {
                        if (md.mesh.elementconnectivity[i] == 0) {
                            md.mesh.elementconnectivity[i] = null;
                        }
                    };

                    for (var i = 2; i < numlayers; i++) {
                        for (var j = (i-1)*md.mesh.nuberofelements2d+1; i <= (i)*md.mesh.numberofelements2d; i++) {
                            md.mesh.elementconnectivity[i] += md.mesh.numberofelements2d;
                        };
                        //md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
                            //=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
                    };

                    md.mesh.elementconnectivity(find(isnan(md.mesh.elementconnectivity)))=0;

                    for (var i = 0; i < md.mesh.elementconnectivity.length; ++i) {
                        if (md.mesh.elementconnectivity[i] === null || typeof md.mesh.elementconectivity[i] === 'undefined') {
                            md.mesh.elementconnectivity[i] = 0;
                        }
                    }
                }

                md.materials=extrude(md.materials,md);
                md.damage=extrude(md.damage,md);
                md.mask=extrude(md.mask,md);
                md.qmu=extrude(md.qmu,md);
                md.basalforcings=extrude(md.basalforcings,md);

                //increase connectivity if less than 25:
                if (md.mesh.average_vertex_connectivity<=25)
                    md.mesh.average_vertex_connectivity=100;
                // }}}

        } //}}}
	//properties
	// {{{
		//Careful here: no other class should be used as default value this is a bug of matlab
		this.mesh             = 0;
		this.mask             = 0;

		this.geometry         = 0;
		this.constants        = 0;
		this.smb              = 0;
		this.basalforcings    = 0;
		this.materials        = 0;
		this.damage           = 0;
		this.friction         = 0;
		this.flowequation     = 0;
		this.timestepping     = 0;
		this.initialization   = 0;
		this.rifts            = 0;
		this.slr              = 0;

		this.debug            = 0;
		this.verbose          = 0;
		this.settings         = 0;
		this.toolkits         = 0;
		this.cluster          = 0;

		this.balancethickness = 0;
		this.stressbalance    = 0;
		this.groundingline    = 0;
		this.hydrology        = 0;
		this.masstransport    = 0;
		this.thermal          = 0;
		this.steadystate      = 0;
		this.trans            = 0;
		this.levelset         = 0;
		this.calving          = 0;
		this.gia              = 0;

		this.autodiff         = 0;
		this.flaim            = 0;
		this.inversion        = 0;
		this.qmu              = 0;

		this.results          = 0;
		this.outputdefinition = 0;
		this.radaroverlay     = 0;
		this.miscellaneous    = 0;
		this.priv             = 0;

		//set default values for fields
		this.setdefaultparameters();

		//}}}
}
