Numericalflux.cpp

Go to the documentation of this file.

 
/*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
 
/*Numericalflux constructors and destructor*/
Numericalflux::Numericalflux(){/*{{{*/
   this->parameters = NULL;
   this->helement   = NULL;
   this->element    = NULL;
   this->hnodes     = NULL;
   this->hvertices  = NULL;
   this->nodes      = NULL;
}
/*}}}*/
Numericalflux::Numericalflux(int numericalflux_id,int i,int index,IoModel* iomodel){/*{{{*/
 
   /* Intermediary */
   int pos1,pos2,pos3,pos4;
   int numnodes;
 
   /*numericalflux constructor data: */
   int numericalflux_elem_ids[2];
   int numericalflux_vertex_ids[2];
   int numericalflux_node_ids[4];
   int numericalflux_type;
   int numericalflux_degree;
 
   /*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 */
      numericalflux_type=BoundaryEnum;
      numericalflux_elem_ids[0]=e1;
   }
   else{
      /* internal edge: connected to 2 elements */
      numericalflux_type=InternalEnum;
      numericalflux_elem_ids[0]=e1;
      numericalflux_elem_ids[1]=e2;
   }
 
   /*FIXME: hardcode element degree for now*/
   this->flux_degree= P1DGEnum;
   //this->flux_degree= P0DGEnum;
 
   /*1: Get vertices ids*/
   numericalflux_vertex_ids[0]=i1;
   numericalflux_vertex_ids[1]=i2;
 
   /*2: Get node ids*/
   if(numericalflux_type==InternalEnum){
      /*Get the column where these ids are located in the index*/
      pos1=pos2=pos3=pos4=UNDEF;
      for(int 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);
 
      /* 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]=3*(e1-1)+pos1;
      numericalflux_node_ids[1]=3*(e1-1)+pos2;
      numericalflux_node_ids[2]=3*(e2-1)+pos3;
      numericalflux_node_ids[3]=3*(e2-1)+pos4;
   }
   else{
      /*Get the column where these ids are located in the index*/
      pos1=pos2=UNDEF;
      for(int 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);
 
      /* 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]=3*(e1-1)+pos1;
      numericalflux_node_ids[1]=3*(e1-1)+pos2;
   }
 
   switch(this->flux_degree){
      case P0DGEnum:
         if(numericalflux_type==InternalEnum) numnodes = 2;
         else numnodes = 1;
         for(int i=0;i<numnodes;i++) numericalflux_node_ids[i] = numericalflux_elem_ids[i]; 
         numericalflux_node_ids[1] = numericalflux_elem_ids[1];
         break;
      case P1DGEnum:
         if(numericalflux_type==InternalEnum) numnodes = 4;
         else numnodes = 2;
         for(int i=0;i<numnodes;i++) numericalflux_node_ids[i] = numericalflux_node_ids[i]; //FIXME: to be improved...
         break;
      default:
         _error_("not supported yet");
 
   }
 
   /*Assign object fields: */
   this->id          = numericalflux_id;
   this->flux_type   = numericalflux_type;
 
   /*Hooks: */
   this->hnodes    = new Hook(numericalflux_node_ids,numnodes);
   this->hvertices = new Hook(&numericalflux_vertex_ids[0],2);
   this->helement  = new Hook(numericalflux_elem_ids,1); // take only the first element for now
 
   /*other fields*/
   this->parameters = NULL;
   this->element    = NULL;
   this->nodes      = NULL;
}
/*}}}*/
Numericalflux::~Numericalflux(){/*{{{*/
   this->parameters=NULL;
   delete helement;
   delete hnodes;
   delete hvertices;
}
/*}}}*/
 
/*Object virtual functions definitions:*/
Object* Numericalflux::copy() {/*{{{*/
 
   Numericalflux* numericalflux=NULL;
 
   numericalflux=new Numericalflux();
 
   /*copy fields: */
   numericalflux->id=this->id;
   numericalflux->flux_type=this->flux_type;
   numericalflux->flux_degree=this->flux_degree;
 
   /*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;
}
/*}}}*/
void    Numericalflux::DeepEcho(void){/*{{{*/
 
   _printf_("Numericalflux:\n");
   _printf_("   id: " << id << "\n");
   _printf_("   flux_type: " << this->flux_type<< "\n");
   _printf_("   flux_degree: " << this->flux_degree<< "\n");
   hnodes->DeepEcho();
   hvertices->DeepEcho();
   helement->DeepEcho();
   _printf_("   parameters\n");
   if(parameters)
    parameters->DeepEcho();
   else
    _printf_("      NULL\n");
}     
/*}}}*/
void    Numericalflux::Echo(void){/*{{{*/
   _printf_("Numericalflux:\n");
   _printf_("   id: " << id << "\n");
   _printf_("   flux_type: " << this->flux_type<< "\n");
   _printf_("   flux_degree: " << this->flux_degree<< "\n");
   hnodes->Echo();
   hvertices->Echo();
   helement->Echo();
   _printf_("   parameters: " << parameters << "\n");
}
/*}}}*/
int     Numericalflux::Id(void){/*{{{*/
   return id;
}
/*}}}*/
void    Numericalflux::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
 
   _assert_(this);
 
   /*ok, marshall operations: */
   MARSHALLING_ENUM(NumericalfluxEnum);
   MARSHALLING(id);
   MARSHALLING(flux_type);
   MARSHALLING(flux_degree);
 
   if(marshall_direction==MARSHALLING_BACKWARD){
      this->hnodes      = new Hook();
      this->hvertices   = new Hook();
      this->helement    = new Hook();
   }
 
   this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
   this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
 
   /*corresponding fields*/
   nodes    =(Node**)this->hnodes->deliverp();
   vertices =(Vertex**)this->hvertices->deliverp();
   element  =(Element*)this->helement->delivers();
 
}
/*}}}*/
int     Numericalflux::ObjectEnum(void){/*{{{*/
 
   return NumericalfluxEnum;
 
}
/*}}}*/
 
/*Load virtual functions definitions:*/
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;
}
/*}}}*/
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;
   }
 
}
/*}}}*/
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;
   }
 
}
/*}}}*/
void  Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/
 
   _assert_(lidlist);
   _assert_(nodes);
 
   int numnodes = this->GetNumberOfNodes();
   for(int i=0;i<numnodes;i++) lidlist[i]=nodes[i]->Lid();
}
/*}}}*/
void  Numericalflux::GetNodesSidList(int* sidlist){/*{{{*/
 
   _assert_(sidlist);
   _assert_(nodes);
 
   int numnodes = this->GetNumberOfNodes();
   for(int i=0;i<numnodes;i++) sidlist[i]=nodes[i]->Sid();
}
/*}}}*/
int   Numericalflux::GetNumberOfNodes(void){/*{{{*/
 
   if(this->flux_degree==P0DGEnum){
      switch(this->flux_type){
         case InternalEnum:
            return 2;
         case BoundaryEnum:
            return 1;
         default:
            _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
      }
   }
   else if(this->flux_degree==P1DGEnum){
      switch(this->flux_type){
         case InternalEnum:
            return 4;
         case BoundaryEnum:
            return 2;
         default:
            _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet");
      }
   }
   else{
      _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
   }
 
}
/*}}}*/
int   Numericalflux::GetNumberOfNodesOneSide(void){/*{{{*/
 
   if(this->flux_degree==P0DGEnum){
      return 1;
   }
   else if(this->flux_degree==P1DGEnum){
      return 2;
   }
   else{
      _error_("Numericalflux " << EnumToStringx(this->flux_degree) << " not supported yet");
   }
 
}
/*}}}*/
bool  Numericalflux::IsPenalty(void){/*{{{*/
   return false;
}
/*}}}*/
void  Numericalflux::PenaltyCreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,IssmDouble kmax){/*{{{*/
 
   /*No stiffness loads applied, do nothing: */
   return;
 
}
/*}}}*/
void  Numericalflux::PenaltyCreatePVector(Vector<IssmDouble>* pf,IssmDouble kmax){/*{{{*/
 
   /*No penalty loads applied, do nothing: */
   return;
 
}
/*}}}*/
void  Numericalflux::ResetHooks(){/*{{{*/
 
   this->nodes      = NULL;
   this->vertices   = NULL;
   this->element    = NULL;
   this->parameters = NULL;
 
   /*Get Element type*/
   this->hnodes->reset();
   this->hvertices->reset();
   this->helement->reset();
 
}
/*}}}*/
void  Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/
   /*Nothing to do :)*/
 
}
/*}}}*/
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]->fsize){
                  if(this->nodes[i]->IsClone())
                   o_nz += 1;
                  else
                   d_nz += 1;
               }
               break;
            case GsetEnum:
               if(nodes[i]->gsize){
                  if(this->nodes[i]->IsClone())
                   o_nz += 1;
                  else
                   d_nz += 1;
               }
               break;
            case SsetEnum:
               if(nodes[i]->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*/
ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/
 
   switch(this->flux_type){
      case InternalEnum:
         return CreateKMatrixAdjointBalancethicknessInternal();
      case BoundaryEnum:
         return CreateKMatrixAdjointBalancethicknessBoundary();
      default:
         _error_("type not supported yet");
   }
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/
 
   ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary();
   if(Ke) Ke->Transpose();
   return Ke;
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/
 
   ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal();
   if (Ke) Ke->Transpose();
   return Ke;
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/
 
   switch(this->flux_type){
      case InternalEnum:
         return CreateKMatrixBalancethicknessInternal();
      case BoundaryEnum:
         return CreateKMatrixBalancethicknessBoundary();
      default:
         _error_("type not supported yet");
   }
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/
 
   /*Initialize Element matrix and return if necessary*/
   Tria* tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
 
   /* Intermediaries*/
   IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
   IssmDouble xyz_list[NUMVERTICES][3];
   IssmDouble normal[2];
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input); 
   Input2* vyaverage_input=tria->GetInput2(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)*/
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
 
   GaussTria* gauss=new GaussTria();
   gauss->GaussEdgeCenter(index1,index2);
   vxaverage_input->GetInputValue(&mean_vx,gauss);
   vyaverage_input->GetInputValue(&mean_vy,gauss);
   delete gauss;
 
   IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
   if(UdotN<=0){
      return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
   }
 
   /*Initialize Element vector and other vectors*/
   int            numnodes = this->GetNumberOfNodes();
   ElementMatrix *Ke       = new ElementMatrix(nodes,numnodes,this->parameters);
   IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
 
   /* Start  looping on the number of gaussian points: */
   gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
 
      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;
 
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
         }
      }
   } 
 
   /*Clean up and return*/
   xDelete<IssmDouble>(basis);
   delete gauss;
   return Ke;
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/
 
   /*Initialize Element matrix and return if necessary*/
   Tria*  tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
 
   /* Intermediaries*/
   IssmDouble A1,A2,Jdet,vx,vy,UdotN;
   IssmDouble xyz_list[NUMVERTICES][3];
   IssmDouble normal[2];
 
   /*Fetch number of nodes for this flux*/
   int numnodes       = this->GetNumberOfNodes();
   int numnodes_plus  = this->GetNumberOfNodesOneSide();
   int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
   _assert_(numnodes==numnodes_plus+numnodes_minus);
 
   /*Initialize variables*/
   ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
   IssmDouble    *basis_plus  = xNew<IssmDouble>(numnodes_plus);
   IssmDouble    *basis_minus = xNew<IssmDouble>(numnodes_minus);
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input); 
   Input2* vyaverage_input=tria->GetInput2(VyEnum); _assert_(vyaverage_input); 
   GetNormal(&normal[0],xyz_list);
 
   /* Start  looping on the number of gaussian points: */
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
   GaussTria* gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
      tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
 
      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);
      A1=gauss->weight*Jdet*UdotN/2;
      A2=gauss->weight*Jdet*fabs(UdotN)/2;
 
      /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
       *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
       *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
       *
       *Term 2 (stabilization)  |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
       *                                      = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
       *     | A++ | A+- |
       * K = |-----------|
       *     | A-+ | A-- |
       *
       *These 4 terms for each expressions are added independently*/
 
      /*First term A++*/
      for(int i=0;i<numnodes_plus;i++){
         for(int j=0;j<numnodes_plus;j++){
            Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
            Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
         }
      }
      /*Second term A+-*/
      for(int i=0;i<numnodes_plus;i++){
         for(int j=0;j<numnodes_minus;j++){
            Ke->values[i*numnodes+numnodes_plus+j] +=  A1*(basis_minus[j]*basis_plus[i]);
            Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
         }
      }
      /*Third term A-+*/
      for(int i=0;i<numnodes_minus;i++){
         for(int j=0;j<numnodes_plus;j++){
            Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
            Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
         }
      }
      /*Fourth term A-+*/
      for(int i=0;i<numnodes_minus;i++){
         for(int j=0;j<numnodes_minus;j++){
            Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
            Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] +=  A2*(basis_minus[j]*basis_minus[i]);
         }
      }
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(basis_plus);
   xDelete<IssmDouble>(basis_minus);
   delete gauss;
   return Ke;
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/
 
   switch(this->flux_type){
      case InternalEnum:
         return CreateKMatrixMasstransportInternal();
      case BoundaryEnum:
         return CreateKMatrixMasstransportBoundary();
      default:
         _error_("type not supported yet");
   }
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/
 
   /*Initialize Element matrix and return if necessary*/
   Tria*  tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
 
   /* Intermediaries*/
   IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy;
   IssmDouble xyz_list[NUMVERTICES][3];
   IssmDouble normal[2];
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
   Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input);
   Input2* vyaverage_input=tria->GetInput2(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)*/
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
 
   GaussTria* gauss=new GaussTria();
   gauss->GaussEdgeCenter(index1,index2);
   vxaverage_input->GetInputValue(&mean_vx,gauss);
   vyaverage_input->GetInputValue(&mean_vy,gauss);
   delete gauss;
 
   IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
   if(UdotN<=0){
      return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
   }
 
   /*Initialize Element vector and other vectors*/
   int            numnodes = this->GetNumberOfNodes();
   ElementMatrix *Ke       = new ElementMatrix(nodes,numnodes,this->parameters);
   IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
 
   /* Start  looping on the number of gaussian points: */
   gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
 
      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;
 
      for(int i=0;i<numnodes;i++){
         for(int j=0;j<numnodes;j++){
            Ke->values[i*numnodes+j]+=DL*basis[i]*basis[j];
         }
      }
   } 
 
   /*Clean up and return*/
   xDelete<IssmDouble>(basis);
   delete gauss;
   return Ke;
}
/*}}}*/
ElementMatrix* Numericalflux::CreateKMatrixMasstransportInternal(void){/*{{{*/
 
   /*Initialize Element matrix and return if necessary*/
   Tria*  tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
 
   /* Intermediaries*/
   IssmDouble A1,A2,Jdet,vx,vy,UdotN;
   IssmDouble xyz_list[NUMVERTICES][3];
   IssmDouble normal[2];
 
   /*Fetch number of nodes for this flux*/
   int numnodes       = this->GetNumberOfNodes();
   int numnodes_plus  = this->GetNumberOfNodesOneSide();
   int numnodes_minus = numnodes_plus; /*For now we are not doing p-adaptive DG*/
   _assert_(numnodes==numnodes_plus+numnodes_minus);
 
   /*Initialize variables*/
   ElementMatrix *Ke = new ElementMatrix(nodes,numnodes,this->parameters);
   IssmDouble    *basis_plus  = xNew<IssmDouble>(numnodes_plus);
   IssmDouble    *basis_minus = xNew<IssmDouble>(numnodes_minus);
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
   Input2* vxaverage_input=tria->GetInput2(VxEnum); _assert_(vxaverage_input); 
   Input2* vyaverage_input=tria->GetInput2(VyEnum); _assert_(vyaverage_input); 
   GetNormal(&normal[0],xyz_list);
 
   /* Start  looping on the number of gaussian points: */
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
   GaussTria* gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      tria->GetSegmentNodalFunctions(&basis_plus[0] ,gauss,index1,index2,tria->FiniteElement());
      tria->GetSegmentNodalFunctions(&basis_minus[0],gauss,index1,index2,tria->FiniteElement());
 
      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);
      A1=gauss->weight*Jdet*dt*UdotN/2;
      A2=gauss->weight*Jdet*dt*fabs(UdotN)/2;
 
      /*Term 1 (numerical flux): {Hv}.[[phi]] = 0.5(H+v+ + H-v-)(phi+n+ + phi-n-)
       *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
       *                                      = v.n/2 (H+phi+ + H-phi+ -H+phi- -H-phi-)
       *
       *Term 2 (stabilization)  |v.n|/2 [[H]].[[phi]] = |v.n|/2 (H+n+ + H-n-)(phi+n+ + phi-n-)
       *                                      = |v.n|/2 (H+phi+ -H-phi+ -H+phi- +H-phi-)
       *     | A++ | A+- |
       * K = |-----------|
       *     | A-+ | A-- |
       *
       *These 4 terms for each expressions are added independently*/
 
      /*First term A++*/
      for(int i=0;i<numnodes_plus;i++){
         for(int j=0;j<numnodes_plus;j++){
            Ke->values[i*numnodes+j] += A1*(basis_plus[j]*basis_plus[i]);
            Ke->values[i*numnodes+j] += A2*(basis_plus[j]*basis_plus[i]);
         }
      }
      /*Second term A+-*/
      for(int i=0;i<numnodes_plus;i++){
         for(int j=0;j<numnodes_minus;j++){
            Ke->values[i*numnodes+numnodes_plus+j] +=  A1*(basis_minus[j]*basis_plus[i]);
            Ke->values[i*numnodes+numnodes_plus+j] += -A2*(basis_minus[j]*basis_plus[i]);
         }
      }
      /*Third term A-+*/
      for(int i=0;i<numnodes_minus;i++){
         for(int j=0;j<numnodes_plus;j++){
            Ke->values[(numnodes_plus+i)*numnodes+j] += -A1*(basis_plus[j]*basis_minus[i]);
            Ke->values[(numnodes_plus+i)*numnodes+j] += -A2*(basis_plus[j]*basis_minus[i]);
         }
      }
      /*Fourth term A-+*/
      for(int i=0;i<numnodes_minus;i++){
         for(int j=0;j<numnodes_minus;j++){
            Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] += -A1*(basis_minus[j]*basis_minus[i]);
            Ke->values[(numnodes_plus+i)*numnodes+numnodes_plus+j] +=  A2*(basis_minus[j]*basis_minus[i]);
         }
      }
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(basis_plus);
   xDelete<IssmDouble>(basis_minus);
   delete gauss;
   return Ke;
}
/*}}}*/
ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/
 
   /*No PVector for the Adjoint*/
   return NULL;
}
/*}}}*/
ElementVector* Numericalflux::CreatePVectorBalancethickness(void){/*{{{*/
 
   switch(this->flux_type){
      case InternalEnum:
         return CreatePVectorBalancethicknessInternal();
      case BoundaryEnum:
         return CreatePVectorBalancethicknessBoundary();
      default:
         _error_("type not supported yet");
   }
}
/*}}}*/
ElementVector* Numericalflux::CreatePVectorBalancethicknessBoundary(void){/*{{{*/
 
   /*Initialize Load Vector and return if necessary*/
   Tria*  tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
 
   /* Intermediaries*/
   IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
   IssmDouble xyz_list[NUMVERTICES][3];
   IssmDouble normal[2];
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   Input2* vxaverage_input = tria->GetInput2(VxEnum);        _assert_(vxaverage_input);
   Input2* vyaverage_input = tria->GetInput2(VyEnum);        _assert_(vyaverage_input);
   Input2* thickness_input = tria->GetInput2(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)*/
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
   GaussTria* gauss=new GaussTria();
   gauss->GaussEdgeCenter(index1,index2);
   vxaverage_input->GetInputValue(&mean_vx,gauss);
   vyaverage_input->GetInputValue(&mean_vy,gauss);
   delete gauss;
   IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
   if(UdotN>0){
      return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
   }
 
   /*Initialize Load Vector */
   int            numnodes = this->GetNumberOfNodes();
   ElementVector *pe       = new ElementVector(nodes,numnodes,this->parameters);
   IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
 
   /* Start  looping on the number of gaussian points: */
   gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
 
      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(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(basis);
   delete gauss;
   return pe;
}
/*}}}*/
ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/
 
   /*Nothing added to PVector*/
   return NULL;
 
}
/*}}}*/
ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/
 
   switch(this->flux_type){
      case InternalEnum:
         return CreatePVectorMasstransportInternal();
      case BoundaryEnum:
         return CreatePVectorMasstransportBoundary();
      default:
         _error_("type not supported yet");
   }
}
/*}}}*/
ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/
 
   /*Initialize Load Vector and return if necessary*/
   Tria* tria=(Tria*)element;
   if(!tria->IsIceInElement()) return NULL;
 
   /* Intermediaries*/
   IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,thickness;
   IssmDouble xyz_list[NUMVERTICES][3];
   IssmDouble normal[2];
 
   /*Retrieve all inputs and parameters*/
   GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
   IssmDouble dt = parameters->FindParam(TimesteppingTimeStepEnum);
   Input2* vxaverage_input    = tria->GetInput2(VxEnum);                        _assert_(vxaverage_input);
   Input2* vyaverage_input    = tria->GetInput2(VyEnum);                        _assert_(vyaverage_input);
   Input2* spcthickness_input = tria->GetInput2(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)*/
   int index1=tria->GetVertexIndex(vertices[0]);
   int index2=tria->GetVertexIndex(vertices[1]);
   GaussTria* gauss=new GaussTria();
   gauss->GaussEdgeCenter(index1,index2);
   vxaverage_input->GetInputValue(&mean_vx,gauss);
   vyaverage_input->GetInputValue(&mean_vy,gauss);
   delete gauss;
   IssmDouble UdotN=mean_vx*normal[0]+mean_vy*normal[1];
   if(UdotN>0){
      return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
   }
 
   /*Initialize Load Vector */
   int            numnodes = this->GetNumberOfNodes();
   ElementVector *pe       = new ElementVector(nodes,numnodes,this->parameters);
   IssmDouble    *basis    = xNew<IssmDouble>(numnodes);
 
   /* Start  looping on the number of gaussian points: */
   gauss=new GaussTria(index1,index2,2);
   for(int ig=gauss->begin();ig<gauss->end();ig++){
 
      gauss->GaussPoint(ig);
 
      tria->GetSegmentNodalFunctions(&basis[0],gauss,index1,index2,tria->FiniteElement());
 
      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(int i=0;i<numnodes;i++) pe->values[i] += DL*basis[i];
   }
 
   /*Clean up and return*/
   xDelete<IssmDouble>(basis);
   delete gauss;
   return pe;
}
/*}}}*/
ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/
 
   /*Nothing added to PVector*/
   return NULL;
 
}
/*}}}*/
void           Numericalflux::GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/
 
   /*Build unit outward pointing vector*/
   IssmDouble vector[2];
 
   vector[0]=xyz_list[1][0] - xyz_list[0][0];
   vector[1]=xyz_list[1][1] - xyz_list[0][1];
 
   IssmDouble norm=sqrt(pow(vector[0],2.0)+pow(vector[1],2.0));
 
   normal[0]= + vector[1]/norm;
   normal[1]= - vector[0]/norm;
}
/*}}}*/