Ignore:
Timestamp:
10/28/13 14:43:03 (11 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 16554

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c/classes/Elements/Penta.cpp

    r16137 r16560  
    3030        this->inputs            = NULL;
    3131        this->parameters        = NULL;
    32         this->results           = NULL;
    3332}
    3433/*}}}*/
    … …  
    3635Penta::~Penta(){
    3736        delete inputs;
    38         delete results;
    3937        this->parameters=NULL;
    4038}
    … …  
    4341Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels)
    4442        :PentaRef(nummodels)
    45         ,ElementHook(nummodels,index+1,6,iomodel) //index+1: material id, iomodel->numberofelements+1: matpar id
    46                                                                       { //i is the element index
    47 
    48         int i;
     43        ,ElementHook(nummodels,index+1,6,iomodel){
     44
    4945        int penta_elements_ids[2];
    5046
    5147        /*Checks in debugging mode*/
    52         /*{{{*/
    5348        _assert_(iomodel->Data(MeshUpperelementsEnum));
    5449        _assert_(iomodel->Data(MeshLowerelementsEnum));
    55         /*}}}*/
    5650
    5751        /*id: */
    … …  
    6963        this->parameters=NULL;
    7064
    71         /*intialize inputs and results: */
     65        /*intialize inputs: */
    7266        this->inputs=new Inputs();
    73         this->results=new Results();
    7467
    7568        /*initialize pointers:*/
    … …  
    112105                penta->inputs=new Inputs();
    113106        }
    114         if(this->results){
    115                 penta->results=(Results*)this->results->Copy();
    116         }
    117         else{
    118                 penta->results=new Results();
    119         }
    120107        /*point parameters: */
    121108        penta->parameters=this->parameters;
    … …  
    161148void Penta::BasalFrictionCreateInput(void){
    162149
    163         /*Constants*/
    164         const int  numdof=NUMVERTICES*NDOF1;
    165 
    166150        /*Intermediaries */
    167         int    count;
    168         IssmDouble basalfriction[NUMVERTICES]={0,0,0,0,0,0};
    169         IssmDouble alpha2,vx,vy;
    170         Friction*  friction=NULL;
    171         GaussPenta* gauss=NULL;
     151        int         count;
     152        IssmDouble  basalfriction[NUMVERTICES];
     153        IssmDouble  alpha2                       ,vx,vy;
     154        Friction   *friction                   = NULL;
     155        GaussPenta *gauss                      = NULL;
    172156
    173157        /* Basal friction can only be found at the base of an ice sheet: */
    … …  
    437421
    438422        /*retrieve parameters: */
    439         ElementMatrix* Ke=NULL;
    440423        int analysis_type;
    441424        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    464447                        break;
    465448                #endif
    466                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
     449                case L2ProjectionBaseAnalysisEnum:
    467450                        return CreateBasalMassMatrix();
    468451                        break;
    … …  
    501484                #endif
    502485                default:
    503                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
     486                        _error_("analysis " << EnumToStringx(analysis_type) << " not supported yet");
    504487        }
    505488
    … …  
    572555
    573556        /*retrieve parameters: */
    574         ElementMatrix* Ke=NULL;
    575557        ElementVector* De=NULL;
    576558        int analysis_type;
    577559        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    578560
    579         /*Checks in debugging {{{*/
     561        /*Checks in debugging*/
    580562        _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    581         /*}}}*/
    582563
    583564        switch(analysis_type){
    … …  
    692673                        break;
    693674                #endif
    694                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    695                         return CreatePVectorSlope();
     675                case L2ProjectionBaseAnalysisEnum:
     676                        return CreatePVectorL2ProjectionBase();
    696677                        break;
    697678                case MasstransportAnalysisEnum:
    … …  
    775756}
    776757/*}}}*/
    777 /*FUNCTION Penta::CreatePVectorSlope {{{*/
    778 ElementVector* Penta::CreatePVectorSlope(void){
     758/*FUNCTION Penta::CreatePVectorL2ProjectionBase {{{*/
     759ElementVector* Penta::CreatePVectorL2ProjectionBase(void){
    779760
    780761        if (!IsOnBed()) return NULL;
    … …  
    782763        /*Call Tria function*/
    783764        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    784         ElementVector* pe=tria->CreatePVectorSlope();
     765        ElementVector* pe=tria->CreatePVectorL2Projection();
    785766        delete tria->material; delete tria;
    786767
    … …  
    855836        _printf_("   inputs\n");
    856837        inputs->DeepEcho();
    857         _printf_("   results\n");
    858         results->DeepEcho();
    859 }
    860 /*}}}*/
    861 /*FUNCTION Penta::DeleteResults {{{*/
    862 void  Penta::DeleteResults(void){
    863 
    864         /*Delete and reinitialize results*/
    865         delete this->results;
    866         this->results=new Results();
    867 
    868838}
    869839/*}}}*/
    … …  
    929899                NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    930900        }
     901        NewTemperatureInput->Configure(this->parameters);
     902        NewPrecipitationInput->Configure(this->parameters);
    931903
    932904        this->inputs->AddInput(NewTemperatureInput);
    … …  
    14171389        /*Recover input*/
    14181390        Input* input=inputs->GetInput(enumtype);
    1419         if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
    1420         int    numnodes = this->NumberofNodes();
     1391        if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
    14211392
    14221393        /* Start looping on the number of vertices: */
    14231394        GaussPenta* gauss=new GaussPenta();
    1424         for (int iv=0;iv<this->NumberofNodes();iv++){
     1395        for(int iv=0;iv<this->NumberofNodes();iv++){
    14251396                gauss->GaussNode(this->element_type,iv);
    14261397                input->GetInputValue(&pvalue[iv],gauss);
    … …  
    14331404
    14341405        Input* input=inputs->GetInput(enumtype);
     1406        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
     1407
     1408        GaussPenta* gauss=new GaussPenta();
     1409        gauss->GaussVertex(this->GetNodeIndex(node));
     1410
     1411        input->GetInputValue(pvalue,gauss);
     1412        delete gauss;
     1413}
     1414/*}}}*/
     1415/*FUNCTION Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype) {{{*/
     1416void Penta::GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype){
     1417
     1418        Input* input=this->material->inputs->GetInput(enumtype);
    14351419        if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
    14361420
    … …  
    15291513                break;
    15301514        case StressbalanceVerticalAnalysisEnum:
    1531                 //GetSolutionFromInputsStressbalanceVert(solution);
    15321515                GetSolutionFromInputsOneDof(solution, VzEnum);
    15331516                break;
    … …  
    15351518        #ifdef _HAVE_THERMAL_
    15361519        case ThermalAnalysisEnum:
    1537                 //GetSolutionFromInputsThermal(solution);
    15381520                GetSolutionFromInputsOneDof(solution, TemperatureEnum);
    15391521                break;
    15401522        case EnthalpyAnalysisEnum:
    1541                 //GetSolutionFromInputsEnthalpy(solution);
    15421523                GetSolutionFromInputsOneDof(solution, EnthalpyEnum);
    15431524                break;
    … …  
    16571638}
    16581639/*}}}*/
    1659 /*FUNCTION Penta::GetVectorFromResults{{{*/
    1660 void  Penta::GetVectorFromResults(Vector<IssmDouble>* vector,int offset,int enum_in,int interp){
    1661 
    1662         /*Get result*/
    1663         ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(offset);
    1664         if(elementresult->InstanceEnum()!=enum_in){
    1665                 _error_("Results of offset "<<offset<<" is "<<EnumToStringx(elementresult->InstanceEnum())<<" when "<<EnumToStringx(enum_in)<<" was expected");
    1666         } 
    1667         if(interp==P1Enum){
    1668                 int vertexpidlist[NUMVERTICES];
    1669                 int connectivity[NUMVERTICES];
    1670                 this->GetVertexSidList(&vertexpidlist[0]);
    1671                 this->GetConnectivityList(&connectivity[0]);
    1672                 elementresult->GetVectorFromResults(vector,&vertexpidlist[0],&connectivity[0],NUMVERTICES);
    1673         }
    1674         else if(interp==P0Enum){
    1675                 elementresult->GetElementVectorFromResults(vector,sid);
    1676         }
    1677         else{
    1678                 _printf_("Interpolation " << EnumToStringx(interp) << " not supported\n");
    1679         }
    1680 }
    1681 /*}}}*/
    16821640/*FUNCTION Penta::GetZcoord {{{*/
    16831641IssmDouble Penta::GetZcoord(GaussPenta* gauss){
    … …  
    16991657        /*Computeportion of the element that is grounded*/
    17001658
    1701         int         normal_orientation;
     1659        int         normal_orientation=0;
    17021660        IssmDouble  s1,s2;
    17031661        IssmDouble  levelset[NUMVERTICES];
    … …  
    17061664        GetInputListOnVertices(&levelset[0],levelsetenum);
    17071665
    1708         if(levelset[0]*levelset[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
     1666        if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
    17091667                /*Portion of the segments*/
    17101668                s1=levelset[2]/(levelset[2]-levelset[1]);
    17111669                s2=levelset[2]/(levelset[2]-levelset[0]);
    17121670
    1713                 if(levelset[2]>0) normal_orientation=1;
     1671                if(levelset[2]>0.) normal_orientation=1;
    17141672                /*New point 1*/
    17151673                xyz_zero[3*normal_orientation+0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
    … …  
    17321690                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5][2]+s2*(xyz_list[3][2]-xyz_list[5][2]);
    17331691        }
    1734         else if(levelset[1]*levelset[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
     1692        else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
    17351693                /*Portion of the segments*/
    17361694                s1=levelset[0]/(levelset[0]-levelset[2]);
    17371695                s2=levelset[0]/(levelset[0]-levelset[1]);
    17381696
    1739                 if(levelset[0]>0) normal_orientation=1;
     1697                if(levelset[0]>0.) normal_orientation=1;
    17401698                /*New point 1*/
    17411699                xyz_zero[3*normal_orientation+0]=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
    … …  
    17581716                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3][2]+s2*(xyz_list[4][2]-xyz_list[3][2]);
    17591717        }
    1760         else if(levelset[0]*levelset[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
     1718        else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
    17611719                /*Portion of the segments*/
    17621720                s1=levelset[1]/(levelset[1]-levelset[0]);
    17631721                s2=levelset[1]/(levelset[1]-levelset[2]);
    17641722
    1765                 if(levelset[1]>0) normal_orientation=1;
     1723                if(levelset[1]>0.) normal_orientation=1;
    17661724                /*New point 0*/
    17671725                xyz_zero[3*normal_orientation+0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
    … …  
    17841742                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4][2]+s2*(xyz_list[5][2]-xyz_list[4][2]);
    17851743        }
    1786         else if(levelset[0]==0 && levelset[1]==0){ //front is on point 0 and 1
     1744        else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
    17871745                xyz_zero[3*0+0]=xyz_list[0][0];
    17881746                xyz_zero[3*0+1]=xyz_list[0][1];
    … …  
    18041762                xyz_zero[3*3+2]=xyz_list[3][2];
    18051763        }
    1806         else if(levelset[0]==0 && levelset[2]==0){ //front is on point 0 and 1
     1764        else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
    18071765                xyz_zero[3*0+0]=xyz_list[2][0];
    18081766                xyz_zero[3*0+1]=xyz_list[2][1];
    … …  
    18241782                xyz_zero[3*3+2]=xyz_list[5][2];
    18251783        }
    1826         else if(levelset[1]==0 && levelset[2]==0){ //front is on point 0 and 1
     1784        else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
    18271785                xyz_zero[3*0+0]=xyz_list[1][0];
    18281786                xyz_zero[3*0+1]=xyz_list[1][1];
    … …  
    18591817}
    18601818/*}}}*/
    1861 /*FUNCTION Penta::InputArtificialNoise{{{*/
    1862 void  Penta::InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max){
    1863 
    1864         Input* input=NULL;
    1865 
    1866         /*Make a copy of the original input: */
    1867         input=(Input*)this->inputs->GetInput(enum_type);
    1868         if(!input)_error_("could not find old input with enum: " << EnumToStringx(enum_type));
    1869 
    1870         /*ArtificialNoise: */
    1871         input->ArtificialNoise(min,max);
    1872 }
    1873 /*}}}*/
    1874 /*FUNCTION Penta::InputCreate(IssmDouble scalar,int enum,int code);{{{*/
    1875 void Penta::InputCreate(IssmDouble scalar,int name,int code){
    1876 
    1877         /*Check that name is an element input*/
    1878         if (!IsInput(name)) return;
    1879 
    1880         if ((code==5) || (code==1)){ //boolean
    1881                 this->inputs->AddInput(new BoolInput(name,reCast<bool,IssmDouble>(scalar)));
    1882         }
    1883         else if ((code==6) || (code==2)){ //integer
    1884                 this->inputs->AddInput(new IntInput(name,reCast<int,IssmDouble>(scalar)));
    1885         }
    1886         else if ((code==7) || (code==3)){ //IssmDouble
    1887                 this->inputs->AddInput(new DoubleInput(name,scalar));
    1888         }
    1889         else _error_("could not recognize nature of vector from code " << code);
    1890 
    1891 }
    1892 /*}}}*/
    1893 /*FUNCTION Penta::InputCreate(IssmDouble* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/
    1894 void Penta::InputCreate(IssmDouble* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements
     1819/*FUNCTION Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{*/
     1820void Penta::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){
    18951821
    18961822        /*Intermediaries*/
    … …  
    19081834                for(i=0;i<6;i++){
    19091835                        _assert_(iomodel->elements);
    1910                         penta_vertex_ids[i]=iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab
     1836                        penta_vertex_ids[i]=iomodel->elements[6*this->sid+i]; //ids for vertices are in the elements array from Matlab
    19111837                }
    19121838
    … …  
    19471873
    19481874                        if (code==5){ //boolean
    1949                                 this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool,IssmDouble>(vector[index])));
     1875                                this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool,IssmDouble>(vector[this->sid])));
    19501876                        }
    19511877                        else if (code==6){ //integer
    1952                                 this->inputs->AddInput(new IntInput(vector_enum,reCast<int,IssmDouble>(vector[index])));
     1878                                this->inputs->AddInput(new IntInput(vector_enum,reCast<int,IssmDouble>(vector[this->sid])));
    19531879                        }
    19541880                        else if (code==7){ //IssmDouble
    1955                                 this->inputs->AddInput(new DoubleInput(vector_enum,vector[index]));
     1881                                this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->sid]));
    19561882                        }
    19571883                        else _error_("could not recognize nature of vector from code " << code);
    … …  
    21472073}
    21482074/*}}}*/
    2149 /*FUNCTION Penta::InputToResult{{{*/
    2150 void  Penta::InputToResult(int enum_type,int step,IssmDouble time){
    2151 
    2152         bool   found = false;
    2153         Input *input = NULL;
    2154 
    2155         /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
    2156         if (enum_type==MaterialsRheologyBbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyBEnum);
    2157         else if (enum_type==MaterialsRheologyZbarEnum) input=this->material->inputs->GetInput(MaterialsRheologyZEnum);
    2158         else input=this->inputs->GetInput(enum_type);
    2159         //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in penta->inputs"); why error out? if the requested input does not exist, we should still
    2160         //try and output whatever we can instead of just failing.
    2161         if(!input)return;
    2162 
    2163         /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result
    2164          * object out of the input, with the additional step and time information: */
    2165         this->results->AddResult(input->SpawnResult(step,time));
    2166         #ifdef _HAVE_CONTROL_
    2167         if(input->ObjectEnum()==ControlInputEnum){
    2168                 if(((ControlInput*)input)->gradient!=NULL) this->results->AddResult(((ControlInput*)input)->SpawnGradient(step,time));
    2169         }
    2170         #endif
    2171 }
    2172 /*}}}*/
    21732075/*FUNCTION Penta::InputUpdateFromConstant(bool value, int name);{{{*/
    21742076void  Penta::InputUpdateFromConstant(bool constant, int name){
    … …  
    22112113        IssmDouble  yts;
    22122114        bool    control_analysis;
    2213         int     num_control_type;
    2214         int     num_cm_responses;
     2115        int     num_control_type,num_responses;
    22152116
    22162117        /*Fetch parameters: */
    … …  
    22182119        iomodel->Constant(&control_analysis,InversionIscontrolEnum);
    22192120        if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
    2220         if(control_analysis) iomodel->Constant(&num_cm_responses,InversionNumCostFunctionsEnum);
     2121        if(control_analysis) iomodel->Constant(&num_responses,InversionNumCostFunctionsEnum);
    22212122
    22222123        /*Checks if debuging*/
    … …  
    22672168                                        }
    22682169                                        break;
    2269                                 case MaterialsRheologyBbarEnum:
    2270                                 case MaterialsRheologyZbarEnum:
    2271                                         /*Material will take care of it*/ break;
     2170                                /*Material will take care of it*/
     2171                                case MaterialsRheologyBbarEnum: break;
     2172                                case DamageDbarEnum:break;
    22722173                                default:
    22732174                                        _error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet");
    … …  
    22772178        #endif
    22782179
    2279         //Need to know the type of approximation for this element
     2180        /*Need to know the type of approximation for this element*/
    22802181        if(iomodel->Data(FlowequationElementEquationEnum)){
    2281                 if (iomodel->Data(FlowequationElementEquationEnum)[index]==SSAApproximationEnum){
    2282                         this->inputs->AddInput(new IntInput(ApproximationEnum,SSAApproximationEnum));
    2283                 }
    2284                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOApproximationEnum){
    2285                         this->inputs->AddInput(new IntInput(ApproximationEnum,HOApproximationEnum));
    2286                 }
    2287                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SSAHOApproximationEnum){
    2288                         this->inputs->AddInput(new IntInput(ApproximationEnum,SSAHOApproximationEnum));
    2289                 }
    2290                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SIAApproximationEnum){
    2291                         this->inputs->AddInput(new IntInput(ApproximationEnum,SIAApproximationEnum));
    2292                 }
    2293                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==L1L2ApproximationEnum){
    2294                         this->inputs->AddInput(new IntInput(ApproximationEnum,L1L2ApproximationEnum));
    2295                 }
    2296                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==FSApproximationEnum){
    2297                         this->inputs->AddInput(new IntInput(ApproximationEnum,FSApproximationEnum));
    2298                 }
    2299                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SSAFSApproximationEnum){
    2300                         this->inputs->AddInput(new IntInput(ApproximationEnum,SSAFSApproximationEnum));
    2301                 }
    2302                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOFSApproximationEnum){
    2303                         this->inputs->AddInput(new IntInput(ApproximationEnum,HOFSApproximationEnum));
    2304                 }
    2305                 else if (iomodel->Data(FlowequationElementEquationEnum)[index]==NoneApproximationEnum){
    2306                         this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum));
    2307                 }
    2308                 else{
    2309                         _error_("Approximation type " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(FlowequationElementEquationEnum)[index])) << " not supported yet");
    2310                 }
     2182                this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
    23112183        }
    23122184
    … …  
    23162188                /*Create inputs and add to DataSetInput*/
    23172189                DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
    2318                 for(i=0;i<num_cm_responses;i++){
    2319                         for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_cm_responses+i];
    2320                         datasetinput->inputs->AddObject(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum));
     2190                for(i=0;i<num_responses;i++){
     2191                        for(j=0;j<6;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(penta_vertex_ids[j]-1)*num_responses+i];
     2192                        datasetinput->AddInput(new PentaInput(InversionCostFunctionsCoefficientsEnum,nodeinputs,P1Enum),reCast<int>(iomodel->Data(InversionCostFunctionsEnum)[i]));
    23212193                }
    23222194
    … …  
    23292201void  Penta::InputUpdateFromSolution(IssmDouble* solution){
    23302202
    2331         int analysis_type;
     2203        int analysis_type,inputenum;
    23322204
    23332205        /*retreive parameters: */
    … …  
    23702242                break;
    23712243        #endif
    2372         case BedSlopeXAnalysisEnum:
    2373                 InputUpdateFromSolutionOneDofCollapsed(solution,BedSlopeXEnum);
    2374                 break;
    2375         case BedSlopeYAnalysisEnum:
    2376                 InputUpdateFromSolutionOneDofCollapsed(solution,BedSlopeYEnum);
    2377                 break;
    2378         case SurfaceSlopeXAnalysisEnum:
    2379                 InputUpdateFromSolutionOneDofCollapsed(solution,SurfaceSlopeXEnum);
    2380                 break;
    2381         case SurfaceSlopeYAnalysisEnum:
    2382                 InputUpdateFromSolutionOneDofCollapsed(solution,SurfaceSlopeYEnum);
     2244        case L2ProjectionBaseAnalysisEnum:
     2245                this->parameters->FindParam(&inputenum,InputToL2ProjectEnum);
     2246                InputUpdateFromSolutionOneDofCollapsed(solution,inputenum);
    23832247                break;
    23842248        case MasstransportAnalysisEnum:
    … …  
    26762540                        }
    26772541                        /*update input*/
    2678                         if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    2679                                 material->inputs->AddInput(new PentaInput(name,values,P1Enum));
    2680                         }
    2681                         else{
    2682                                 this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    2683                         }
     2542                        this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    26842543                        return;
    26852544
    … …  
    26892548                        }
    26902549                        /*update input*/
    2691                         if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
    2692                                 material->inputs->AddInput(new PentaInput(name,values,P1Enum));
    2693                         }
    2694                         else{
    2695                                 this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    2696                         }
     2550                        this->inputs->AddInput(new PentaInput(name,values,P1Enum));
    26972551                        return;
    26982552
    … …  
    27442598                                name==SurfaceEnum ||
    27452599                                name==BedEnum ||
     2600                                name==BathymetryEnum ||
    27462601                                name==SurfaceSlopeXEnum ||
    27472602                                name==SurfaceSlopeYEnum ||
    … …  
    27822637                                name==QmuTemperatureEnum ||
    27832638                                name==QmuMeltingEnum ||
     2639                                name==QmuMaskGroundediceLevelsetEnum ||
     2640                                name==QmuMaskIceLevelsetEnum ||
    27842641                                name==GiaWEnum ||
    27852642                                name==GiadWdtEnum ||
    … …  
    28542711}
    28552712/*}}}*/
    2856 /*FUNCTION Penta::ListResultsInfo{{{*/
    2857 void Penta::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,IssmDouble** in_resultstimes,int** in_resultssteps,int* in_num_results){
    2858 
    2859         /*Intermediaries*/
    2860         int     i;
    2861         int     numberofresults = 0;
    2862         int     *resultsenums   = NULL;
    2863         int     *resultssizes   = NULL;
    2864         IssmDouble  *resultstimes   = NULL;
    2865         int     *resultssteps   = NULL;
    2866 
    2867         /*Checks*/
    2868         _assert_(in_num_results);
    2869 
    2870         /*Count number of results*/
    2871         for(i=0;i<this->results->Size();i++){
    2872                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2873                 numberofresults++;
    2874         }
    2875 
    2876         if(numberofresults){
    2877 
    2878                 /*Allocate output*/
    2879                 resultsenums=xNew<int>(numberofresults);
    2880                 resultssizes=xNew<int>(numberofresults);
    2881                 resultstimes=xNew<IssmDouble>(numberofresults);
    2882                 resultssteps=xNew<int>(numberofresults);
    2883 
    2884                 /*populate enums*/
    2885                 for(i=0;i<this->results->Size();i++){
    2886                         ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2887                         resultsenums[i]=elementresult->InstanceEnum();
    2888                         resultstimes[i]=elementresult->GetTime();
    2889                         resultssteps[i]=elementresult->GetStep();
    2890                         if(elementresult->ObjectEnum()==PentaP1ElementResultEnum){
    2891                                 resultssizes[i]=P1Enum;
    2892                         }
    2893                         else{
    2894                                 resultssizes[i]=P0Enum;
    2895                         }
    2896                 }
    2897         }
    2898 
    2899         /*Assign output pointers:*/
    2900         *in_num_results=numberofresults;
    2901         *in_resultsenums=resultsenums;
    2902         *in_resultssizes=resultssizes;
    2903         *in_resultstimes=resultstimes;
    2904         *in_resultssteps=resultssteps;
    2905 
    2906 }/*}}}*/
    29072713/*FUNCTION Penta::MinEdgeLength{{{*/
    29082714IssmDouble Penta::MinEdgeLength(IssmDouble xyz_list[6][3]){
    … …  
    29572763        if(found)*pvalue=value;
    29582764        return found;
    2959 }
    2960 /*}}}*/
    2961 /*FUNCTION Penta::PatchFill{{{*/
    2962 void  Penta::PatchFill(int* pcount, Patch* patch){
    2963 
    2964         int i,count;
    2965         int vertices_ids[6];
    2966 
    2967         /*recover pointer: */
    2968         count=*pcount;
    2969 
    2970         /*will be needed later: */
    2971         for(i=0;i<6;i++) vertices_ids[i]=vertices[i]->Id(); //vertices id start at column 3 of the patch.
    2972 
    2973         for(i=0;i<this->results->Size();i++){
    2974                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2975 
    2976                 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand
    2977                  *it to the result object, to fill the rest: */
    2978                 patch->fillelementinfo(count,this->sid+1,vertices_ids,6);
    2979                 elementresult->PatchFill(count,patch);
    2980 
    2981                 /*increment counter: */
    2982                 count++;
    2983         }
    2984 
    2985         /*Assign output pointers:*/
    2986         *pcount=count;
    2987 }/*}}}*/
    2988 /*FUNCTION Penta::PatchSize{{{*/
    2989 void  Penta::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
    2990 
    2991         int     i;
    2992         int     numrows       = 0;
    2993         int     numnodes      = 0;
    2994         int     temp_numnodes = 0;
    2995 
    2996         /*Go through all the results objects, and update the counters: */
    2997         for (i=0;i<this->results->Size();i++){
    2998                 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
    2999                 /*first, we have one more result: */
    3000                 numrows++;
    3001                 /*now, how many vertices and how many nodal values for this result? :*/
    3002                 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
    3003                 if(temp_numnodes>numnodes)numnodes=temp_numnodes;
    3004         }
    3005 
    3006         /*Assign output pointers:*/
    3007         *pnumrows=numrows;
    3008         *pnumvertices=NUMVERTICES;
    3009         *pnumnodes=numnodes;
    30102765}
    30112766/*}}}*/
    … …  
    30622817}
    30632818/*}}}*/
    3064 /*FUNCTION Penta::ReduceMatrixFS {{{*/
    3065 void Penta::ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){
    3066 
    3067         int    i,j;
    3068         IssmDouble Kii[24][24];
    3069         IssmDouble Kib[24][3];
    3070         IssmDouble Kbb[3][3];
    3071         IssmDouble Kbi[3][24];
    3072         IssmDouble Kbbinv[3][3];
    3073         IssmDouble Kright[24][24];
    3074 
    3075         /*Create the four matrices used for reduction */
    3076         for(i=0;i<24;i++){
    3077                 for(j=0;j<24;j++){
    3078                         Kii[i][j]=*(Ke_temp+27*i+j);
    3079                 }
    3080                 for(j=0;j<3;j++){
    3081                         Kib[i][j]=*(Ke_temp+27*i+j+24);
    3082                 }
    3083         }
    3084         for(i=0;i<3;i++){
    3085                 for(j=0;j<24;j++){
    3086                         Kbi[i][j]=*(Ke_temp+27*(i+24)+j);
    3087                 }
    3088                 for(j=0;j<3;j++){
    3089                         Kbb[i][j]=*(Ke_temp+27*(i+24)+j+24);
    3090                 }
    3091         }
    3092 
    3093         /*Inverse the matrix corresponding to bubble part Kbb */
    3094         Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
    3095 
    3096         /*Multiply matrices to create the reduce matrix Ke_reduced */
    3097         TripleMultiply(&Kib[0][0],24,3,0,
    3098                                 &Kbbinv[0][0],3,3,0,
    3099                                 &Kbi[0][0],3,24,0,
    3100                                 &Kright[0][0],0);
    3101 
    3102         /*Affect value to the reduced matrix */
    3103         for(i=0;i<24;i++) for(j=0;j<24;j++) *(Ke_reduced+24*i+j)=Kii[i][j]-Kright[i][j];
    3104 }
    3105 /*}}}*/
    3106 /*FUNCTION Penta::ReduceVectorFS {{{*/
    3107 void Penta::ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){
    3108 
    3109         int    i,j;
    3110         IssmDouble Pi[24];
    3111         IssmDouble Pb[3];
    3112         IssmDouble Kbb[3][3];
    3113         IssmDouble Kib[24][3];
    3114         IssmDouble Kbbinv[3][3];
    3115         IssmDouble Pright[24];
    3116 
    3117         /*Create the four matrices used for reduction */
    3118         for(i=0;i<24;i++) Pi[i]=*(Pe_temp+i);
    3119         for(i=0;i<3;i++) Pb[i]=*(Pe_temp+i+24);
    3120         for(j=0;j<3;j++){
    3121                 for(i=0;i<24;i++){
    3122                         Kib[i][j]=*(Ke_temp+3*i+j);
    3123                 }
    3124                 for(i=0;i<3;i++){
    3125                         Kbb[i][j]=*(Ke_temp+3*(i+24)+j);
    3126                 }
    3127         }
    3128 
    3129         /*Inverse the matrix corresponding to bubble part Kbb */
    3130         Matrix3x3Invert(&Kbbinv[0][0], &Kbb[0][0]);
    3131 
    3132         /*Multiply matrices to create the reduce matrix Ke_reduced */
    3133         TripleMultiply(&Kib[0][0],24,3,0,
    3134                                 &Kbbinv[0][0],3,3,0,
    3135                                 &Pb[0],3,1,0,&Pright[0],0);
    3136 
    3137         /*Affect value to the reduced matrix */
    3138         for(i=0;i<24;i++) *(Pe_reduced+i)=Pi[i]-Pright[i];
    3139 }
    3140 /*}}}*/
    3141 /*FUNCTION Penta::RequestedOutput{{{*/
    3142 void Penta::RequestedOutput(int output_enum,int step,IssmDouble time){
    3143         if(IsInput(output_enum)){
    3144                 /*just transfer this input to results, and we are done: */
    3145                 InputToResult(output_enum,step,time);
    3146         }
    3147         else{
    3148                 /*this input does not exist, compute it, and then transfer to results: */
     2819/*FUNCTION Penta::ResultInterpolation{{{*/
     2820void Penta::ResultInterpolation(int* pinterpolation,int output_enum){
     2821
     2822        Input* input=this->inputs->GetInput(output_enum);
     2823
     2824        /*If this input is not already in Inputs, maybe it needs to be computed?*/
     2825        if(!input){
    31492826                switch(output_enum){
    3150                         case BasalFrictionEnum:
    3151 
    3152                                 /*create input: */
    3153                                 BasalFrictionCreateInput();
    3154 
    3155                                 /*transfer to results :*/
    3156                                 InputToResult(output_enum,step,time);
    3157 
    3158                                 /*erase input: */
    3159                                 inputs->DeleteInput(output_enum);
     2827                        case ViscousHeatingEnum:
     2828                                this->ViscousHeatingCreateInput();
     2829                                input=this->inputs->GetInput(output_enum);
    31602830                                break;
    3161                         case ViscousHeatingEnum:
    3162 
    3163                                 /*create input: */
    3164                                 ViscousHeatingCreateInput();
    3165 
    3166                                 /*transfer to results :*/
    3167                                 InputToResult(output_enum,step,time);
    3168 
    3169                                 /*erase input: */
    3170                                 inputs->DeleteInput(output_enum);
     2831                        case StressTensorxxEnum:
     2832                                this->ComputeStressTensor();
     2833                                input=this->inputs->GetInput(output_enum);
    31712834                                break;
    3172 
    3173                         case StressTensorEnum:
     2835                        case StressTensorxyEnum:
    31742836                                this->ComputeStressTensor();
    3175                                 InputToResult(StressTensorxxEnum,step,time);
    3176                                 InputToResult(StressTensorxyEnum,step,time);
    3177                                 InputToResult(StressTensorxzEnum,step,time);
    3178                                 InputToResult(StressTensoryyEnum,step,time);
    3179                                 InputToResult(StressTensoryzEnum,step,time);
    3180                                 InputToResult(StressTensorzzEnum,step,time);
     2837                                input=this->inputs->GetInput(output_enum);
    31812838                                break;
    3182 
     2839                        case StressTensorxzEnum:
     2840                                this->ComputeStressTensor();
     2841                                input=this->inputs->GetInput(output_enum);
     2842                                break;
     2843                        case StressTensoryyEnum:
     2844                                this->ComputeStressTensor();
     2845                                input=this->inputs->GetInput(output_enum);
     2846                                break;
     2847                        case StressTensoryzEnum:
     2848                                this->ComputeStressTensor();
     2849                                input=this->inputs->GetInput(output_enum);
     2850                                break;
     2851                        case StressTensorzzEnum:
     2852                                this->ComputeStressTensor();
     2853                                input=this->inputs->GetInput(output_enum);
     2854                                break;
    31832855                        default:
    3184                                 /*do nothing, no need to derail the computation because one of the outputs requested cannot be found: */
    3185                                 break;
    3186                 }
    3187         }
     2856                                _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
     2857                }
     2858        }
     2859
     2860        /*Assign output pointer*/
     2861        *pinterpolation = input->GetResultInterpolation();
     2862
     2863}
     2864/*}}}*/
     2865/*FUNCTION Penta::ResultToVector{{{*/
     2866void Penta::ResultToVector(Vector<IssmPDouble>* vector,int output_enum){
     2867
     2868        Input* input=this->inputs->GetInput(output_enum);
     2869        if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
     2870
     2871        switch(input->GetResultInterpolation()){
     2872                case P0Enum:
     2873                        _error_("not implemented...");
     2874                        break;
     2875                case P1Enum:{
     2876                        IssmDouble  values[NUMVERTICES];
     2877                        IssmPDouble pvalues[NUMVERTICES];
     2878                        int         connectivity[NUMVERTICES];
     2879                        int         sidlist[NUMVERTICES];
     2880
     2881                        this->GetVertexSidList(&sidlist[0]);
     2882                        this->GetConnectivityList(&connectivity[0]);
     2883                        this->GetInputListOnVertices(&values[0],output_enum);
     2884                        for(int i=0;i<NUMVERTICES;i++) pvalues[i] = reCast<IssmPDouble>(values[i])/reCast<IssmPDouble>(connectivity[i]);
     2885
     2886                        vector->SetValues(NUMVERTICES,&sidlist[0],&pvalues[0],ADD_VAL);
     2887                        break;
     2888                                        }
     2889                default:
     2890                        _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet");
     2891        }
     2892
    31882893
    31892894}
    … …  
    33403045
    33413046        /*Spawn material*/
    3342         tria->material=NULL;
    33433047        tria->material=(Material*)this->material->copy();
    33443048        delete tria->material->inputs;
    … …  
    35153219        bool       dakota_analysis;
    35163220        bool       isFS;
    3517         IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat;
    35183221        int        numnodes;
    35193222        int*       penta_node_ids = NULL;
    … …  
    35233226        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    35243227        iomodel->Constant(&isFS,FlowequationIsFSEnum);
    3525         iomodel->Constant(&beta,MaterialsBetaEnum);
    3526         iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
    3527         iomodel->Constant(&referencetemperature,ConstantsReferencetemperatureEnum);
    3528         iomodel->Constant(&meltingpoint,MaterialsMeltingpointEnum);
    3529         iomodel->Constant(&latentheat,MaterialsLatentheatEnum);
    35303228
    35313229        /*Checks if debuging*/
    3532         /*{{{*/
    35333230        _assert_(iomodel->elements);
    3534         /*}}}*/
    35353231
    35363232        /*Recover element type*/
    … …  
    36883384                case StressbalanceAnalysisEnum:
    36893385
    3690                         /*default vx,vy and vz: either observation or 0 */
    3691                         if(!iomodel->Data(VxEnum)){
    3692                                 for(i=0;i<6;i++)nodeinputs[i]=0;
    3693                                 this->inputs->AddInput(new PentaInput(VxEnum,nodeinputs,P1Enum));
    3694                                 if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVxEnum,nodeinputs,P1Enum));
    3695                         }
    3696                         if(!iomodel->Data(VyEnum)){
    3697                                 for(i=0;i<6;i++)nodeinputs[i]=0;
    3698                                 this->inputs->AddInput(new PentaInput(VyEnum,nodeinputs,P1Enum));
    3699                                 if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVyEnum,nodeinputs,P1Enum));
    3700                         }
    3701                         if(!iomodel->Data(VzEnum)){
    3702                                 for(i=0;i<6;i++)nodeinputs[i]=0;
    3703                                 this->inputs->AddInput(new PentaInput(VzEnum,nodeinputs,P1Enum));
    3704                                 if(dakota_analysis) this->inputs->AddInput(new PentaInput(QmuVzEnum,nodeinputs,P1Enum));
    3705                         }
    3706                         if(!iomodel->Data(PressureEnum)){
    3707                                 for(i=0;i<6;i++)nodeinputs[i]=0;
    3708                                 if(dakota_analysis){
    3709                                         this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum));
    3710                                         this->inputs->AddInput(new PentaInput(QmuPressureEnum,nodeinputs,P1Enum));
    3711                                 }
    3712                                 if(isFS){
    3713                                         this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum));
    3714                                         this->inputs->AddInput(new PentaInput(PressurePicardEnum,nodeinputs,P1Enum));
    3715                                 }
    3716                         }
    37173386                        if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==HOFSApproximationEnum){
    37183387                                /*Create VzHO and VzFS Enums*/
    … …  
    37443413                        }
    37453414                        break;
    3746 
    3747                 case ThermalAnalysisEnum:
    3748                         /*Initialize mesh velocity*/
    3749                         for(i=0;i<6;i++)nodeinputs[i]=0;
    3750                         this->inputs->AddInput(new PentaInput(VxMeshEnum,nodeinputs,P1Enum));
    3751                         this->inputs->AddInput(new PentaInput(VyMeshEnum,nodeinputs,P1Enum));
    3752                         this->inputs->AddInput(new PentaInput(VzMeshEnum,nodeinputs,P1Enum));
    3753                         if(dakota_analysis){
    3754                                 this->inputs->AddInput(new PentaInput(QmuVxMeshEnum,nodeinputs,P1Enum));
    3755                                 this->inputs->AddInput(new PentaInput(QmuVyMeshEnum,nodeinputs,P1Enum));
    3756                                 this->inputs->AddInput(new PentaInput(QmuVzMeshEnum,nodeinputs,P1Enum));
    3757                         }
    3758                         break;
    3759 
    3760                 case EnthalpyAnalysisEnum:
    3761                         /*Initialize mesh velocity*/
    3762                         for(i=0;i<6;i++)nodeinputs[i]=0;
    3763                         this->inputs->AddInput(new PentaInput(VxMeshEnum,nodeinputs,P1Enum));
    3764                         this->inputs->AddInput(new PentaInput(VyMeshEnum,nodeinputs,P1Enum));
    3765                         this->inputs->AddInput(new PentaInput(VzMeshEnum,nodeinputs,P1Enum));
    3766                         if (iomodel->Data(TemperatureEnum) && iomodel->Data(WaterfractionEnum) && iomodel->Data(PressureEnum)) {
    3767                                 for(i=0;i<6;i++){
    3768                                         if(iomodel->Data(TemperatureEnum)[penta_vertex_ids[i]-1] < meltingpoint-beta*iomodel->Data(PressureEnum)[penta_vertex_ids[i]-1]){
    3769                                                 nodeinputs[i]=heatcapacity*(iomodel->Data(TemperatureEnum)[penta_vertex_ids[i]-1]-referencetemperature);
    3770                                         }
    3771                                         else nodeinputs[i]=heatcapacity*
    3772                                          (meltingpoint-beta*iomodel->Data(PressureEnum)[penta_vertex_ids[i]-1]-referencetemperature)
    3773                                                 +latentheat*iomodel->Data(WaterfractionEnum)[penta_vertex_ids[i]-1];
    3774                                 }
    3775                                 this->inputs->AddInput(new PentaInput(EnthalpyEnum,nodeinputs,P1Enum));
    3776                         }
    3777                         else _error_("temperature and waterfraction required for the enthalpy solution");
    3778                         break;
    3779 
    37803415                default:
    37813416                        /*No update for other solution types*/
    … …  
    37863421/*FUNCTION Penta::ViscousHeatingCreateInput {{{*/
    37873422void Penta::ViscousHeatingCreateInput(void){
    3788 
    3789         /*Constants*/
    3790         const int  numdof=NUMVERTICES*NDOF1;
    37913423
    37923424        /*Intermediaries*/
    … …  
    37953427        IssmDouble xyz_list[NUMVERTICES][3];
    37963428        IssmDouble epsilon[6];
    3797         IssmDouble viscousheating[NUMVERTICES]={0,0,0,0,0,0};
     3429        IssmDouble viscousheating[NUMVERTICES];
    37983430        IssmDouble thickness;
    3799         GaussPenta *gauss=NULL;
    3800 
    3801         /*Initialize Element vector*/
    3802         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    38033431
    38043432        /*Retrieve all inputs and parameters*/
    … …  
    38103438
    38113439        /*loop over vertices: */
    3812         gauss=new GaussPenta();
     3440        GaussPenta* gauss=new GaussPenta();
    38133441        for (int iv=0;iv<NUMVERTICES;iv++){
    38143442                gauss->GaussVertex(iv);
    … …  
    38603488}
    38613489/*}}}*/
     3490/*FUNCTION Penta::IceVolumeAboveFloatation {{{*/
     3491IssmDouble Penta::IceVolumeAboveFloatation(void){
     3492
     3493        /*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
     3494        IssmDouble rho_ice,rho_water;
     3495        IssmDouble base,bed,surface,bathymetry;
     3496        IssmDouble xyz_list[NUMVERTICES][3];
     3497
     3498        if(NoIceInElement() || IsFloating() || !IsOnBed())return 0;
     3499
     3500        rho_ice=matpar->GetRhoIce();
     3501        rho_water=matpar->GetRhoWater();
     3502        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     3503
     3504        /*First calculate the area of the base (cross section triangle)
     3505         * http://en.wikipedia.org/wiki/Pentangle
     3506         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     3507        base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     3508
     3509        /*Now get the average height above floatation*/
     3510        Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
     3511        Input* bed_input        = inputs->GetInput(BedEnum);        _assert_(bed_input);
     3512        Input* bathymetry_input = inputs->GetInput(BathymetryEnum); _assert_(bathymetry_input);
     3513        surface_input->GetInputAverage(&surface);
     3514        bed_input->GetInputAverage(&bed);
     3515        bathymetry_input->GetInputAverage(&bathymetry);
     3516
     3517        /*Return: */
     3518        return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
     3519}
     3520/*}}}*/
    38623521/*FUNCTION Penta::MinVel{{{*/
    38633522void  Penta::MinVel(IssmDouble* pminvel){
    … …  
    39243583}
    39253584/*}}}*/
     3585/*FUNCTION Penta::MassFlux {{{*/
     3586IssmDouble Penta::MassFlux( IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){
     3587
     3588        IssmDouble mass_flux=0;
     3589
     3590        if(!IsOnBed()) return mass_flux;
     3591
     3592        /*Depth Averaging Vx and Vy*/
     3593        this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
     3594        this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
     3595
     3596        /*Spawn Tria element from the base of the Penta: */
     3597        Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
     3598        mass_flux=tria->MassFlux(x1,y1,x2,y2,segment_id);
     3599        delete tria->material; delete tria;
     3600
     3601        /*Delete Vx and Vy averaged*/
     3602        this->inputs->DeleteInput(VxAverageEnum);
     3603        this->inputs->DeleteInput(VyAverageEnum);
     3604
     3605        /*clean up and return*/
     3606        return mass_flux;
     3607}
     3608/*}}}*/
    39263609/*FUNCTION Penta::MaxAbsVx{{{*/
    39273610void  Penta::MaxAbsVx(IssmDouble* pmaxabsvx){
    … …  
    40023685                        *presponse=this->material->GetBbar();
    40033686                        break;
    4004                 case MaterialsRheologyZbarEnum:
    4005                         *presponse=this->material->GetZbar();
     3687                case DamageDbarEnum:
     3688                        *presponse=this->material->GetDbar();
    40063689                        break;
    40073690                case VelEnum:
    … …  
    40883771        /*Intermediaries */
    40893772        int        stabilization;
    4090         int        i,j,found=0;
     3773        int        i,j;
    40913774        IssmDouble Jdet,u,v,w,um,vm,wm;
    40923775        IssmDouble h,hx,hy,hz,vx,vy,vz,vel;
    … …  
    41073790        IssmDouble D[3][3];
    41083791        IssmDouble K[3][3]={0.0};
    4109         Tria*      tria=NULL;
    41103792        GaussPenta *gauss=NULL;
    41113793
    … …  
    41423824
    41433825                /*Conduction: */ 
    4144                 /*Need to change that depending on enthalpy value -> cold or temperate ice: */ 
    41453826                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    41463827
    41473828                enthalpy_input->GetInputValue(&enthalpy, gauss);
    41483829                pressure_input->GetInputValue(&pressure, gauss);
    4149                 kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
    4150                 D_scalar_conduct=gauss->weight*Jdet*kappa;
     3830                kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); _assert_(kappa>0.);
     3831                D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice;
    41513832                if(reCast<bool,IssmDouble>(dt)) D_scalar_conduct=D_scalar_conduct*dt;
    41523833
    … …  
    41913872                }
    41923873
    4193                 /*Artifficial diffusivity*/
     3874                /*Artificial diffusivity*/
    41943875                if(stabilization==1){
    41953876                        /*Build K: */
    … …  
    41973878                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
    41983879                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
     3880
    41993881                        K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
    42003882                        K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
    42013883                        K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
     3884
    42023885                        D_scalar_stab=gauss->weight*Jdet;
    42033886                        if(reCast<bool,IssmDouble>(dt)) D_scalar_stab=D_scalar_stab*dt;
    … …  
    43244007        /*Intermediaries */
    43254008        int        stabilization;
    4326         int        i,j,found=0;
     4009        int        i,j;
    43274010        IssmDouble Jdet,u,v,w,um,vm,wm,vel;
    43284011        IssmDouble h,hx,hy,hz,vx,vy,vz;
    … …  
    43404023        IssmDouble D[3][3];
    43414024        IssmDouble K[3][3]={0.0};
    4342         Tria*      tria=NULL;
    43434025        GaussPenta *gauss=NULL;
    43444026
    … …  
    43894071
    43904072                /*Advection: */
    4391 
    43924073                GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss);
    43934074                GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    … …  
    45444225
    45454226        /*Intermediaries*/
    4546         int    i,j,found=0;
    4547         int    friction_type,stabilization;
     4227        int        i;
     4228        int        stabilization;
    45484229        IssmDouble Jdet,phi,dt;
    45494230        IssmDouble rho_ice,heatcapacity;
    … …  
    45544235        IssmDouble u,v,w;
    45554236        IssmDouble scalar_def,scalar_transient;
    4556         IssmDouble temperature_list[NUMVERTICES];
    45574237        IssmDouble xyz_list[NUMVERTICES][3];
    45584238        IssmDouble L[numdof];
    … …  
    46004280                scalar_def=phi/rho_ice*Jdet*gauss->weight;
    46014281                if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
    4602 
     4282               
     4283                /*TODO: add -beta*laplace T_m(p)*/
    46034284                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
    46044285
    … …  
    46194300                        enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
    46204301                        kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
    4621                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
     4302                        kappa/=rho_ice;
     4303                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
    46224304
    46234305                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    … …  
    47024384        IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    47034385        IssmDouble Jdet2d,dt;
    4704         IssmDouble rho_ice,heatcapacity,geothermalflux_value;
    4705         IssmDouble basalfriction,alpha2,vx,vy;
     4386        IssmDouble rho_ice,heatcapacity;
     4387        IssmDouble geothermalflux_value, heatflux_value;
     4388        IssmDouble basalfriction,alpha2,vx,vy,vz;
    47064389        IssmDouble scalar,enthalpy,enthalpyup;
    47074390        IssmDouble pressure,pressureup;
    … …  
    47394422        gauss=new GaussPenta(0,1,2,2);
    47404423        gaussup=new GaussPenta(3,4,5,2);
     4424
    47414425        for(int ig=gauss->begin();ig<gauss->end();ig++){
    47424426
    … …  
    47494433                enthalpy_input->GetInputValue(&enthalpy,gauss);
    47504434                pressure_input->GetInputValue(&pressure,gauss);
    4751                 if(enthalpy>=matpar->PureIceEnthalpy(pressure)){
     4435                watercolumn_input->GetInputValue(&watercolumn,gauss);
     4436
     4437                if((watercolumn<=0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
     4438                        /* the above check is equivalent to
     4439                         NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
     4440                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     4441
     4442                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4443                        vx_input->GetInputValue(&vx,gauss);
     4444                        vy_input->GetInputValue(&vy,gauss);
     4445                        vz_input->GetInputValue(&vz,gauss);
     4446                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
     4447
     4448                        heatflux_value=(basalfriction+geothermalflux_value)/(rho_ice);
     4449                        scalar=gauss->weight*Jdet2d*heatflux_value;
     4450                        if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
     4451                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     4452                }
     4453                else if(enthalpy>=matpar->PureIceEnthalpy(pressure)){
     4454                        /* check positive thickness of temperate basal ice layer */
    47524455                        enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    47534456                        pressure_input->GetInputValue(&pressureup,gaussup);
    47544457                        if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
    4755                                 // temperate ice has positive thickness: grad enthalpy*n=0.
     4458                                // TODO: temperate ice has positive thickness: grad enthalpy*n=0.
    47564459                        }
    47574460                        else{
    4758                                 // only base temperate, set Dirichlet BCs in Penta::UpdateThermalBasalConstraints()
     4461                                // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
    47594462                        }
    47604463                }
    4761                 else{
    4762                         watercolumn_input->GetInputValue(&watercolumn,gauss);
    4763                         if(watercolumn==0.){
    4764                                 /*add geothermal heat flux and basal friction*/
    4765                                 geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4766                                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    4767                                 vx_input->GetInputValue(&vx,gauss);
    4768                                 vy_input->GetInputValue(&vy,gauss);
    4769                                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    4770 
    4771                                 scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
    4772                                 if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
    4773 
    4774                                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    4775                         }
    4776                         else{  /*do nothing (water layer acts as insulation)*/  }
     4464                else {
     4465                        // base cold, but watercolumn positive. Set base to h_pmp.
    47774466                }
    47784467        }
    … …  
    48134502
    48144503        /*Intermediaries*/
    4815         int    i,j,found=0;
    4816         int    friction_type,stabilization;
     4504        int        i;
     4505        int        stabilization;
    48174506        IssmDouble Jdet,phi,dt;
    48184507        IssmDouble rho_ice,heatcapacity;
    … …  
    48224511        IssmDouble u,v,w;
    48234512        IssmDouble scalar_def,scalar_transient;
    4824         IssmDouble temperature_list[NUMVERTICES];
    48254513        IssmDouble xyz_list[NUMVERTICES][3];
    48264514        IssmDouble L[numdof];
    … …  
    48444532        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    48454533        Input* temperature_input=NULL;
    4846         if (reCast<bool,IssmDouble>(dt)) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
    4847         if (stabilization==2) diameter=MinEdgeLength(xyz_list);
     4534        if(reCast<bool,IssmDouble>(dt)) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
     4535        if(stabilization==2) diameter=MinEdgeLength(xyz_list);
    48484536
    48494537        /* Start  looping on the number of gaussian points: */
    … …  
    49574645
    49584646        /*Intermediaries */
    4959         int        i,j;
    4960         int        analysis_type;
    4961         IssmDouble xyz_list[NUMVERTICES][3];
    4962         IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    4963         IssmDouble Jdet2d,dt;
    4964         IssmDouble rho_ice,heatcapacity,geothermalflux_value;
    4965         IssmDouble basalfriction,alpha2,vx,vy;
    4966         IssmDouble basis[NUMVERTICES];
    4967         IssmDouble scalar;
    4968         Friction*  friction=NULL;
     4647        int         i,j;
     4648        int         analysis_type;
     4649        IssmDouble  xyz_list[NUMVERTICES][3];
     4650        IssmDouble  xyz_list_tria[NUMVERTICES2D][3]={0.0};
     4651        IssmDouble  Jdet2d,dt;
     4652        IssmDouble  rho_ice,heatcapacity,geothermalflux_value;
     4653        IssmDouble  basalfriction,alpha2,vx,vy;
     4654        IssmDouble  basis[NUMVERTICES];
     4655        IssmDouble  scalar;
     4656        Friction*   friction=NULL;
    49694657        GaussPenta* gauss=NULL;
    49704658
    … …  
    50174705}
    50184706/*}}}*/
    5019 /*FUNCTION Penta::GetSolutionFromInputsThermal{{{*/
    5020 void  Penta::GetSolutionFromInputsThermal(Vector<IssmDouble>* solution){
    5021 
    5022         const int    numdof=NDOF1*NUMVERTICES;
    5023 
    5024         int         i;
    5025         int        *doflist = NULL;
    5026         IssmDouble  values[numdof];
    5027         IssmDouble  temp;
    5028         GaussPenta *gauss = NULL;
    5029 
    5030         /*Get dof list: */
    5031         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    5032         Input* t_input=inputs->GetInput(TemperatureEnum); _assert_(t_input);
    5033 
    5034         gauss=new GaussPenta();
    5035         for(i=0;i<NUMVERTICES;i++){
    5036                 /*Recover temperature*/
    5037                 gauss->GaussVertex(i);
    5038                 t_input->GetInputValue(&temp,gauss);
    5039                 values[i]=temp;
    5040         }
    5041 
    5042         /*Add value to global vector*/
    5043         solution->SetValues(numdof,doflist,values,INS_VAL);
    5044 
    5045         /*Free ressources:*/
    5046         delete gauss;
    5047         xDelete<int>(doflist);
    5048 }
    5049 /*}}}*/
    5050 /*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{*/
    5051 void  Penta::GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solution){
    5052 
    5053         const int    numdof=NDOF1*NUMVERTICES;
    5054 
    5055         int*        doflist=NULL;
    5056         IssmDouble  values[numdof];
    5057         IssmDouble  enthalpy;
    5058         GaussPenta *gauss=NULL;
    5059 
    5060         /*Get dof list: */
    5061         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    5062         Input* h_input=inputs->GetInput(EnthalpyEnum); _assert_(h_input);
    5063 
    5064         gauss=new GaussPenta();
    5065         for(int i=0;i<NUMVERTICES;i++){
    5066                 /*Recover temperature*/
    5067                 gauss->GaussVertex(i);
    5068                 h_input->GetInputValue(&enthalpy,gauss);
    5069                 values[i]=enthalpy;
    5070         }
    5071 
    5072         /*Add value to global vector*/
    5073         solution->SetValues(numdof,doflist,values,INS_VAL);
    5074 
    5075         /*Free ressources:*/
    5076         delete gauss;
    5077         xDelete<int>(doflist);
    5078 }
    5079 /*}}}*/
    50804707/*FUNCTION Penta::InputUpdateFromSolutionThermal {{{*/
    50814708void  Penta::InputUpdateFromSolutionThermal(IssmDouble* solution){
    … …  
    50904717        IssmDouble  B_average,s_average;
    50914718        int        *doflist = NULL;
    5092         IssmDouble  pressure[numdof];
     4719        bool        hack    = false;
    50934720
    50944721        /*Get dof list: */
    … …  
    50984725        for(i=0;i<numdof;i++){
    50994726                values[i]=solution[doflist[i]];
    5100                 //GetInputListOnVertices(&pressure[0],PressureEnum);
    5101                 //if(values[i]>matpar->TMeltingPoint(pressure[i])) values[i]=matpar->TMeltingPoint(pressure[i]);
    5102                 //if(values[i]<matpar->TMeltingPoint(pressure[i])-50) values[i]=matpar->TMeltingPoint(pressure[i])-50;
    51034727
    51044728                /*Check solution*/
    … …  
    51074731                //if(values[i]>275)    _printf_("temperature > 275°K found in solution vector (Paterson's rheology associated is negative)\n");
    51084732        }
     4733
     4734        if(hack){
     4735                /*Force temperature between [Tpmp-50 Tpmp] to disable penalties*/
     4736                IssmDouble pressure[numdof];
     4737                GetInputListOnVertices(&pressure[0],PressureEnum);
     4738                for(i=0;i<numdof;i++){
     4739                        if(values[i]>matpar->TMeltingPoint(pressure[i]))     values[i]=matpar->TMeltingPoint(pressure[i]);
     4740                        if(values[i]<matpar->TMeltingPoint(pressure[i])-50.) values[i]=matpar->TMeltingPoint(pressure[i])-50.;
     4741                }
     4742        }
     4743
    51094744
    51104745        /*Get all inputs and parameters*/
    … …  
    52164851                        case LliboutryDuvalEnum:
    52174852                                B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,
    5218                                                         (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0,
    5219                                                         material->GetN());
     4853                                                                                                                                 (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0,
     4854                                                                                                                                 material->GetN(),
     4855                                                                                                                                 matpar->GetBeta(), matpar->GetReferenceTemperature(), matpar->GetHeatCapacity(), matpar->GetLatentHeat());
    52204856                                for(i=0;i<numdof;i++) B[i]=B_average;
    52214857                                this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
    … …  
    52334869}
    52344870/*}}}*/
    5235 /*FUNCTION Penta::UpdateThermalBasalConstraints{{{*/
    5236 void  Penta::UpdateThermalBasalConstraints(void){
     4871/*FUNCTION Penta::UpdateBasalConstraintsEnthalpy{{{*/
     4872void  Penta::UpdateBasalConstraintsEnthalpy(void){
    52374873
    52384874        /*Intermediary*/
    5239         bool        isenthalpy,isdynamicbasalspc,istemperatelayer;
     4875        bool        isenthalpy,isdynamicbasalspc,setspc;
    52404876        int         numindices, numindicesup;
    5241         IssmDouble  h_pmp,pressure, pressureup;
    5242     IssmDouble  enthalpy, enthalpyup;
     4877        IssmDouble  pressure, pressureup;
     4878        IssmDouble  h_pmp, enthalpy, enthalpyup;
     4879        IssmDouble  watercolumn;
    52434880        int        *indices = NULL, *indicesup = NULL;
    52444881
    … …  
    52494886        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    52504887        if(!isenthalpy) return;
    5251         //parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
    5252         isdynamicbasalspc = true;
     4888        parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum);
    52534889        if(!isdynamicbasalspc) return;
    52544890
    52554891        /*Fetch indices of basal & surface nodes for this finite element*/
    52564892        BasalNodeIndices(&numindices,&indices,this->element_type);
    5257     SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
    5258     _assert_(numindices==numindicesup);
     4893        SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type);
     4894        _assert_(numindices==numindicesup);
    52594895
    52604896        /*Get parameters and inputs: */
    52614897        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    5262     Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
    5263 
    5264         /*if there is a temperate layer of positive thickness, set enthalpy=h_pmp at that node*/
     4898        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
     4899        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); //_assert_(watercolumn_input);
     4900
     4901        /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/
    52654902        GaussPenta* gauss=new GaussPenta();
    5266     GaussPenta* gaussup=new GaussPenta();
     4903        GaussPenta* gaussup=new GaussPenta();
    52674904        for(int i=0;i<numindices;i++){
    52684905                gauss->GaussNode(this->element_type,indices[i]);
    5269         gaussup->GaussNode(this->element_type,indicesup[i]); // TODO: check: are the nodes corresponding?
    5270 
    5271                 /*Check wether there is a temperate layer at the base or not -> TODO: Johannes:) */
    5272         /*check if node is temperate, if not, return*/
    5273         enthalpy_input->GetInputValue(&enthalpy, gauss);
     4906                gaussup->GaussNode(this->element_type,indicesup[i]);
     4907
     4908                /*Check wether there is a temperate layer at the base or not */
     4909                /*check if node is temperate, if not, continue*/
     4910                enthalpy_input->GetInputValue(&enthalpy, gauss);
    52744911                pressure_input->GetInputValue(&pressure, gauss);
    5275         if (enthalpy<matpar->PureIceEnthalpy(pressure)){
    5276           // TODO: reset, if necessary, all spcs to non-valid
    5277           continue;
    5278         }
    5279         /*check if upper node is temperate. if yes, then we have a temperate layer of positive thickness. if not, continue.*/
    5280         enthalpy_input->GetInputValue(&enthalpyup, gaussup);
    5281                 pressure_input->GetInputValue(&pressureup, gaussup);   
    5282         istemperatelayer = false;
    5283         if (enthalpyup>=matpar->PureIceEnthalpy(pressureup))
    5284           istemperatelayer=true;
    5285 
    5286                 /*Add Dirichlet constraint to this node if there is no layer of temperate ice with positive thickness*/
    5287                 if(!istemperatelayer){
     4912                watercolumn_input->GetInputValue(&watercolumn,gauss);
     4913                setspc = false;
     4914                if (enthalpy>=matpar->PureIceEnthalpy(pressure)){               
     4915                        /*check if upper node is temperate, too.
     4916                                if yes, then we have a temperate layer of positive thickness and reset the spc.
     4917                                if not, apply dirichlet BC.*/
     4918                        enthalpy_input->GetInputValue(&enthalpyup, gaussup);
     4919                        pressure_input->GetInputValue(&pressureup, gaussup);   
     4920                        setspc=((enthalpyup<matpar->PureIceEnthalpy(pressureup)) && (watercolumn>=0.))?true:false;
     4921                }
     4922
     4923                if (setspc) {
    52884924                        /*Calculate enthalpy at pressure melting point */
    52894925                        h_pmp=matpar->PureIceEnthalpy(pressure);
    5290 
    52914926                        /*Apply Dirichlet condition (dof = 0 here, since there is only one degree of freedom per node)*/
    52924927                        nodes[indices[i]]->ApplyConstraint(1,h_pmp);
    52934928                }
    5294         else {
    5295           /*remove spc*/
    5296           nodes[indices[i]]->DofInFSet(0);
    5297         }
     4929                else {
     4930                        /*remove spc*/
     4931                        nodes[indices[i]]->DofInFSet(0);
     4932                }
    52984933        }
    52994934
    53004935        /*Free ressources:*/
    53014936        xDelete<int>(indices);
    5302     xDelete<int>(indicesup);
     4937        xDelete<int>(indicesup);
     4938        delete gauss;
     4939        delete gaussup;
    53034940}
    53044941/*}}}*/
    … …  
    53064943void Penta::ComputeBasalMeltingrate(void){
    53074944        /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/
     4945        /* melting rate is positive when melting, negative when refreezing*/
    53084946
    53094947        /* Intermediaries */
    … …  
    53224960        IssmDouble  vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES];
    53234961        IssmDouble  geothermalflux[NUMVERTICES];
    5324         IssmDouble  dt,meltingrate_enthalpy;
     4962        IssmDouble  dt, yts;
     4963        IssmDouble  melting_overshoot,meltingrate_enthalpy;
     4964        IssmDouble  lambda,heating;
    53254965        Friction   *friction  = NULL;
    53264966
    … …  
    53344974        /*Fetch parameters and inputs */
    53354975        latentheat=matpar->GetLatentHeat();
    5336         rho_ice=this->matpar->GetRhoIce();
     4976        rho_ice=matpar->GetRhoIce();
    53374977        GetInputListOnVertices(&vx[0],VxEnum);
    53384978        GetInputListOnVertices(&vy[0],VyEnum);
    … …  
    53594999                checkpositivethickness=true;
    53605000
     5001                _assert_(watercolumn[iv]>=0.);
     5002
    53615003                /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
    53625004                meltingrate_enthalpy=0.;
     5005                heating=0.;
    53635006                if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv]))){
    53645007                        /*ensure that no ice is at T<Tm(p), if water layer present*/
    53655008                        enthalpy[iv]=matpar->PureIceEnthalpy(pressure[iv]);
    5366                         //meltingrate_enthalpy=0.;
    53675009                }
    53685010                else if(enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv])){
    53695011                        /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/
    5370                         meltingrate_enthalpy=0.;
    5371                         checkpositivethickness=false;
    5372                 }
    5373                 else {/*do nothing, go to next check*/}
     5012                        checkpositivethickness=false; // cold base, skip next test
     5013                }
     5014                else {/*we have a temperate base, go to next test*/}
    53745015
    53755016                if(checkpositivethickness){
    … …  
    53775018                        istemperatelayer=false;
    53785019                        if(enthalpy[iv+3]>=matpar->PureIceEnthalpy(pressure[iv+3])) istemperatelayer=true;
    5379                         if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.;
     5020                        if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
    53805021                        else{
    53815022                                enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss);
    5382                                 kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure[iv]);
     5023                                kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy[iv],enthalpy[iv+NUMVERTICES2D], pressure[iv],pressure[iv+NUMVERTICES2D]); _assert_(kappa>0.);
    53835024                                for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
    53845025                        }
    … …  
    53945035
    53955036                        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]);
    5396                         meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux[iv]))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent
     5037                        // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
     5038                        heating=(heatflux+basalfriction+geothermalflux[iv]);
     5039                        meltingrate_enthalpy=heating/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent
    53975040                }
    53985041
    53995042                /*Update water column, basal meltingrate*/
    5400                 basalmeltingrate[iv]+=meltingrate_enthalpy;
    54015043                this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    5402                 if(reCast<bool,IssmDouble>(dt))
    5403                  watercolumn[iv]+=dt*meltingrate_enthalpy;
    5404                 else
    5405                  watercolumn[iv]+=meltingrate_enthalpy;
     5044                if(reCast<bool,IssmDouble>(dt)){
     5045                        if(watercolumn[iv]+meltingrate_enthalpy*dt<0.){                         
     5046                                melting_overshoot=watercolumn[iv]+meltingrate_enthalpy*dt;
     5047                                lambda=melting_overshoot/(meltingrate_enthalpy*dt); _assert_(lambda>0); _assert_(lambda<1);
     5048                                basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy;
     5049                                watercolumn[iv]=0.;
     5050                                yts=365*24*60*60;
     5051                                enthalpy[iv]+=dt/yts*lambda*heating;
     5052                        }
     5053                        else{
     5054                                basalmeltingrate[iv]=meltingrate_enthalpy;
     5055                                watercolumn[iv]+=dt*meltingrate_enthalpy;
     5056                        }
     5057                }
     5058                else{
     5059                        basalmeltingrate[iv]=meltingrate_enthalpy;
     5060                        watercolumn[iv]+=meltingrate_enthalpy;
     5061                }         
    54065062        } 
    54075063        /*feed updated variables back into model*/
    … …  
    54175073/*FUNCTION Penta::DrainWaterfraction{{{*/
    54185074void Penta::DrainWaterfraction(void){
    5419    
     5075
    54205076    /*Intermediaries*/
    54215077        bool isenthalpy;
    54225078        IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES];
     5079        IssmDouble watercolumnbase[NUMVERTICES];
    54235080        IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
    54245081        IssmDouble latentheat, dt;
     5082        IssmDouble dw, dwc;
     5083        Penta *pentabase = NULL;
    54255084
    54265085        /*Check wether enthalpy is activated*/
    54275086        parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    54285087        if(!isenthalpy) return;       
    5429 
    5430         GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
    5431         GetInputListOnVertices(&pressure[0],PressureEnum);
     5088       
     5089        /*get basal element, needed for basal watercolumn*/
     5090        pentabase=this->GetBasalElement();
     5091       
     5092        this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
     5093        this->GetInputListOnVertices(&pressure[0],PressureEnum);
     5094        pentabase->GetInputListOnVertices(&watercolumnbase[0], WatercolumnEnum);
     5095
    54325096        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    54335097        latentheat=matpar->GetLatentHeat();
    … …  
    54355099        for(int iv=0;iv<NUMVERTICES;iv++){
    54365100                matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]);
    5437 
     5101                dw=DrainageFunctionWaterfraction(waterfraction[iv], dt);
    54385102                /*drain water fraction & update enthalpy*/
    5439                 waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt);
     5103                waterfraction[iv]-=dw;
    54405104                matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);       
    5441         }
     5105                /*add drained water to watercolumn at base*/
     5106                dwc=dw*this->IceVolume();
     5107                watercolumnbase[iv%NUMVERTICES2D]+=dwc;
     5108        }
     5109
    54425110        /*feed updated results back into model*/
    54435111        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
    54445112        this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
    5445         // this->inputs->AddInput(new PentaInput(TemperatureEnum,temperature,P1Enum));    // temperature should not change during drainage...
     5113        pentabase->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumnbase,P1Enum));
     5114
    54465115}
    54475116/*}}}*/
    … …  
    54595128                input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
    54605129        }
    5461         else if(enum_type==MaterialsRheologyZbarEnum){
     5130        else if(enum_type==DamageDbarEnum){
    54625131                if(!IsOnBed()) return;
    5463                 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum);
    5464         }
    5465 
     5132                input=(Input*)material->inputs->GetInput(DamageDEnum);
     5133        }
    54665134        else{
    54675135                input=inputs->GetInput(enum_type);
    … …  
    54825150                input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
    54835151        }
    5484         else if(enum_type==MaterialsRheologyZbarEnum){
    5485                 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum);
     5152        else if(enum_type==DamageDbarEnum){
     5153                input=(Input*)material->inputs->GetInput(DamageDEnum);
    54865154        }
    54875155        else{
    … …  
    55045172                input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
    55055173        }
    5506         else if(enum_type==MaterialsRheologyZbarEnum){
    5507                 input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum);
     5174        else if(enum_type==DamageDbarEnum){
     5175                input=(Input*)material->inputs->GetInput(DamageDEnum);
    55085176        }
    55095177        else{
    … …  
    55175185        grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
    55185186        ((ControlInput*)input)->SetGradient(grad_input);
     5187
     5188}/*}}}*/
     5189/*FUNCTION Penta::ControlToVectors{{{*/
     5190void Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){
     5191
     5192        Input* input=NULL;
     5193
     5194        if(control_enum==MaterialsRheologyBbarEnum){
     5195                input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum);
     5196        }
     5197        else if(control_enum==DamageDbarEnum){
     5198                input=(Input*)material->inputs->GetInput(DamageDEnum);
     5199        }
     5200        else{
     5201                input=inputs->GetInput(control_enum);
     5202        }
     5203        if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
     5204        if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
     5205
     5206        int         sidlist[NUMVERTICES];
     5207        int         connectivity[NUMVERTICES];
     5208        IssmPDouble values[NUMVERTICES];
     5209        IssmPDouble gradients[NUMVERTICES];
     5210        IssmDouble  value,gradient;
     5211
     5212        this->GetConnectivityList(&connectivity[0]);
     5213        this->GetVertexSidList(&sidlist[0]);
     5214
     5215        GaussPenta* gauss=new GaussPenta();
     5216        for (int iv=0;iv<NUMVERTICES;iv++){
     5217                gauss->GaussVertex(iv);
     5218
     5219                ((ControlInput*)input)->GetInputValue(&value,gauss);
     5220                ((ControlInput*)input)->GetGradientValue(&gradient,gauss);
     5221
     5222                values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
     5223                gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
     5224        }
     5225        delete gauss;
     5226
     5227        vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
     5228        vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
    55195229
    55205230}/*}}}*/
    … …  
    55515261                case MaticeEnum:
    55525262                        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    5553                         break;
    5554                 case MatdamageiceEnum:
    5555                         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    5556                         this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
     5263                        this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    55575264                        break;
    55585265                default:
    … …  
    55675274        /*Delete averaged fields*/
    55685275        this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    5569         this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum);
     5276        this->material->inputs->DeleteInput(DamageDbarEnum);
    55705277
    55715278        /*clean up and return*/
    … …  
    55865293        IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    55875294        IssmDouble eps1[3],eps2[3];
    5588         IssmDouble phi[NUMVERTICES];
    55895295        IssmDouble dphi[3][NUMVERTICES];
    55905296        GaussPenta *gauss=NULL;
    … …  
    57555461ElementVector* Penta::CreatePVectorAdjointHO(void){
    57565462
    5757 
    57585463        /*Nothing to be done if not on surface*/
    57595464        if (!IsOnSurface()) return NULL;
    … …  
    57835488        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    57845489        this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    5785         this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
     5490        this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    57865491        Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    57875492        Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    … …  
    58185523                for(resp=0;resp<num_responses;resp++){
    58195524
    5820                         weights_input->GetInputValue(&weight,gauss,resp);
     5525                        weights_input->GetInputValue(&weight,gauss,responses[resp]);
    58215526
    58225527                        switch(responses[resp]){
    … …  
    59795684        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    59805685        this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    5981         this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
     5686        this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    59825687        Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    59835688        Input* vx_input      = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
    … …  
    60145719                for(resp=0;resp<num_responses;resp++){
    60155720
    6016                         weights_input->GetInputValue(&weight,gauss,resp);
     5721                        weights_input->GetInputValue(&weight,gauss,responses[resp]);
    60175722
    60185723                        switch(responses[resp]){
    … …  
    61605865        /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
    61615866
    6162         int              i,approximation;
    6163         Tria*            tria=NULL;
     5867        int   approximation;
     5868        Tria* tria=NULL;
    61645869
    61655870        /*If on water, skip grad (=0): */
    61665871        if(NoIceInElement())return;
    6167 
     5872                                       
    61685873        /*First deal with ∂/∂alpha(KU-F)*/
    61695874        switch(control_type){
    … …  
    62145919
    62155920        /*Now deal with ∂J/∂alpha*/
    6216         int        *responses = NULL;
    6217         int         num_responses,resp;
     5921        int *responses = NULL;
     5922        int num_responses,resp;
    62185923        this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    6219         this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
     5924        this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    62205925
    62215926        for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
    … …  
    62325937                        if(IsOnBed()){
    62335938                                tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6234                                 tria->GradjDragGradient(gradient,resp,control_index);
     5939                                tria->GradjDragGradient(gradient,control_index);
    62355940                                delete tria->material; delete tria;
    62365941                        }
    … …  
    62395944                        if(IsOnBed()){
    62405945                                tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6241                                 tria->GradjBGradient(gradient,resp,control_index);
     5946                                tria->GradjBGradient(gradient,control_index);
    62425947                                delete tria->material; delete tria;
    62435948                        }
    … …  
    64326137        /*Depth Average B*/
    64336138        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
     6139        this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    64346140
    64356141        /*Collapse element to the base*/
    … …  
    64406146        /*delete Average B*/
    64416147        this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     6148        this->material->inputs->DeleteInput(DamageDbarEnum);
    64426149
    64436150} /*}}}*/
    … …  
    64486155        if (!IsOnBed()) return;
    64496156
    6450         /*Depth Average B*/
     6157        /*Depth Average B and D*/
    64516158        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
     6159        this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    64526160
    64536161        /*Collapse element to the base*/
    … …  
    64586166        /*delete Average B*/
    64596167        this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     6168        this->material->inputs->DeleteInput(DamageDbarEnum);
    64606169} /*}}}*/
    64616170/*FUNCTION Penta::GradjBbarFS {{{*/
    … …  
    64656174        if (!IsOnBed()) return;
    64666175
    6467         /*Depth Average B*/
     6176        /*Depth Average B and D*/
    64686177        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
     6178        this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    64696179
    64706180        /*Collapse element to the base*/
    … …  
    64756185        /*delete Average B*/
    64766186        this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     6187        this->material->inputs->DeleteInput(DamageDbarEnum);
    64776188} /*}}}*/
    64786189/*FUNCTION Penta::InputControlUpdate{{{*/
    … …  
    64946205                        input=(Input*)material->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
    64956206                }
    6496                 else if(control_type[i]==MaterialsRheologyZbarEnum){
     6207                else if(control_type[i]==DamageDbarEnum){
    64976208                        if (!IsOnBed()) goto cleanup_and_return;
    6498                         input=(Input*)material->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input);
     6209                        input=(Input*)material->inputs->GetInput(DamageDEnum); _assert_(input);
    64996210                }
    65006211                else{
    … …  
    65116222                        this->InputExtrude(MaterialsRheologyBEnum,MaterialsEnum);
    65126223                }
    6513                 else if(control_type[i]==MaterialsRheologyZbarEnum){
    6514                         this->InputExtrude(MaterialsRheologyZEnum,MaterialsEnum);
     6224                else if(control_type[i]==DamageDbarEnum){
     6225                        this->InputExtrude(DamageDEnum,MaterialsEnum);
    65156226                }
    65166227        }
    … …  
    65286239        int*         pdoflist=NULL;
    65296240        IssmDouble   FSreconditioning;
    6530         GaussPenta  *gauss;
    65316241
    65326242        /*Fetch number of nodes and dof for this finite element*/
    … …  
    66226332/*}}}*/
    66236333/*FUNCTION Penta::SurfaceAverageVelMisfit {{{*/
    6624 IssmDouble Penta::SurfaceAverageVelMisfit(int weight_index){
     6334IssmDouble Penta::SurfaceAverageVelMisfit(void){
    66256335
    66266336        int    approximation;
    … …  
    66456355                 * and compute SurfaceAverageVelMisfit*/
    66466356                tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6647                 J=tria->SurfaceAverageVelMisfit(weight_index);
     6357                J=tria->SurfaceAverageVelMisfit();
    66486358                delete tria->material; delete tria;
    66496359                return J;
    … …  
    66526362
    66536363                tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    6654                 J=tria->SurfaceAverageVelMisfit(weight_index);
     6364                J=tria->SurfaceAverageVelMisfit();
    66556365                delete tria->material; delete tria;
    66566366                return J;
    … …  
    66596369/*}}}*/
    66606370/*FUNCTION Penta::SurfaceAbsVelMisfit {{{*/
    6661 IssmDouble Penta::SurfaceAbsVelMisfit(int weight_index){
     6371IssmDouble Penta::SurfaceAbsVelMisfit(void){
    66626372
    66636373        int    approximation;
    … …  
    66826392                 * and compute SurfaceAbsVelMisfit*/
    66836393                tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6684                 J=tria->SurfaceAbsVelMisfit(weight_index);
     6394                J=tria->SurfaceAbsVelMisfit();
    66856395                delete tria->material; delete tria;
    66866396                return J;
    … …  
    66896399
    66906400                tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    6691                 J=tria->SurfaceAbsVelMisfit(weight_index);
     6401                J=tria->SurfaceAbsVelMisfit();
    66926402                delete tria->material; delete tria;
    66936403                return J;
    … …  
    66966406/*}}}*/
    66976407/*FUNCTION Penta::SurfaceLogVelMisfit {{{*/
    6698 IssmDouble Penta::SurfaceLogVelMisfit(int weight_index){
     6408IssmDouble Penta::SurfaceLogVelMisfit(void){
    66996409
    67006410        int    approximation;
    … …  
    67196429                 * and compute SurfaceLogVelMisfit*/
    67206430                tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6721                 J=tria->SurfaceLogVelMisfit(weight_index);
     6431                J=tria->SurfaceLogVelMisfit();
    67226432                delete tria->material; delete tria;
    67236433                return J;
    … …  
    67266436
    67276437                tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    6728                 J=tria->SurfaceLogVelMisfit(weight_index);
     6438                J=tria->SurfaceLogVelMisfit();
    67296439                delete tria->material; delete tria;
    67306440                return J;
    … …  
    67336443/*}}}*/
    67346444/*FUNCTION Penta::SurfaceLogVxVyMisfit {{{*/
    6735 IssmDouble Penta::SurfaceLogVxVyMisfit(int weight_index){
     6445IssmDouble Penta::SurfaceLogVxVyMisfit(void){
    67366446
    67376447        IssmDouble J;
    … …  
    67586468                 * and compute SurfaceLogVxVyMisfit*/
    67596469                tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6760                 J=tria->SurfaceLogVxVyMisfit(weight_index);
     6470                J=tria->SurfaceLogVxVyMisfit();
    67616471                delete tria->material; delete tria;
    67626472                return J;
    … …  
    67656475
    67666476                tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    6767                 J=tria->SurfaceLogVxVyMisfit(weight_index);
     6477                J=tria->SurfaceLogVxVyMisfit();
    67686478                delete tria->material; delete tria;
    67696479                return J;
    … …  
    67726482/*}}}*/
    67736483/*FUNCTION Penta::SurfaceRelVelMisfit {{{*/
    6774 IssmDouble Penta::SurfaceRelVelMisfit(int weight_index){
     6484IssmDouble Penta::SurfaceRelVelMisfit(void){
    67756485
    67766486        int    approximation;
    … …  
    67956505                 * and compute SurfaceRelVelMisfit*/
    67966506                tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    6797                 J=tria->SurfaceRelVelMisfit(weight_index);
     6507                J=tria->SurfaceRelVelMisfit();
    67986508                delete tria->material; delete tria;
    67996509                return J;
    … …  
    68026512
    68036513                tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    6804                 J=tria->SurfaceRelVelMisfit(weight_index);
     6514                J=tria->SurfaceRelVelMisfit();
    68056515                delete tria->material; delete tria;
    68066516                return J;
    … …  
    68096519/*}}}*/
    68106520/*FUNCTION Penta::ThicknessAbsGradient{{{*/
    6811 IssmDouble Penta::ThicknessAbsGradient(int weight_index){
     6521IssmDouble Penta::ThicknessAbsGradient(void){
    68126522
    68136523        _error_("Not implemented yet");
    … …  
    68156525/*}}}*/
    68166526/*FUNCTION Penta::ThicknessAbsMisfit {{{*/
    6817 IssmDouble Penta::ThicknessAbsMisfit(int weight_index){
     6527IssmDouble Penta::ThicknessAbsMisfit(void){
    68186528
    68196529        int    approximation;
    … …  
    68296539
    68306540        tria=(Tria*)SpawnTria(0);
    6831         J=tria->ThicknessAbsMisfit(weight_index);
     6541        J=tria->ThicknessAbsMisfit();
    68326542        delete tria->material; delete tria;
    68336543        return J;
    … …  
    68356545/*}}}*/
    68366546/*FUNCTION Penta::DragCoefficientAbsGradient{{{*/
    6837 IssmDouble Penta::DragCoefficientAbsGradient(int weight_index){
     6547IssmDouble Penta::DragCoefficientAbsGradient(void){
    68386548
    68396549        IssmDouble J;
    … …  
    68446554
    68456555        tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1
    6846         J=tria->DragCoefficientAbsGradient(weight_index);
     6556        J=tria->DragCoefficientAbsGradient();
    68476557        delete tria->material; delete tria;
    68486558        return J;
    … …  
    68506560/*}}}*/
    68516561/*FUNCTION Penta::RheologyBbarAbsGradient{{{*/
    6852 IssmDouble Penta::RheologyBbarAbsGradient(int weight_index){
     6562IssmDouble Penta::RheologyBbarAbsGradient(void){
    68536563
    68546564        IssmDouble J;
    … …  
    68596569
    68606570        tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1
    6861         J=tria->RheologyBbarAbsGradient(weight_index);
     6571        J=tria->RheologyBbarAbsGradient();
    68626572        delete tria->material; delete tria;
    68636573        return J;
    … …  
    75927302        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    75937303        delete Ke1; delete Ke2;
    7594        
     7304
    75957305        /*Compute HO Matrix with P1 element type\n");*/
    75967306        this->element_type=P1Enum;
    … …  
    77247434                case MaticeEnum:
    77257435                        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    7726                         break;
    7727                 case MatdamageiceEnum:
    7728                         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    7729                         this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
     7436                        this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    77307437                        break;
    77317438                default:
    … …  
    77407447        /*Delete averaged fields*/
    77417448        this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    7742         this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum);
     7449        this->material->inputs->DeleteInput(DamageDbarEnum);
    77437450
    77447451        /*clean up and return*/
    … …  
    79257632        int         i,j;
    79267633        IssmDouble  Jdet,viscosity;
    7927         IssmDouble  epsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
    79287634        IssmDouble  xyz_list[NUMVERTICES][3];
    79297635        IssmDouble  B[3][numdof2d];
    … …  
    80207726
    80217727        /*Intermediaries */
    8022         int         i,j;
    80237728        int         approximation;
    80247729        IssmDouble  xyz_list[NUMVERTICES][3];
    … …  
    80657770
    80667771                D_scalar=2*newviscosity*gauss->weight*Jdet;
    8067                 for (i=0;i<5;i++) D[i*5+i]=D_scalar;
     7772                for(int i=0;i<5;i++) D[i*5+i]=D_scalar;
    80687773
    80697774                TripleMultiply(B,5,numdof,1,
    … …  
    82217926
    82227927        /*Intermediaries */
    8223         int        i,j,approximation;
     7928        int        i,j;
    82247929        IssmDouble Jdet,viscosity,FSreconditioning,diameter,rigidity;
    82257930        IssmDouble xyz_list[NUMVERTICES][3];
    82267931        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    8227         IssmDouble B[8][24];
    8228         IssmDouble B_prime[8][24];
    8229         IssmDouble B_stab[3][NUMVERTICES];
    8230         IssmDouble D_scalar,D_scalar_stab;
    8231         IssmDouble D[8][8]={0.0};
    8232         IssmDouble D_stab[3][3]={0.0};
    8233         IssmDouble Ke_temp[24][24]={0.0}; //for the six nodes
    8234         IssmDouble Ke_temp_stab[6][6]={0.0}; //for the six nodes
    82357932        GaussPenta *gauss=NULL;
    82367933
    82377934        /*Stabilization*/
    8238         bool       stabilization = true;
     7935        IssmDouble D_scalar;
    82397936        IssmDouble dbasis[3][6];
    82407937        IssmDouble dmu[3];
    … …  
    82817978
    82827979                D_scalar=gauss->weight*Jdet;
    8283 
    82847980
    82857981                /*Add stabilization*/
    … …  
    83338029
    83348030        /*Intermediaries */
    8335         int        i,j,approximation;
     8031        int        i,approximation;
    83368032        IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
    83378033        IssmDouble xyz_list[NUMVERTICES][3];
    … …  
    84088104        int         analysis_type,approximation;
    84098105        IssmDouble  alpha2,Jdet2d;
    8410         IssmDouble  FSreconditioning,viscosity;
    8411         IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    84128106        IssmDouble  xyz_list[NUMVERTICES][3];
    84138107        IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
    … …  
    84358129        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    84368130        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    8437         parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    84388131        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    84398132        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    … …  
    84528145                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    84538146                GetLFS(BFriction,gauss);
    8454 
    8455                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8456                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    84578147
    84588148                friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    … …  
    85988288
    85998289        /*Intermediaries */
    8600         int         i,j;
    8601         int         approximation;
     8290        int         i,approximation;
    86028291        IssmDouble  viscosity,Jdet;
    86038292        IssmDouble  FSreconditioning;
    … …  
    86768365        IssmDouble  xyz_list[NUMVERTICES][3];
    86778366        IssmDouble  basis[6]; //for the six nodes of the penta
    8678         Tria*       tria=NULL;
    86798367        Friction*   friction=NULL;
    86808368        GaussPenta  *gauss=NULL;
    … …  
    88398527        IssmDouble  xyz_list[NUMVERTICES][3];
    88408528        IssmDouble  basis[6]; //for the six nodes of the penta
    8841         Tria*       tria=NULL;
    88428529        Friction*   friction=NULL;
    88438530        GaussPenta  *gauss=NULL;
    … …  
    91378824
    91388825        /*Intermediaries*/
    9139         int         i,j;
    91408826        IssmDouble  Jdet;
    91418827        IssmDouble  slope[3]; //do not put 2! this goes into GetInputDerivativeValue, which addresses slope[3] also!
    … …  
    91708856                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
    91718857
    9172                 for(i=0;i<numnodes;i++){
     8858                for(int i=0;i<numnodes;i++){
    91738859                        pe->values[i*NDOF2+0]+= -driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
    91748860                        pe->values[i*NDOF2+1]+= -driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
    … …  
    93829068        /*Intermediaries*/
    93839069        int        i,j;
    9384         int        approximation;
    93859070        IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab,viscosity;
    93869071        IssmDouble forcex,forcey,forcez,diameter,FSreconditioning;
    93879072        IssmDouble xyz_list[NUMVERTICES][3];
    93889073        IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    9389         IssmDouble l1l6[6]; //for the six nodes and the bubble
    9390         IssmDouble dh1dh6[3][NUMVERTICES];
    93919074        GaussPenta *gauss=NULL;
    93929075
    93939076        /*Stabilization*/
    9394         bool       stabilization = true;
    93959077        IssmDouble dbasis[3][6];
    93969078        IssmDouble dmu[3];
    … …  
    93989080        IssmDouble dnodalbasis[6][6][3];
    93999081        IssmDouble SW[6][4][4];
    9400         int p,q,ii;
     9082        int p,ii;
    94019083        int c=3; //index of pressure
    94029084
    … …  
    94329114
    94339115                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    9434                 GetNodalFunctionsP1(&l1l6[0], gauss);
    94359116                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    94369117                material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    … …  
    94739154
    94749155        /*Intermediaries*/
    9475         int        i,j;
     9156        int        i;
    94769157        int        approximation;
    94779158        IssmDouble Jdet,gravity,rho_ice;
    … …  
    97869467                case MaticeEnum:
    97879468                        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    9788                         break;
    9789                 case MatdamageiceEnum:
    9790                         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
    9791                         this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
     9469                        this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum,MaterialsEnum);
    97929470                        break;
    97939471                default:
    … …  
    98029480        /*Delete averaged inputs*/
    98039481        this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
    9804         this->material->inputs->DeleteInput(MaterialsRheologyZbarEnum);
     9482        this->material->inputs->DeleteInput(DamageDbarEnum);
    98059483
    98069484        /*clean up and return*/
    … …  
    98109488/*FUNCTION Penta::CreateJacobianStressbalanceHO{{{*/
    98119489ElementMatrix* Penta::CreateJacobianStressbalanceHO(void){
    9812 
    9813         /*Constants*/
    9814         const int    numdof=NDOF2*NUMVERTICES;
    98159490
    98169491        /*Intermediaries */
    … …  
    98769551
    98779552        /*Intermediaries */
    9878         int        i,j,approximation;
     9553        int        i,j;
    98799554        IssmDouble xyz_list[NUMVERTICES][3];
    98809555        IssmDouble Jdet;
    … …  
    100329707                values[i*NDOF2+0]=vx;
    100339708                values[i*NDOF2+1]=vy;
    10034         }
    10035 
    10036         /*Add value to global vector*/
    10037         solution->SetValues(numdof,doflist,values,INS_VAL);
    10038 
    10039         /*Free ressources:*/
    10040         delete gauss;
    10041         xDelete<int>(doflist);
    10042 }
    10043 /*}}}*/
    10044 /*FUNCTION Penta::GetSolutionFromInputsStressbalanceVert{{{*/
    10045 void  Penta::GetSolutionFromInputsStressbalanceVert(Vector<IssmDouble>* solution){
    10046 
    10047         const int    numdof=NDOF1*NUMVERTICES;
    10048 
    10049         int          i;
    10050         int*         doflist=NULL;
    10051         IssmDouble       vz;
    10052         IssmDouble       values[numdof];
    10053         GaussPenta*  gauss=NULL;
    10054 
    10055         /*Get dof list: */
    10056         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    10057         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    10058 
    10059         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    10060         /*P1 element only for now*/
    10061         gauss=new GaussPenta();
    10062         for(i=0;i<NUMVERTICES;i++){
    10063                 /*Recover vz */
    10064                 gauss->GaussVertex(i);
    10065                 vz_input->GetInputValue(&vz,gauss);
    10066                 values[i]=vz;
    100679709        }
    100689710
    … …  
    1041210054void  Penta::InputUpdateFromSolutionStressbalanceSSAFS(IssmDouble* solution){
    1041310055
    10414         const int    numdofm=NDOF2*NUMVERTICES;
    10415         const int    numdofs=NDOF3*NUMVERTICES;
    10416         const int    numdof2d=NDOF2*NUMVERTICES2D;
    10417         const int    numdofpressure=NDOF1*NUMVERTICES;
     10056        const int    numdofSSA = NDOF2*NUMVERTICES;
     10057        const int    numdof2d  = NDOF2*NUMVERTICES2D;
     10058        const int    numdofFSv = NDOF3*NUMVERTICES;
     10059        const int    numdofFSp = NDOF1*NUMVERTICES;
    1041810060
    1041910061        int     i;
    1042010062        IssmDouble  FSreconditioning;
    10421         IssmDouble  SSA_values[numdofm];
    10422         IssmDouble  FS_values[numdofs];
    10423         IssmDouble  Pressure_values[numdofs];
     10063        IssmDouble  SSA_values[numdofSSA];
     10064        IssmDouble  FS_values[numdofFSv+numdofFSp];
    1042410065        IssmDouble  vx[NUMVERTICES];
    1042510066        IssmDouble  vy[NUMVERTICES];
    … …  
    1043010071        IssmDouble  pressure[NUMVERTICES];
    1043110072        IssmDouble  xyz_list[NUMVERTICES][3];
    10432         int*    doflistm        = NULL;
    10433         int*    doflists        = NULL;
    10434         int*    doflistpressure = NULL;
    10435         Penta   *penta          = NULL;
     10073        int   *doflistSSA = NULL;
     10074        int   *doflistFSv = NULL;
     10075        int   *doflistFSp = NULL;
     10076        Penta *penta      = NULL;
     10077
     10078        /*Prepare coordinate system list*/
     10079        int* cs_list = xNew<int>(2*NUMVERTICES);
     10080        for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
     10081        for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
    1043610082
    1043710083        /*OK, we have to add results of this element for SSA
    … …  
    1044010086
    1044110087        /*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
    10442         penta->GetDofList(&doflistm,SSAApproximationEnum,GsetEnum);
    10443         GetDofList(&doflists,FSvelocityEnum,GsetEnum);
    10444         GetDofListPressure(&doflistpressure,GsetEnum);
     10088        penta->GetDofList(&doflistSSA,SSAApproximationEnum,GsetEnum);
     10089        GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
     10090        GetDofListPressure(&doflistFSp,GsetEnum);
    1044510091        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1044610092
    … …  
    1045010096        /*Use the dof list to index into the solution vector: */
    1045110097        for(i=0;i<numdof2d;i++){
    10452                 SSA_values[i]=solution[doflistm[i]];
    10453                 SSA_values[i+numdof2d]=solution[doflistm[i]];
    10454         }
    10455         for(i=0;i<numdofs;i++)FS_values[i]=solution[doflists[i]];
    10456         for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
     10098                SSA_values[i]=solution[doflistSSA[i]];
     10099                SSA_values[i+numdof2d]=solution[doflistSSA[i]];
     10100        }
     10101        for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
     10102        for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
    1045710103
    1045810104        /*Transform solution in Cartesian Space*/
    1045910105        TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
    10460         TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
     10106        TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    1046110107
    1046210108        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    1046310109        for(i=0;i<NUMVERTICES;i++){
    10464                 vx[i]=FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
    10465                 vy[i]=FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
    10466                 vzFS[i]=FS_values[i*NDOF3+2];
    10467                 pressure[i]=Pressure_values[i*NDOF1]*FSreconditioning;
     10110                vx[i]       = FS_values[i*NDOF3+0]+SSA_values[i*NDOF2+0];
     10111                vy[i]       = FS_values[i*NDOF3+1]+SSA_values[i*NDOF2+1];
     10112                vzFS[i]     = FS_values[i*NDOF3+2];
     10113                pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
    1046810114
    1046910115                /*Check solution*/
    1047010116                if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    1047110117                if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    10472                 if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
     10118                if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
    1047310119                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    1047410120        }
    … …  
    1048810134        /*Now Compute vel*/
    1048910135        for(i=0;i<NUMVERTICES;i++) {
    10490                 vz[i]=vzSSA[i]+vzFS[i];
    10491                 vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     10136                vz[i]  = vzSSA[i]+vzFS[i];
     10137                vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1049210138        }
    1049310139
    … …  
    1050810154
    1050910155        /*Free ressources:*/
    10510         xDelete<int>(doflistm);
    10511         xDelete<int>(doflists);
    10512         xDelete<int>(doflistpressure);
     10156        xDelete<int>(doflistSSA);
     10157        xDelete<int>(doflistFSv);
     10158        xDelete<int>(doflistFSp);
     10159        xDelete<int>(cs_list);
    1051310160}
    1051410161/*}}}*/
    … …  
    1067010317void  Penta::InputUpdateFromSolutionStressbalanceHOFS(IssmDouble* solution){
    1067110318
    10672         const int    numdofp=NDOF2*NUMVERTICES;
    10673         const int    numdofs=NDOF3*NUMVERTICES;
    10674         const int    numdofpressure=NDOF1*NUMVERTICES;
     10319        const int    numdofHO  = NDOF2*NUMVERTICES;
     10320        const int    numdofFSv = NDOF3*NUMVERTICES;
     10321        const int    numdofFSp = NDOF1*NUMVERTICES;
    1067510322
    1067610323        int        i;
    10677         IssmDouble HO_values[numdofp];
    10678         IssmDouble FS_values[numdofs];
    10679         IssmDouble Pressure_values[numdofpressure];
     10324        IssmDouble HO_values[numdofHO];
     10325        IssmDouble FS_values[numdofFSv+numdofFSp];
    1068010326        IssmDouble vx[NUMVERTICES];
    1068110327        IssmDouble vy[NUMVERTICES];
    … …  
    1068710333        IssmDouble xyz_list[NUMVERTICES][3];
    1068810334        IssmDouble FSreconditioning;
    10689         int*       doflistp        = NULL;
    10690         int*       doflists        = NULL;
    10691         int*       doflistpressure = NULL;
    10692         Penta      *penta          = NULL;
    10693 
    10694         /*OK, we have to add results of this element for HO
    10695          * and results from the penta at base for SSA. Now recover results*/
    10696         penta=GetBasalElement();
     10335        int*       doflistHO        = NULL;
     10336        int*       doflistFSv        = NULL;
     10337        int*       doflistFSp = NULL;
     10338
     10339        /*Prepare coordinate system list*/
     10340        int* cs_list = xNew<int>(2*NUMVERTICES);
     10341        for(i=0;i<NUMVERTICES;i++) cs_list[i]             = XYZEnum;
     10342        for(i=0;i<NUMVERTICES;i++) cs_list[NUMVERTICES+i] = PressureEnum;
    1069710343
    1069810344        /*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
    10699         GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
    10700         GetDofList(&doflists,FSvelocityEnum,GsetEnum);
    10701         GetDofListPressure(&doflistpressure,GsetEnum);
     10345        GetDofList(&doflistHO,HOApproximationEnum,GsetEnum);
     10346        GetDofList(&doflistFSv,FSvelocityEnum,GsetEnum);
     10347        GetDofListPressure(&doflistFSp,GsetEnum);
    1070210348        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1070310349
    … …  
    1070610352
    1070710353        /*Use the dof list to index into the solution vector: */
    10708         for(i=0;i<numdofp;i++) HO_values[i]=solution[doflistp[i]];
    10709         for(i=0;i<numdofs;i++) FS_values[i]=solution[doflists[i]];
    10710         for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
     10354        for(i=0;i<numdofHO;i++)  HO_values[i]=solution[doflistHO[i]];
     10355        for(i=0;i<numdofFSv;i++) FS_values[i]=solution[doflistFSv[i]];
     10356        for(i=0;i<numdofFSp;i++) FS_values[numdofFSv+i]=solution[doflistFSp[i]];
    1071110357
    1071210358        /*Transform solution in Cartesian Space*/
    1071310359        TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
    10714         TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZEnum);
     10360        TransformSolutionCoord(&FS_values[0],this->nodes,2*NUMVERTICES,cs_list);
    1071510361
    1071610362        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    1071710363        for(i=0;i<NUMVERTICES;i++){
    10718                 vx[i]=FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
    10719                 vy[i]=FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
    10720                 vzFS[i]=FS_values[i*NDOF3+2];
    10721                 pressure[i]=Pressure_values[i*NDOF1]*FSreconditioning;
     10364                vx[i]       = FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
     10365                vy[i]       = FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
     10366                vzFS[i]     = FS_values[i*NDOF3+2];
     10367                pressure[i] = FS_values[NUMVERTICES*NDOF3+i]*FSreconditioning;
    1072210368
    1072310369                /*Check solution*/
    1072410370                if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
    1072510371                if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
    10726                 if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
     10372                if(xIsNan<IssmDouble>(vzFS[i]))     _error_("NaN found in solution vector");
    1072710373                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    1072810374        }
    … …  
    1074210388        /*Now Compute vel*/
    1074310389        for(i=0;i<NUMVERTICES;i++) {
    10744                 vz[i]=vzHO[i]+vzFS[i];
    10745                 vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     10390                vz[i]  = vzHO[i]+vzFS[i];
     10391                vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1074610392        }
    1074710393
    … …  
    1076210408
    1076310409        /*Free ressources:*/
    10764         xDelete<int>(doflistp);
    10765         xDelete<int>(doflists);
    10766         xDelete<int>(doflistpressure);
     10410        xDelete<int>(doflistHO);
     10411        xDelete<int>(doflistFSv);
     10412        xDelete<int>(doflistFSp);
     10413        xDelete<int>(cs_list);
    1076710414}
    1076810415/*}}}*/
    … …  
    1094210589        int*         pdoflist=NULL;
    1094310590        IssmDouble   FSreconditioning;
    10944         GaussPenta  *gauss;
    1094510591
    1094610592        /*Fetch number of nodes and dof for this finite element*/
    … …  
    1095110597
    1095210598        /*Initialize values*/
    10953         IssmDouble* vvalues  = xNew<IssmDouble>(vnumdof);
    10954         IssmDouble* pvalues  = xNew<IssmDouble>(pnumdof);
     10599        IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
    1095510600        IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
    1095610601        IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
    1095710602        IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
    1095810603        IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
    10959         IssmDouble* pressure = xNew<IssmDouble>(vnumnodes);
     10604        IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
     10605
     10606        /*Prepare coordinate system list*/
     10607        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     10608        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     10609        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    1096010610
    1096110611        /*Get dof list: */
    … …  
    1096410614
    1096510615        /*Use the dof list to index into the solution vector: */
    10966         for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
    10967         for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
     10616        for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
     10617        for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
    1096810618
    1096910619        /*Transform solution in Cartesian Space*/
    10970         TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYZEnum);
     10620        TransformSolutionCoord(&values[0],nodes,vnumnodes+pnumdof,cs_list);
    1097110621
    1097210622        /*Ok, we have vx and vy in values, fill in all arrays: */
    1097310623        for(i=0;i<vnumnodes;i++){
    10974                 vx[i] = vvalues[i*NDOF3+0];
    10975                 vy[i] = vvalues[i*NDOF3+1];
    10976                 vz[i] = vvalues[i*NDOF3+2];
     10624                vx[i] = values[i*NDOF3+0];
     10625                vy[i] = values[i*NDOF3+1];
     10626                vz[i] = values[i*NDOF3+2];
    1097710627                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    1097810628                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    … …  
    1098010630        }
    1098110631        for(i=0;i<pnumnodes;i++){
    10982                 pressure[i] = pvalues[i];
     10632                pressure[i] = values[vnumdof+i];
    1098310633                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    1098410634        }
    … …  
    1098610636        /*Recondition pressure and compute vel: */
    1098710637        this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    10988         for(i=0;i<pnumnodes;i++) pressure[i]=pressure[i]*FSreconditioning;
    10989         for(i=0;i<vnumnodes;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     10638        for(i = 0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
     10639        for(i = 0;i<vnumnodes;i++) vel[i]      = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1099010640
    1099110641        /*Now, we have to move the previous inputs  to old
    … …  
    1100910659        xDelete<IssmDouble>(vy);
    1101010660        xDelete<IssmDouble>(vx);
    11011         xDelete<IssmDouble>(vvalues);
    11012         xDelete<IssmDouble>(pvalues);
     10661        xDelete<IssmDouble>(values);
    1101310662        xDelete<int>(vdoflist);
    1101410663        xDelete<int>(pdoflist);
     10664        xDelete<int>(cs_list);
    1101510665}
    1101610666/*}}}*/
    … …  
    1127810928        IssmDouble h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],r[NUMVERTICES];
    1127910929        IssmDouble melting[NUMVERTICES],phi[NUMVERTICES];
    11280         bool       grounded[NUMVERTICES],floating[NUMVERTICES];
    1128110930
    1128210931        if(!IsOnBed()) return;
    … …  
    1130310952                                b[i]        = r[i];
    1130410953                                s[i]        = b[i]+h[i];
    11305                                 floating[i] = false;
    11306                                 grounded[i] = true;
    1130710954                        }
    1130810955                }
    … …  
    1131610963                                        s[i]        = (1-density)*h[i];
    1131710964                                        b[i]        = -density*h[i];
    11318                                         floating[i] = true;
    11319                                         grounded[i] = false;
    1132010965                                }
    1132110966                                else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){
    1132210967                                        s[i]        = (1-density)*h[i];
    1132310968                                        b[i]        = -density*h[i];
    11324                                         floating[i] = true;
    11325                                         grounded[i] = false;
    1132610969                                }
    1132710970                                else{
Note: See TracChangeset for help on using the changeset viewer.