source: issm/trunk-jpl/src/c/classes/Loads/Icefront.cpp@ 15298

/*!\file Icefront.c
 * \brief: implementation of the Icefront object
 */

/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../classes.h"
#include "shared/shared.h"
/*}}}*/ 

/*Load macros*/
#define NUMVERTICESSEG 2
#define NUMVERTICESQUA 4

/*Icefront constructors and destructor*/
/*FUNCTION Icefront::Icefront() {{{*/
Icefront::Icefront(){

        this->inputs=NULL;
        this->parameters=NULL;

        this->hnodes=NULL;
        this->nodes= NULL;
        this->hvertices=NULL;
        this->vertices= NULL;
        this->helement=NULL;
        this->element= NULL;
        this->hmatpar=NULL;
        this->matpar= NULL;
}
/*}}}*/
/*FUNCTION Icefront::Icefront(int id, int i, IoModel* iomodel,int analysis_type) {{{*/
Icefront::Icefront(int icefront_id,int i, IoModel* iomodel,int in_icefront_type, int in_analysis_type){

        int segment_width;
        int element;
        int numnodes; 
        int numvertices; 
        int dim;
        int numberofelements;

        /*icefront constructor data: */
        int  icefront_eid;
        int  icefront_mparid;
        int  icefront_node_ids[NUMVERTICESQUA]; //initialize with largest size
        int  icefront_vertex_ids[NUMVERTICESQUA]; //initialize with largest size
        int  icefront_fill;

        /*find parameters: */
        iomodel->Constant(&dim,MeshDimensionEnum);
        iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);

        /*First, retrieve element index and element type: */
        if (dim==2){
                segment_width=4;
        }
        else{
                segment_width=6;
        }
        _assert_(iomodel->Data(DiagnosticIcefrontEnum));
        element=reCast<int,IssmDouble>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)-1); //element is in the penultimate column (node1 node2 ... elem fill)

        /*Build ids for hook constructors: */
        icefront_eid=reCast<int,IssmDouble>( *(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+segment_width-2)); //matlab indexing
        icefront_mparid=numberofelements+1; //matlab indexing

        if (in_icefront_type==MacAyeal2dIceFrontEnum || in_icefront_type==MacAyeal3dIceFrontEnum){
                icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
                icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
                icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
                icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
        }
        else if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
                icefront_node_ids[0]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
                icefront_node_ids[1]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
                icefront_node_ids[2]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
                icefront_node_ids[3]=iomodel->nodecounter+reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
                icefront_vertex_ids[0]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+0));
                icefront_vertex_ids[1]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+1));
                icefront_vertex_ids[2]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+2));
                icefront_vertex_ids[3]=reCast<int>(*(iomodel->Data(DiagnosticIcefrontEnum)+segment_width*i+3));
        }
        else _error_("in_icefront_type " << EnumToStringx(in_icefront_type) << " not supported yet!");

        if (in_icefront_type==PattynIceFrontEnum || in_icefront_type==StokesIceFrontEnum){
                numnodes=4;
                numvertices=4;
        }
        else{
                numnodes=2;
                numvertices=2;
        }

        /*Fill*/
        icefront_fill=reCast<int>(iomodel->Data(DiagnosticIcefrontEnum)[segment_width*i+segment_width-1]);

        /*Ok, we have everything to build the object: */
        this->id=icefront_id;
        this->analysis_type=in_analysis_type;

        /*Hooks: */
        this->hnodes=new Hook(icefront_node_ids,numnodes);
        this->hvertices=new Hook(icefront_vertex_ids,numvertices);
        this->helement=new Hook(&icefront_eid,1);
        this->hmatpar=new Hook(&icefront_mparid,1);

        //intialize  and add as many inputs per element as requested: 
        this->inputs=new Inputs();
        this->inputs->AddInput(new IntInput(FillEnum,icefront_fill));
        this->inputs->AddInput(new IntInput(IceFrontTypeEnum,in_icefront_type));

        //parameters and hooked fields: we still can't point to them, they may not even exist. Configure will handle this.
        this->parameters = NULL;
        this->nodes      = NULL;
        this->vertices   = NULL;
        this->element    = NULL;
        this->matpar     = NULL;
}

/*}}}*/
/*FUNCTION Icefront::~Icefront() {{{*/
Icefront::~Icefront(){
        delete inputs;
        this->parameters=NULL;
        delete hnodes;
        delete hvertices;
        delete helement;
        delete hmatpar;
}
/*}}}*/

/*Object virtual functions definitions:*/
/*FUNCTION Icefront::Echo {{{*/
void Icefront::Echo(void){
        _printf_("Icefront:\n");
        _printf_("   id: " << id << "\n");
        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
        hnodes->Echo();
        hvertices->Echo();
        helement->Echo();
        hmatpar->Echo();
        _printf_("   parameters: " << parameters << "\n");
        if(parameters)parameters->Echo();
        _printf_("   inputs: " << inputs << "\n");
        if(inputs)inputs->Echo();
}
/*}}}*/
/*FUNCTION Icefront::DeepEcho{{{*/
void Icefront::DeepEcho(void){

        _printf_("Icefront:\n");
        _printf_("   id: " << id << "\n");
        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
        hnodes->DeepEcho();
        hvertices->DeepEcho();
        helement->DeepEcho();
        hmatpar->DeepEcho();
        _printf_("   parameters: " << parameters << "\n");
        if(parameters)parameters->DeepEcho();
        _printf_("   inputs: " << inputs << "\n");
        if(inputs)inputs->DeepEcho();
}
/*}}}*/
/*FUNCTION Icefront::Id {{{*/
int    Icefront::Id(void){ return id; }
/*}}}*/
/*FUNCTION Icefront::ObjectEnum{{{*/
int Icefront::ObjectEnum(void){

        return IcefrontEnum;

}
/*}}}*/
/*FUNCTION Icefront::copy {{{*/
Object* Icefront::copy() {

        Icefront* icefront=NULL;

        icefront=new Icefront();

        /*copy fields: */
        icefront->id=this->id;
        icefront->analysis_type=this->analysis_type;
        if(this->inputs){
                icefront->inputs=(Inputs*)this->inputs->Copy();
        }
        else{
                icefront->inputs=new Inputs();
        }
        /*point parameters: */
        icefront->parameters=this->parameters;

        /*now deal with hooks and objects: */
        icefront->hnodes    = (Hook*)this->hnodes->copy();
        icefront->hvertices = (Hook*)this->hvertices->copy();
        icefront->helement  = (Hook*)this->helement->copy();
        icefront->hmatpar   = (Hook*)this->hmatpar->copy();

        /*corresponding fields*/
        icefront->nodes    = (Node**)icefront->hnodes->deliverp();
        icefront->vertices = (Vertex**)icefront->hvertices->deliverp();
        icefront->element  = (Element*)icefront->helement->delivers();
        icefront->matpar   = (Matpar*)icefront->hmatpar->delivers();

        return icefront;

}
/*}}}*/

/*Load virtual functions definitions:*/
/*FUNCTION Icefront::Configure {{{*/
void  Icefront::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){

        /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective 
         * datasets, using internal ids and offsets hidden in hooks: */
        hnodes->configure((DataSet*)nodesin);
        hvertices->configure((DataSet*)verticesin);
        helement->configure((DataSet*)elementsin);
        hmatpar->configure((DataSet*)materialsin);

        /*Initialize hooked fields*/
        this->nodes    = (Node**)hnodes->deliverp();
        this->vertices = (Vertex**)hvertices->deliverp();
        this->element  = (Element*)helement->delivers();
        this->matpar   = (Matpar*)hmatpar->delivers();

        /*point parameters to real dataset: */
        this->parameters=parametersin;
}
/*}}}*/
/*FUNCTION Icefront::SetCurrentConfiguration {{{*/
void  Icefront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
}
/*}}}*/
/*FUNCTION Icefront::CreateKMatrix {{{*/
void  Icefront::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){

        /*No stiffness loads applied, do nothing: */
        return;

}
/*}}}*/
/*FUNCTION Icefront::CreatePVector {{{*/
void  Icefront::CreatePVector(Vector<IssmDouble>* pf){

        /*Checks in debugging mode*/
        _assert_(nodes);
        _assert_(element);
        _assert_(matpar);

        /*Retrieve parameters: */
        ElementVector* pe=NULL;
        int analysis_type;
        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);

        /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */
        switch(analysis_type){
                #ifdef _HAVE_DIAGNOSTIC_
                case DiagnosticHorizAnalysisEnum:
                        pe=CreatePVectorDiagnosticHoriz();
                        break;
                #endif
                #ifdef _HAVE_CONTROL_
                case AdjointHorizAnalysisEnum:
                        pe=CreatePVectorAdjointHoriz();
                        break;
                #endif
                default:
                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
        }

        /*Add to global Vector*/
        if(pe){
                pe->AddToGlobal(pf);
                delete pe;
        }
}
/*}}}*/
/*FUNCTION Icefront::CreateJacobianMatrix{{{*/
void  Icefront::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){
        this->CreateKMatrix(Jff,NULL);
}
/*}}}*/
/*FUNCTION Icefront::GetNodesSidList{{{*/
void Icefront::GetNodesSidList(int* sidlist){

        int type;
        inputs->GetInputValue(&type,IceFrontTypeEnum);
        _assert_(sidlist);
        _assert_(nodes);

        switch(type){
                case MacAyeal2dIceFrontEnum:
                case MacAyeal3dIceFrontEnum:
                        for(int i=0;i<NUMVERTICESSEG;i++) sidlist[i]=nodes[i]->Sid();
                        return;
#ifdef _HAVE_3D_
                case PattynIceFrontEnum:
                case StokesIceFrontEnum:
                        for(int i=0;i<NUMVERTICESQUA;i++) sidlist[i]=nodes[i]->Sid();
                        return;
#endif
                default:
                        _error_("Icefront type " << EnumToStringx(type) << " not supported yet");
        }
}
/*}}}*/
/*FUNCTION Icefront::GetNumberOfNodes{{{*/
int Icefront::GetNumberOfNodes(void){

        int type;
        inputs->GetInputValue(&type,IceFrontTypeEnum);

        switch(type){
                case MacAyeal2dIceFrontEnum:
                        return NUMVERTICESSEG;
#ifdef _HAVE_3D_
                case MacAyeal3dIceFrontEnum:
                        return NUMVERTICESSEG;
                case PattynIceFrontEnum:
                        return NUMVERTICESQUA;
                case StokesIceFrontEnum:
                        return NUMVERTICESQUA;
#endif
                default:
                        _error_("Icefront type " << EnumToStringx(type) << " not supported yet");
        }

}
/*}}}*/
/*FUNCTION Icefront::IsPenalty{{{*/
bool Icefront::IsPenalty(void){
        return false;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreateKMatrix {{{*/
void  Icefront::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs, IssmDouble kmax){
        /*do nothing: */
        return;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreatePVector{{{*/
void  Icefront::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){
        /*do nothing: */
        return;
}
/*}}}*/
/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{*/
void  Icefront::PenaltyCreateJacobianMatrix(Matrix<IssmDouble>* Jff,IssmDouble kmax){
        this->PenaltyCreateKMatrix(Jff,NULL,kmax);
}
/*}}}*/
/*FUNCTION Icefront::SetwiseNodeConnectivity{{{*/
void Icefront::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int set1_enum,int set2_enum){

        /*Output */
        int d_nz = 0;
        int o_nz = 0;

        /*Loop over all nodes*/
        for(int i=0;i<this->GetNumberOfNodes();i++){

                if(!flags[this->nodes[i]->Sid()]){

                        /*flag current node so that no other element processes it*/
                        flags[this->nodes[i]->Sid()]=true;

                        /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
                        switch(set2_enum){
                                case FsetEnum:
                                        if(nodes[i]->indexing.fsize){
                                                if(this->nodes[i]->IsClone())
                                                 o_nz += 1;
                                                else
                                                 d_nz += 1;
                                        }
                                        break;
                                case GsetEnum:
                                        if(nodes[i]->indexing.gsize){
                                                if(this->nodes[i]->IsClone())
                                                 o_nz += 1;
                                                else
                                                 d_nz += 1;
                                        }
                                        break;
                                case SsetEnum:
                                        if(nodes[i]->indexing.ssize){
                                                if(this->nodes[i]->IsClone())
                                                 o_nz += 1;
                                                else
                                                 d_nz += 1;
                                        }
                                        break;
                                default: _error_("not supported");
                        }
                }
        }

        /*Assign output pointers: */
        *pd_nz=d_nz;
        *po_nz=o_nz;
}
/*}}}*/
/*FUNCTION Icefront::InAnalysis{{{*/
bool Icefront::InAnalysis(int in_analysis_type){
        if (in_analysis_type==this->analysis_type)return true;
        else return false;
}
/*}}}*/

/*Update virtual functions definitions:*/
/*FUNCTION Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVector(IssmDouble* vector, int name, int type){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVector(int* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVector(int* vector, int name, int type){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVector(bool* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVector(bool* vector, int name, int type){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type) {{{*/
void  Icefront::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVectorDakota(int* vector, int name, int type){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type) {{{*/
void  Icefront::InputUpdateFromVectorDakota(bool* vector, int name, int type){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(IssmDouble constant, int name) {{{*/
void  Icefront::InputUpdateFromConstant(IssmDouble constant, int name){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(int constant, int name) {{{*/
void  Icefront::InputUpdateFromConstant(int constant, int name){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromConstant(bool constant, int name) {{{*/
void  Icefront::InputUpdateFromConstant(bool constant, int name){
        /*Nothing updated yet*/
}
/*}}}*/
/*FUNCTION Icefront::InputUpdateFromSolution{{{*/
void  Icefront::InputUpdateFromSolution(IssmDouble* solution){
        /*Nothing updated yet*/
}
/*}}}*/

/*Icefront numerics: */
#ifdef _HAVE_DIAGNOSTIC_
/*FUNCTION Icefront::CreatePVectorDiagnosticHoriz {{{*/
ElementVector* Icefront::CreatePVectorDiagnosticHoriz(void){

        int type;
        inputs->GetInputValue(&type,IceFrontTypeEnum);

        switch(type){
                case MacAyeal2dIceFrontEnum:
                        return CreatePVectorDiagnosticMacAyeal2d();
        #ifdef _HAVE_3D_
                case MacAyeal3dIceFrontEnum:
                        return CreatePVectorDiagnosticMacAyeal3d();
                case PattynIceFrontEnum:
                        return CreatePVectorDiagnosticPattyn();
                case StokesIceFrontEnum:
                        return CreatePVectorDiagnosticStokes();
            #endif
                default:
                        _error_("Icefront type " << EnumToStringx(type) << " not supported yet");
        }
}
/*}}}*/
/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal2d(void){

        /*Intermediary*/
        int         ig,index1,index2,fill;
        IssmDouble  Jdet;
        IssmDouble  thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity;
        IssmDouble  water_pressure,air_pressure,surface_under_water,base_under_water;
        IssmDouble  xyz_list[NUMVERTICESSEG][3];
        IssmDouble  normal[2];
        GaussTria *gauss;

        /*return of element is on water*/
        Tria* tria=((Tria*)element);
        if(tria->IsOnWater()) return NULL;

        /*Fetch number of nodes and dof for this finite element*/
        //int numnodes = this->NumberofNodes();
        int numnodes = 2; 
        int numdof   = numnodes*NDOF2;

        /*Initialize Element vector and vectors*/
        ElementVector* pe=new ElementVector(nodes,NUMVERTICESSEG,this->parameters,MacAyealApproximationEnum);
        IssmDouble*    basis = xNew<IssmDouble>(numnodes);

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESSEG);
        Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
        Input* bed_input      =tria->inputs->GetInput(BedEnum);       _assert_(bed_input);
        inputs->GetInputValue(&fill,FillEnum);
        rho_water=matpar->GetRhoWater();
        rho_ice  =matpar->GetRhoIce();
        gravity  =matpar->GetG();
        GetSegmentNormal(&normal[0],xyz_list);

        /*Start looping on Gaussian points*/
        index1=tria->GetNodeIndex(nodes[0]);
        index2=tria->GetNodeIndex(nodes[1]);
        gauss=new GaussTria(index1,index2,3);

        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                thickness_input->GetInputValue(&thickness,gauss);
                bed_input->GetInputValue(&bed,gauss);

                switch(fill){
                        case WaterEnum:
                                surface_under_water=min(0.,thickness+bed); // 0 if the top of the glacier is above water level
                                base_under_water=min(0.,bed);              // 0 if the bottom of the glacier is above water level
                                water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
                                break;
                        case AirEnum:
                                water_pressure=0;
                                break;
                        case IceEnum:
                                water_pressure=-1.0/2.0*gravity*rho_ice*pow(thickness,2); // we are facing a wall of ice. use water_pressure to cancel the lithostatic pressure.
                                break;
                        default:
                                _error_("fill type " << EnumToStringx(fill) << " not supported yet");
                }
                ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
                air_pressure=0;
                pressure = ice_pressure + water_pressure + air_pressure;

                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
                tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2);

                for (int i=0;i<numnodes;i++){
                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
                }
        }

        /*Transform load vector*/
        TransformLoadVectorCoord(pe,nodes,NUMVERTICESSEG,XYEnum);

        /*Clean up and return*/
        xDelete<IssmDouble>(basis);
        delete gauss;
        return pe;
}
/*}}}*/
#endif

#ifdef _HAVE_CONTROL_
/*FUNCTION Icefront::CreatePVectorAdjointHoriz {{{*/
ElementVector* Icefront::CreatePVectorAdjointHoriz(void){

        /*No load vector applied to the adjoint*/
        return NULL;
}
/*}}}*/
#endif
#ifdef _HAVE_3D_
/*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal3d{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticMacAyeal3d(void){

        Icefront *icefront = NULL;
        Penta    *penta    = NULL;
        Tria     *tria     = NULL;

        /*Cast element onto Penta*/
        penta   =(Penta*)this->element;

        /*Return if not on bed*/
        if(!penta->IsOnBed() || penta->IsOnWater()) return NULL;

        /*Spawn Tria and call MacAyeal2d*/
        tria    =(Tria*)penta->SpawnTria(0,1,2);
        icefront=(Icefront*)this->copy();
        icefront->element=tria;
        icefront->inputs->AddInput(new IntInput(IceFrontTypeEnum,MacAyeal2dIceFrontEnum));
        ElementVector* pe=icefront->CreatePVectorDiagnosticMacAyeal2d();

        /*clean-up and return*/
        delete tria->material;
        delete tria;
        delete icefront;
        return pe;
}
/*}}}*/
/*FUNCTION Icefront::CreatePVectorDiagnosticPattyn{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticPattyn(void){

        /*Constants*/
        const int numdofs = NUMVERTICESQUA *NDOF2;

        /*Intermediaries*/
        int         i,j,ig,index1,index2,index3,index4;
        int         fill;
        IssmDouble      surface,pressure,ice_pressure,rho_water,rho_ice,gravity;
        IssmDouble      water_pressure,air_pressure;
        IssmDouble      Jdet,z_g;
        IssmDouble      xyz_list[NUMVERTICESQUA][3];
        IssmDouble      normal[3];
        IssmDouble      l1l4[4];
        GaussPenta *gauss = NULL;

        Penta* penta=(Penta*)element;

        /*Initialize Element vector and return if necessary*/
        if(penta->IsOnWater()) return NULL;
        ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,PattynApproximationEnum);

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
        Input* surface_input  =penta->inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
        inputs->GetInputValue(&fill,FillEnum);
        rho_water=matpar->GetRhoWater();
        rho_ice  =matpar->GetRhoIce();
        gravity  =matpar->GetG();
        GetQuadNormal(&normal[0],xyz_list);

        /*Identify which nodes are in the quad: */
        index1=element->GetNodeIndex(nodes[0]);
        index2=element->GetNodeIndex(nodes[1]);
        index3=element->GetNodeIndex(nodes[2]);
        index4=element->GetNodeIndex(nodes[3]);

        /* Start  looping on the number of gaussian points: */
        IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
        IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
        if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,10); //refined in vertical because of the sea level discontinuity
        else                 gauss=new GaussPenta(index1,index2,index3,index4,3,3);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
                penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
                z_g=penta->GetZcoord(gauss);
                surface_input->GetInputValue(&surface,gauss);

                switch(fill){
                        case WaterEnum:
                                water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
                                break;
                        case AirEnum:
                                water_pressure=0;
                                break;
                        default:
                                _error_("fill type " << EnumToStringx(fill) << " not supported yet");
                }
                ice_pressure=rho_ice*gravity*(surface-z_g);
                air_pressure=0;
                pressure = ice_pressure + water_pressure + air_pressure;

                for(i=0;i<NUMVERTICESQUA;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
        }

        /*Transform load vector*/
        TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYEnum);

        /*Clean up and return*/
        delete gauss;
        return pe;
}
/*}}}*/
/*FUNCTION Icefront::CreatePVectorDiagnosticStokes{{{*/
ElementVector* Icefront::CreatePVectorDiagnosticStokes(void){

        /*Constants*/
        const int numdofs = NUMVERTICESQUA *NDOF4;

        /*Intermediaries*/
        int         i,j,ig,index1,index2,index3,index4;
        int         fill;
        IssmDouble      pressure,rho_water,gravity;
        IssmDouble      water_pressure,air_pressure;
        IssmDouble      Jdet,z_g;
        IssmDouble      xyz_list[NUMVERTICESQUA][3];
        IssmDouble      normal[3];
        IssmDouble      l1l4[4];
        GaussPenta *gauss = NULL;

        Penta* penta=(Penta*)element;

        /*Initialize Element vector and return if necessary*/
        if(penta->IsOnWater()) return NULL;
        ElementVector* pe=new ElementVector(nodes,NUMVERTICESQUA,this->parameters,StokesApproximationEnum);

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICESQUA);
        inputs->GetInputValue(&fill,FillEnum);
        rho_water=matpar->GetRhoWater();
        gravity  =matpar->GetG();
        GetQuadNormal(&normal[0],xyz_list);

        /*Identify which nodes are in the quad: */
        index1=element->GetNodeIndex(nodes[0]);
        index2=element->GetNodeIndex(nodes[1]);
        index3=element->GetNodeIndex(nodes[2]);
        index4=element->GetNodeIndex(nodes[3]);

        /* Start  looping on the number of gaussian points: */
        IssmDouble zmax=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
        IssmDouble zmin=xyz_list[0][2]; for(i=1;i<NUMVERTICESQUA;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
        if(zmax>0 && zmin<0) gauss=new GaussPenta(index1,index2,index3,index4,3,30); //refined in vertical because of the sea level discontinuity
        else                 gauss=new GaussPenta(index1,index2,index3,index4,3,3);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                penta->GetQuadNodalFunctions(l1l4,gauss,index1,index2,index3,index4);
                penta->GetQuadJacobianDeterminant(&Jdet,xyz_list,gauss);
                z_g=penta->GetZcoord(gauss);

                switch(fill){
                        case WaterEnum:
                                water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
                                break;
                        case AirEnum:
                                water_pressure=0;
                                break;
                        default:
                                _error_("fill type " << EnumToStringx(fill) << " not supported yet");
                }
                air_pressure=0;
                pressure = water_pressure + air_pressure; //no ice pressure fore Stokes

                for(i=0;i<NUMVERTICESQUA;i++){
                        for(j=0;j<NDOF4;j++){
                                if(j<3)  pe->values[i*NDOF4+j]+=Jdet*gauss->weight*pressure*l1l4[i]*normal[j];
                                else     pe->values[i*NDOF4+j]+=0; //pressure term
                        }
                }
        }

        /*Transform load vector*/
        TransformLoadVectorCoord(pe,nodes,NUMVERTICESQUA,XYZPEnum);

        /*Clean up and return*/
        delete gauss;
        return pe;
}
/*}}}*/
#endif
/*FUNCTION Icefront::GetDofList {{{*/
void  Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){

        int numberofdofs=0;
        int count=0;
        int type;
        int numberofnodes=2;

        /*output: */
        int* doflist=NULL;

        /*recover type: */
        inputs->GetInputValue(&type,IceFrontTypeEnum);

        /*Some checks for debugging*/
        _assert_(nodes);

        /*How many nodes? :*/
        if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum)
         numberofnodes=2;
        else 
         numberofnodes=4;

        /*Figure out size of doflist: */
        for(int i=0;i<numberofnodes;i++){
                numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
        }

        /*Allocate: */
        doflist=xNew<int>(numberofdofs);

        /*Populate: */
        count=0;
        for(int i=0;i<numberofnodes;i++){
                nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
                count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
        }

        /*Assign output pointers:*/
        *pdoflist=doflist;
}
/*}}}*/
/*FUNCTION Icefront::GetSegmentNormal {{{*/
void Icefront:: GetSegmentNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){

        /*Build unit outward pointing vector*/
        const int numnodes=NUMVERTICESSEG;
        IssmDouble vector[2];
        IssmDouble norm;

        vector[0]=xyz_list[1][0] - xyz_list[0][0];
        vector[1]=xyz_list[1][1] - xyz_list[0][1];

        norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));

        normal[0]= + vector[1]/norm;
        normal[1]= - vector[0]/norm;
}
/*}}}*/
/*FUNCTION Icefront::GetQuadNormal {{{*/
void Icefront:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){

        /*Build unit outward pointing vector*/
        IssmDouble AB[3];
        IssmDouble AC[3];
        IssmDouble norm;

        AB[0]=xyz_list[1][0] - xyz_list[0][0];
        AB[1]=xyz_list[1][1] - xyz_list[0][1];
        AB[2]=xyz_list[1][2] - xyz_list[0][2];
        AC[0]=xyz_list[2][0] - xyz_list[0][0];
        AC[1]=xyz_list[2][1] - xyz_list[0][1];
        AC[2]=xyz_list[2][2] - xyz_list[0][2];

        cross(normal,AB,AC);
        norm=sqrt(pow(normal[0],2.0)+pow(normal[1],2.0)+pow(normal[2],2.0));

        for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
}
/*}}}*/