source: issm/trunk/src/c/classes/Loads/Numericalflux.cpp@ 16560

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

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

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

/*Load macros*/
#define NUMVERTICES 2
#define NUMNODES_INTERNAL 4
#define NUMNODES_BOUNDARY 2

/*Numericalflux constructors and destructor*/
/*FUNCTION Numericalflux::Numericalflux(){{{*/
Numericalflux::Numericalflux(){
        this->inputs     = NULL;
        this->parameters = NULL;
        this->helement   = NULL;
        this->element    = NULL;
        this->hnodes     = NULL;
        this->hvertices  = NULL;
        this->nodes      = NULL;
}
/*}}}*/
/*FUNCTION Numericalflux::Numericalflux(int id, int i, IoModel* iomodel, int analysis_type) {{{*/
Numericalflux::Numericalflux(int numericalflux_id,int i,int index,IoModel* iomodel, int in_analysis_type){

        /* Intermediary */
        int  j;
        int  pos1,pos2,pos3,pos4;
        int  num_nodes;

        /*numericalflux constructor data: */
        int   numericalflux_elem_ids[2];
        int   numericalflux_vertex_ids[2];
        int   numericalflux_node_ids[4];
        int   numericalflux_type;

        /*Get edge*/
        int i1 = iomodel->faces[4*index+0];
        int i2 = iomodel->faces[4*index+1];
        int e1 = iomodel->faces[4*index+2];
        int e2 = iomodel->faces[4*index+3];

        /*First, see wether this is an internal or boundary edge (if e2=-1)*/
        if(e2==-1){
                /* Boundary edge, only one element */
                num_nodes=2;
                numericalflux_type=BoundaryEnum;
                numericalflux_elem_ids[0]=e1;
        }
        else{
                /* internal edge: connected to 2 elements */
                 num_nodes=4;
                numericalflux_type=InternalEnum;
                numericalflux_elem_ids[0]=e1;
                numericalflux_elem_ids[1]=e2;
        }

        /*1: Get vertices ids*/
        numericalflux_vertex_ids[0]=i1;
        numericalflux_vertex_ids[1]=i2;

        /*2: Get node ids*/
        if (numericalflux_type==InternalEnum){

                /*Now, we must get the nodes of the 4 nodes located on the edge*/

                /*2: Get the column where these ids are located in the index*/
                pos1=pos2=pos3=pos4=UNDEF;
                for(j=0;j<3;j++){
                        if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
                        if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
                        if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1;
                        if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1;
                }
                _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF);

                /*3: We have the id of the elements and the position of the vertices in the index
                 * we can compute their dofs!*/
                numericalflux_node_ids[0]=iomodel->nodecounter+3*(e1-1)+pos1;
                numericalflux_node_ids[1]=iomodel->nodecounter+3*(e1-1)+pos2;
                numericalflux_node_ids[2]=iomodel->nodecounter+3*(e2-1)+pos3;
                numericalflux_node_ids[3]=iomodel->nodecounter+3*(e2-1)+pos4;
        }
        else{

                /*2: Get the column where these ids are located in the index*/
                pos1=pos2=UNDEF;
                for(j=0;j<3;j++){
                        if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1;
                        if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1;
                }
                _assert_(pos1!=UNDEF && pos2!=UNDEF);

                /*3: We have the id of the elements and the position of the vertices in the index
                 * we can compute their dofs!*/
                numericalflux_node_ids[0]=iomodel->nodecounter+3*(e1-1)+pos1;
                numericalflux_node_ids[1]=iomodel->nodecounter+3*(e1-1)+pos2;
        }

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

        /*Hooks: */
        this->hnodes    =new Hook(numericalflux_node_ids,num_nodes);
        this->hvertices =new Hook(&numericalflux_vertex_ids[0],2);
        this->helement  =new Hook(numericalflux_elem_ids,1); // take only the first element for now

        //intialize  and add as many inputs per element as requested: 
        this->inputs=new Inputs();
        this->inputs->AddInput(new IntInput(NumericalfluxTypeEnum,numericalflux_type));

        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
        this->parameters=NULL;
        this->element=NULL;
        this->nodes=NULL;
}
/*}}}*/
/*FUNCTION Numericalflux::~Numericalflux(){{{*/
Numericalflux::~Numericalflux(){
        delete inputs;
        this->parameters=NULL;
        delete helement;
        delete hnodes;
        delete hvertices;
}
/*}}}*/

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

        _printf_("Numericalflux:\n");
        _printf_("   id: " << id << "\n");
        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
        hnodes->DeepEcho();
        hvertices->DeepEcho();
        helement->DeepEcho();
        _printf_("   parameters\n");
        if(parameters)
         parameters->DeepEcho();
        else
         _printf_("      NULL\n");
        _printf_("   inputs\n");
        inputs->DeepEcho();

}               
/*}}}*/
/*FUNCTION Numericalflux::Id {{{*/
int    Numericalflux::Id(void){
        return id;
}
/*}}}*/
/*FUNCTION Numericalflux::ObjectEnum{{{*/
int Numericalflux::ObjectEnum(void){

        return NumericalfluxEnum;

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

        Numericalflux* numericalflux=NULL;

        numericalflux=new Numericalflux();

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

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

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

        return numericalflux;
}
/*}}}*/

/*Load virtual functions definitions:*/
/*FUNCTION Numericalflux::Configure {{{*/
void  Numericalflux::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);

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

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

}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrix {{{*/
void  Numericalflux::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){

        /*recover some parameters*/
        ElementMatrix* Ke=NULL;
        int analysis_type;
        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);

        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
        switch(analysis_type){
                case MasstransportAnalysisEnum:
                        Ke=CreateKMatrixMasstransport();
                        break;
                case BalancethicknessAnalysisEnum:
                        Ke=CreateKMatrixBalancethickness();
                        break;
                case AdjointBalancethicknessAnalysisEnum:
                        Ke=CreateKMatrixAdjointBalancethickness();
                        break;
                default:
                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
        }

        /*Add to global matrix*/
        if(Ke){
                Ke->AddToGlobal(Kff,Kfs);
                delete Ke;
        }

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

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

        switch(analysis_type){
                case MasstransportAnalysisEnum:
                        pe=CreatePVectorMasstransport();
                        break;
                case BalancethicknessAnalysisEnum:
                        pe=CreatePVectorBalancethickness();
                        break;
                case AdjointBalancethicknessAnalysisEnum:
                        pe=CreatePVectorAdjointBalancethickness();
                        break;
                default:
                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
        }

        /*Add to global matrix*/
        if(pe){
                pe->AddToGlobal(pf);
                delete pe;
        }

}
/*}}}*/
/*FUNCTION Numericalflux::GetNodesSidList{{{*/
void Numericalflux::GetNodesSidList(int* sidlist){

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

        switch(type){
                case InternalEnum:
                        for(int i=0;i<NUMNODES_INTERNAL;i++) sidlist[i]=nodes[i]->Sid();
                        return;
                case BoundaryEnum:
                        for(int i=0;i<NUMNODES_BOUNDARY;i++) sidlist[i]=nodes[i]->Sid();
                        return;
                default:
                        _error_("Numericalflux type " << EnumToStringx(type) << " not supported yet");
        }
}
/*}}}*/
/*FUNCTION Numericalflux::GetNodesLidList{{{*/
void Numericalflux::GetNodesLidList(int* lidlist){

        int type;
        inputs->GetInputValue(&type,NumericalfluxTypeEnum);
        _assert_(lidlist);
        _assert_(nodes);

        switch(type){
                case InternalEnum:
                        for(int i=0;i<NUMNODES_INTERNAL;i++) lidlist[i]=nodes[i]->Lid();
                        return;
                case BoundaryEnum:
                        for(int i=0;i<NUMNODES_BOUNDARY;i++) lidlist[i]=nodes[i]->Lid();
                        return;
                default:
                        _error_("Numericalflux type " << EnumToStringx(type) << " not supported yet");
        }
}
/*}}}*/
/*FUNCTION Numericalflux::GetNumberOfNodes{{{*/
int Numericalflux::GetNumberOfNodes(void){

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

        switch(type){
                case InternalEnum:
                        return NUMNODES_INTERNAL;
                case BoundaryEnum:
                        return NUMNODES_BOUNDARY;
                default:
                        _error_("Numericalflux type " << EnumToStringx(type) << " not supported yet");
        }

}
/*}}}*/
/*FUNCTION Numericalflux::IsPenalty{{{*/
bool Numericalflux::IsPenalty(void){
        return false;
}
/*}}}*/
/*FUNCTION Numericalflux::PenaltyCreateKMatrix {{{*/
void  Numericalflux::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){

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

}
/*}}}*/
/*FUNCTION Numericalflux::PenaltyCreatePVector{{{*/
void  Numericalflux::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){

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

}
/*}}}*/
/*FUNCTION Numericalflux::InAnalysis{{{*/
bool Numericalflux::InAnalysis(int in_analysis_type){
        if (in_analysis_type==this->analysis_type) return true;
        else return false;
}
/*}}}*/
/*FUNCTION Numericalflux::SetwiseNodeConnectivity{{{*/
void Numericalflux::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,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]->Lid()]){

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

                        int counter=0;
                        while(flagsindices[counter]>=0) counter++;
                        flagsindices[counter]=this->nodes[i]->Lid();

                        /*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;
}
/*}}}*/

/*Numericalflux management*/
/*FUNCTION Numericalflux::CreateKMatrixMasstransport{{{*/
ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){

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

        switch(type){
                case InternalEnum:
                        return CreateKMatrixMasstransportInternal();
                case BoundaryEnum:
                        return CreateKMatrixMasstransportBoundary();
                default:
                        _error_("type not supported yet");
        }
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixMasstransportInternal {{{*/
ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){

        /* constants*/
        const int numdof=NDOF1*NUMNODES_INTERNAL;

        /* Intermediaries*/
        int        i,j,ig,index1,index2;
        IssmDouble     DL1,DL2,Jdet,dt,vx,vy,UdotN;
        IssmDouble     xyz_list[NUMVERTICES][3];
        IssmDouble     normal[2];
        IssmDouble     B[numdof];
        IssmDouble     Bprime[numdof];
        IssmDouble     Ke_g1[numdof][numdof];
        IssmDouble     Ke_g2[numdof][numdof];
        GaussTria *gauss;

        /*Initialize Element matrix and return if necessary*/
        Tria*  tria=(Tria*)element;
        if(tria->NoIceInElement()) return NULL;
        ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
        GetNormal(&normal[0],xyz_list);

        /* Start  looping on the number of gaussian points: */
        index1=tria->GetNodeIndex(nodes[0]);
        index2=tria->GetNodeIndex(nodes[1]);
        gauss=new GaussTria(index1,index2,2);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
                tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);

                vxaverage_input->GetInputValue(&vx,gauss);
                vyaverage_input->GetInputValue(&vy,gauss);
                UdotN=vx*normal[0]+vy*normal[1];
                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
                DL1=gauss->weight*Jdet*dt*UdotN/2;
                DL2=gauss->weight*Jdet*dt*fabs(UdotN)/2;

                TripleMultiply(&B[0],1,numdof,1,
                                        &DL1,1,1,0,
                                        &Bprime[0],1,numdof,0,
                                        &Ke_g1[0][0],0);
                TripleMultiply(&B[0],1,numdof,1,
                                        &DL2,1,1,0,
                                        &B[0],1,numdof,0,
                                        &Ke_g2[0][0],0);

                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
        }

        /*Clean up and return*/
        delete gauss;
        return Ke;
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixMasstransportBoundary {{{*/
ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){

        /* constants*/
        const int numdof=NDOF1*NUMNODES_BOUNDARY;

        /* Intermediaries*/
        int        i,j,ig,index1,index2;
        IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN;
        IssmDouble     xyz_list[NUMVERTICES][3];
        IssmDouble     normal[2];
        IssmDouble     L[numdof];
        IssmDouble     Ke_g[numdof][numdof];
        GaussTria *gauss;

        /*Initialize Element matrix and return if necessary*/
        ElementMatrix* Ke = NULL;
        Tria*  tria=(Tria*)element;
        if(tria->NoIceInElement()) return NULL;

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input);
        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
        GetNormal(&normal[0],xyz_list);

        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
        index1=tria->GetNodeIndex(nodes[0]);
        index2=tria->GetNodeIndex(nodes[1]);

        gauss=new GaussTria();
        gauss->GaussEdgeCenter(index1,index2);
        vxaverage_input->GetInputValue(&mean_vx,gauss);
        vyaverage_input->GetInputValue(&mean_vy,gauss);
        delete gauss;

        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
        if (UdotN<=0){
                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
        }
        else{
                Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
        }

        /* Start  looping on the number of gaussian points: */
        gauss=new GaussTria(index1,index2,2);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);

                vxaverage_input->GetInputValue(&vx,gauss);
                vyaverage_input->GetInputValue(&vy,gauss);
                UdotN=vx*normal[0]+vy*normal[1];
                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
                DL=gauss->weight*Jdet*dt*UdotN;

                TripleMultiply(&L[0],1,numdof,1,
                                        &DL,1,1,0,
                                        &L[0],1,numdof,0,
                                        &Ke_g[0][0],0);

                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
        } 

        /*Clean up and return*/
        delete gauss;
        return Ke;
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixBalancethickness{{{*/
ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){

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

        switch(type){
                case InternalEnum:
                        return CreateKMatrixBalancethicknessInternal();
                case BoundaryEnum:
                        return CreateKMatrixBalancethicknessBoundary();
                default:
                        _error_("type not supported yet");
        }
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixBalancethicknessInternal {{{*/
ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){

        /* constants*/
        const int numdof=NDOF1*NUMNODES_INTERNAL;

        /* Intermediaries*/
        int        i,j,ig,index1,index2;
        IssmDouble     DL1,DL2,Jdet,vx,vy,UdotN;
        IssmDouble     xyz_list[NUMVERTICES][3];
        IssmDouble     normal[2];
        IssmDouble     B[numdof];
        IssmDouble     Bprime[numdof];
        IssmDouble     Ke_g1[numdof][numdof];
        IssmDouble     Ke_g2[numdof][numdof];
        GaussTria *gauss;

        /*Initialize Element matrix and return if necessary*/
        Tria*  tria=(Tria*)element;
        if(tria->NoIceInElement()) return NULL;
        ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters);

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
        GetNormal(&normal[0],xyz_list);

        /* Start  looping on the number of gaussian points: */
        index1=tria->GetNodeIndex(nodes[0]);
        index2=tria->GetNodeIndex(nodes[1]);
        gauss=new GaussTria(index1,index2,2);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                tria->GetSegmentBFlux(&B[0],gauss,index1,index2);
                tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2);

                vxaverage_input->GetInputValue(&vx,gauss);
                vyaverage_input->GetInputValue(&vy,gauss);
                UdotN=vx*normal[0]+vy*normal[1];
                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
                DL1=gauss->weight*Jdet*UdotN/2;
                DL2=gauss->weight*Jdet*fabs(UdotN)/2;

                TripleMultiply(&B[0],1,numdof,1,
                                        &DL1,1,1,0,
                                        &Bprime[0],1,numdof,0,
                                        &Ke_g1[0][0],0);
                TripleMultiply(&B[0],1,numdof,1,
                                        &DL2,1,1,0,
                                        &B[0],1,numdof,0,
                                        &Ke_g2[0][0],0);

                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g1[i][j];
                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g2[i][j];
        }

        /*Clean up and return*/
        delete gauss;
        return Ke;
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixBalancethicknessBoundary {{{*/
ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){

        /* constants*/
        const int numdof=NDOF1*NUMNODES_BOUNDARY;

        /* Intermediaries*/
        int        i,j,ig,index1,index2;
        IssmDouble     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN;
        IssmDouble     xyz_list[NUMVERTICES][3];
        IssmDouble     normal[2];
        IssmDouble     L[numdof];
        IssmDouble     Ke_g[numdof][numdof];
        GaussTria *gauss;

        /*Initialize Element matrix and return if necessary*/
        ElementMatrix* Ke = NULL;
        Tria*  tria=(Tria*)element;
        if(tria->NoIceInElement()) return NULL;

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);
        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);
        GetNormal(&normal[0],xyz_list);

        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
        index1=tria->GetNodeIndex(nodes[0]);
        index2=tria->GetNodeIndex(nodes[1]);

        gauss=new GaussTria();
        gauss->GaussEdgeCenter(index1,index2);
        vxaverage_input->GetInputValue(&mean_vx,gauss);
        vyaverage_input->GetInputValue(&mean_vy,gauss);
        delete gauss;

        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
        if (UdotN<=0){
                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
        }
        else{
                Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters);
        }

        /* Start  looping on the number of gaussian points: */
        gauss=new GaussTria(index1,index2,2);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);

                vxaverage_input->GetInputValue(&vx,gauss);
                vyaverage_input->GetInputValue(&vy,gauss);
                UdotN=vx*normal[0]+vy*normal[1];
                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
                DL=gauss->weight*Jdet*UdotN;

                TripleMultiply(&L[0],1,numdof,1,
                                        &DL,1,1,0,
                                        &L[0],1,numdof,0,
                                        &Ke_g[0][0],0);

                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j];
        } 

        /*Clean up and return*/
        delete gauss;
        return Ke;
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancethickness{{{*/
ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){

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

        switch(type){
                case InternalEnum:
                        return CreateKMatrixAdjointBalancethicknessInternal();
                case BoundaryEnum:
                        return CreateKMatrixAdjointBalancethicknessBoundary();
                default:
                        _error_("type not supported yet");
        }
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancethicknessInternal {{{*/
ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){

        ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
        if (Ke) Ke->Transpose();
        return Ke;
}
/*}}}*/
/*FUNCTION Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary {{{*/
ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){

        ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
        if(Ke) Ke->Transpose();
        return Ke;
}
/*}}}*/
/*FUNCTION Numericalflux::CreatePVectorMasstransport{{{*/
ElementVector* Numericalflux::CreatePVectorMasstransport(void){

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

        switch(type){
                case InternalEnum:
                        return CreatePVectorMasstransportInternal();
                case BoundaryEnum:
                        return CreatePVectorMasstransportBoundary();
                default:
                        _error_("type not supported yet");
        }
}
/*}}}*/
/*FUNCTION Numericalflux::CreatePVectorMasstransportInternal{{{*/
ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){

        /*Nothing added to PVector*/
        return NULL;

}
/*}}}*/
/*FUNCTION Numericalflux::CreatePVectorMasstransportBoundary{{{*/
ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){

        /* constants*/
        const int numdof=NDOF1*NUMNODES_BOUNDARY;

        /* Intermediaries*/
        int        i,ig,index1,index2;
        IssmDouble     DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness;
        IssmDouble     xyz_list[NUMVERTICES][3];
        IssmDouble     normal[2];
        IssmDouble     L[numdof];
        GaussTria *gauss;

        /*Initialize Load Vector and return if necessary*/
        ElementVector* pe = NULL;
        Tria*  tria=(Tria*)element;
        if(tria->NoIceInElement()) return NULL;

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
        Input* vxaverage_input   =tria->inputs->GetInput(VxEnum);                     _assert_(vxaverage_input); 
        Input* vyaverage_input   =tria->inputs->GetInput(VyEnum);                     _assert_(vyaverage_input);
        Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input);
        GetNormal(&normal[0],xyz_list);

        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
        index1=tria->GetNodeIndex(nodes[0]);
        index2=tria->GetNodeIndex(nodes[1]);

        gauss=new GaussTria();
        gauss->GaussEdgeCenter(index1,index2);
        vxaverage_input->GetInputValue(&mean_vx,gauss);
        vyaverage_input->GetInputValue(&mean_vy,gauss);
        delete gauss;

        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
        if (UdotN>0){
                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
        }
        else{
                pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
        }

        /* Start  looping on the number of gaussian points: */
        gauss=new GaussTria(index1,index2,2);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);

                vxaverage_input->GetInputValue(&vx,gauss);
                vyaverage_input->GetInputValue(&vy,gauss);
                spcthickness_input->GetInputValue(&thickness,gauss);
                if(xIsNan<IssmDouble>(thickness)) _error_("Cannot weakly apply constraint because NaN was provided");

                UdotN=vx*normal[0]+vy*normal[1];
                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
                DL= - gauss->weight*Jdet*dt*UdotN*thickness;

                for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
        }

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

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

        switch(type){
                case InternalEnum:
                        return CreatePVectorBalancethicknessInternal();
                case BoundaryEnum:
                        return CreatePVectorBalancethicknessBoundary();
                default:
                        _error_("type not supported yet");
        }
}
/*}}}*/
/*FUNCTION Numericalflux::CreatePVectorBalancethicknessInternal{{{*/
ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){

        /*Nothing added to PVector*/
        return NULL;

}
/*}}}*/
/*FUNCTION Numericalflux::CreatePVectorBalancethicknessBoundary{{{*/
ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){

        /* constants*/
        const int numdof=NDOF1*NUMNODES_BOUNDARY;

        /* Intermediaries*/
        int        i,ig,index1,index2;
        IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
        IssmDouble xyz_list[NUMVERTICES][3];
        IssmDouble normal[2];
        IssmDouble L[numdof];
        GaussTria *gauss;

        /*Initialize Load Vector and return if necessary*/
        ElementVector* pe = NULL;
        Tria*  tria=(Tria*)element;
        if(tria->NoIceInElement()) return NULL;

        /*Retrieve all inputs and parameters*/
        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); 
        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input);
        Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
        GetNormal(&normal[0],xyz_list);

        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
        index1=tria->GetNodeIndex(nodes[0]);
        index2=tria->GetNodeIndex(nodes[1]);

        gauss=new GaussTria();
        gauss->GaussEdgeCenter(index1,index2);
        vxaverage_input->GetInputValue(&mean_vx,gauss);
        vyaverage_input->GetInputValue(&mean_vy,gauss);
        delete gauss;
        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
        if (UdotN>0){
                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
        }
        else{
                pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters);
        }

        /* Start  looping on the number of gaussian points: */
        gauss=new GaussTria(index1,index2,2);
        for(ig=gauss->begin();ig<gauss->end();ig++){

                gauss->GaussPoint(ig);

                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);

                vxaverage_input->GetInputValue(&vx,gauss);
                vyaverage_input->GetInputValue(&vy,gauss);
                thickness_input->GetInputValue(&thickness,gauss);

                UdotN=vx*normal[0]+vy*normal[1];
                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
                DL= - gauss->weight*Jdet*UdotN*thickness;

                for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
        }

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

        /*No PVector for the Adjoint*/
        return NULL;
}
/*}}}*/
/*FUNCTION Numericalflux::GetNormal {{{*/
void Numericalflux:: GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){

        /*Build unit outward pointing vector*/
        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;
}
/*}}}*/