Changeset 24091


Ignore:
Timestamp:
07/15/19 12:50:06 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: more flexible masstransport

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

    r23629 r24091  
    624624void           DamageEvolutionAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    625625
    626         int domaintype;
    627         IssmDouble  max_damage;
    628         int                     *doflist = NULL;
    629 
    630         element->FindParam(&domaintype,DomainTypeEnum);
    631 
    632626        /*Fetch number of nodes and dof for this finite element*/
    633627        int numnodes = element->GetNumberOfNodes();
    634628
    635629        /*Fetch dof list and allocate solution vector*/
     630        int *doflist = NULL;
    636631        element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
    637632        IssmDouble* newdamage = xNew<IssmDouble>(numnodes);
    638633
    639634        /*Get user-supplied max_damage: */
    640         element->FindParam(&max_damage,DamageMaxDamageEnum);
     635        IssmDouble max_damage = element->FindParam(DamageMaxDamageEnum);
    641636
    642637        /*Use the dof list to index into the solution vector: */
    … …  
    652647
    653648        /*Get all inputs and parameters*/
     649        int domaintype;
     650        element->FindParam(&domaintype,DomainTypeEnum);
    654651        if(domaintype==Domain2DhorizontalEnum){
    655652                element->AddInput(DamageDbarEnum,newdamage,element->GetElementType());

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r24045 r24091  
    55#include "../modules/modules.h"
    66
     7#define FINITEELEMENT P1Enum
     8
    79/*Model processing*/
    810void MasstransportAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    … …  
    1416        /*Do not add constraints in DG,  they are weakly imposed*/
    1517        if(stabilization!=3){
    16                 IoModelToConstraintsx(constraints,iomodel,"md.masstransport.spcthickness",MasstransportAnalysisEnum,P1Enum);
     18                IoModelToConstraintsx(constraints,iomodel,"md.masstransport.spcthickness",MasstransportAnalysisEnum,FINITEELEMENT);
    1719        }
    1820
    … …  
    2527
    2628        /*Intermediaries*/
    27         int element;
    2829        int penpair_ids[2];
    2930        int count=0;
    … …  
    3536
    3637        /*Loads only in DG*/
    37         if (stabilization==3){
     38        if(stabilization==3){
    3839
    3940                /*Get faces and elements*/
    … …  
    4546
    4647                        /*Get left and right elements*/
    47                         element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     48                        int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
    4849
    4950                        /*Now, if this element is not in the partition, pass: */
    5051                        if(!iomodel->my_elements[element]) continue;
    51 
    52                         /* Add load */
    5352                        loads->AddObject(new Numericalflux(i+1,i,i,iomodel));
    5453                }
    … …  
    102101        if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    103102        if(stabilization!=3){
    104                 ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,P1Enum,isamr);
     103                ::CreateNodes(nodes,iomodel,MasstransportAnalysisEnum,FINITEELEMENT,isamr);
    105104        }
    106105        else{
    … …  
    130129
    131130        /*Finite element type*/
    132         finiteelement = P1Enum;
     131        finiteelement = FINITEELEMENT;
    133132        if(stabilization==3){
    134133                finiteelement = P1DGEnum;
     134                //finiteelement = P0DGEnum;
    135135        }
    136136
    … …  
    260260                        Ke = CreateKMatrixCG(basalelement);
    261261                        break;
     262                case P0DGEnum:
    262263                case P1DGEnum:
    263264                        Ke = CreateKMatrixDG(basalelement);
    … …  
    549550                        pe = CreatePVectorCG(basalelement);
    550551                        break;
     552                case P0DGEnum:
    551553                case P1DGEnum:
    552554                        pe = CreatePVectorDG(basalelement);
    … …  
    706708        element->GetVerticesCoordinates(&xyz_list);
    707709        element->FindParam(&dt,TimesteppingTimeStepEnum);
    708         Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
    709         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
    710         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);          _assert_(ms_input);
    711         Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
    712         Input* thickness_input     = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
     710        Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
     711        Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
     712        Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                      _assert_(ms_input);
     713        Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
     714        Input* thickness_input  = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
    713715
    714716        /* Start  looping on the number of gaussian points: */
    … …  
    803805void           MasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    804806
    805         int        i,hydroadjustment,domaintype;
    806         int*       doflist=NULL;
    807         IssmDouble rho_ice,rho_water,minthickness;
    808         Element*   basalelement=NULL;
    809         bool        isgroundingline=0;
    810 
    811         element->FindParam(&domaintype,DomainTypeEnum);
    812         if(domaintype!=Domain2DhorizontalEnum){
    813                 if(!element->IsOnBase()) return;
    814                 basalelement=element->SpawnBasalElement();
    815         }
    816         else{
    817                 basalelement = element;
    818         }
    819 
    820         /*Fetch number of nodes for this finite element*/
    821         int numnodes = basalelement->GetNumberOfNodes();
     807        /*Only update if on base*/
     808        if(!element->IsOnBase()) return;
     809
     810        /*Get basal element*/
     811        int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
     812   Element* basalelement=element;
     813        if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
     814
     815        /*Fetch number of nodes and dof for this finite element*/
     816        int numnodes    = basalelement->GetNumberOfNodes();
     817   int numvertices = basalelement->GetNumberOfVertices();
     818
     819   /*Keep old thickness for later*/
     820        IssmDouble* oldthickness = xNew<IssmDouble>(numvertices);
     821   basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
    822822
    823823        /*Fetch dof list and allocate solution vector*/
    824         basalelement->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
    825         IssmDouble* newthickness   = xNew<IssmDouble>(numnodes);
    826         IssmDouble* cumdeltathickness = xNew<IssmDouble>(numnodes);
    827         IssmDouble* deltathickness    = xNew<IssmDouble>(numnodes);
    828         IssmDouble* newbase        = xNew<IssmDouble>(numnodes);
    829         IssmDouble* bed            = xNew<IssmDouble>(numnodes);
    830         IssmDouble* newsurface     = xNew<IssmDouble>(numnodes);
    831         IssmDouble* oldthickness   = xNew<IssmDouble>(numnodes);
    832         IssmDouble* oldbase        = xNew<IssmDouble>(numnodes);
    833         IssmDouble* oldsurface     = xNew<IssmDouble>(numnodes);
    834         IssmDouble* phi            = xNew<IssmDouble>(numnodes);
    835         IssmDouble* sealevel       = xNew<IssmDouble>(numnodes);
     824        int *doflist = NULL;
     825        element->GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
     826        IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
    836827
    837828        /*Use the dof list to index into the solution vector: */
    838         basalelement->FindParam(&minthickness,MasstransportMinThicknessEnum);
    839         for(i=0;i<numnodes;i++){
     829   IssmDouble minthickness = basalelement->FindParam(MasstransportMinThicknessEnum);
     830        for(int i=0;i<numnodes;i++){
    840831                newthickness[i]=solution[doflist[i]];
     832                /*Check solution*/
    841833                if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
    842834                if(xIsInf<IssmDouble>(newthickness[i])) _error_("Inf found in solution vector");
    843                 /*Constrain thickness to be at least 1m*/
    844                 if(newthickness[i]<minthickness) newthickness[i]=minthickness;
    845         }
     835      if(newthickness[i]<minthickness) newthickness[i]=minthickness;
     836        }
     837
     838        element->AddBasalInput(ThicknessEnum,newthickness,element->GetElementType());
     839
     840        /*Free ressources:*/
     841        xDelete<IssmDouble>(newthickness);
     842        xDelete<int>(doflist);
     843
     844        /*Now, we need to do some "processing"*/
     845   newthickness  = xNew<IssmDouble>(numvertices);
     846        IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices);
     847        IssmDouble* deltathickness    = xNew<IssmDouble>(numvertices);
     848        IssmDouble* newbase           = xNew<IssmDouble>(numvertices);
     849        IssmDouble* bed               = xNew<IssmDouble>(numvertices);
     850        IssmDouble* newsurface        = xNew<IssmDouble>(numvertices);
     851        IssmDouble* oldbase           = xNew<IssmDouble>(numvertices);
     852        IssmDouble* oldsurface        = xNew<IssmDouble>(numvertices);
     853        IssmDouble* phi               = xNew<IssmDouble>(numvertices);
     854        IssmDouble* sealevel          = xNew<IssmDouble>(numvertices);
    846855
    847856        /*Get previous base, thickness, surfac and current sealevel and bed:*/
    848         basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum);
    849         basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum);
    850         basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum);
    851         basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum);
    852         basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum);
    853         basalelement->GetInputListOnNodes(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
    854 
     857   basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum);
     858        basalelement->GetInputListOnVertices(&oldbase[0],BaseEnum);
     859        basalelement->GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
     860        basalelement->GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
     861        basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum);
     862        basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
     863
     864   /*Do we do grounding line migration?*/
     865        bool isgroundingline;
    855866        element->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
    856         if(isgroundingline) basalelement->GetInputListOnNodes(&bed[0],BedEnum);
     867        if(isgroundingline) basalelement->GetInputListOnVertices(&bed[0],BedEnum);
    857868
    858869        /*What is the delta thickness forcing the sea-level rise core: cumulated over time, hence the +=:*/
    859         for(i=0;i<numnodes;i++){
     870        for(int i=0;i<numvertices;i++){
    860871                cumdeltathickness[i]+=newthickness[i]-oldthickness[i];
    861872                deltathickness[i]=newthickness[i]-oldthickness[i];
    … …  
    863874
    864875        /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
     876   int hydroadjustment;
    865877        basalelement->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
    866         rho_ice   = basalelement->FindParam(MaterialsRhoIceEnum);
    867         rho_water = basalelement->FindParam(MaterialsRhoSeawaterEnum);
    868 
    869         for(i=0;i<numnodes;i++) {
     878        IssmDouble rho_ice   = basalelement->FindParam(MaterialsRhoIceEnum);
     879        IssmDouble rho_water = basalelement->FindParam(MaterialsRhoSeawaterEnum);
     880
     881        for(int i=0;i<numvertices;i++) {
    870882                if (phi[i]>0.){ //this is an ice sheet: just add thickness to base.
    871883                        /*Update! actually, the bed has moved too!:*/
    872884                        if(isgroundingline){
    873885                                newsurface[i] = bed[i]+newthickness[i]; //surface = bed + newthickness
    874                                 newbase[i]     = bed[i];                 //new base at new bed
     886                                newbase[i]    = bed[i];                 //new base at new bed
    875887                        }
    876888                        else{
    … …  
    881893                else{ //this is an ice shelf: hydrostatic equilibrium*/
    882894                        if(hydroadjustment==AbsoluteEnum){
    883                                 newsurface[i] = newthickness[i]*(1-rho_ice/rho_water)+sealevel[i];
    884                                 newbase[i]     = newthickness[i]*(-rho_ice/rho_water)+sealevel[i];
     895                                newsurface[i] = newthickness[i]*(1.-rho_ice/rho_water)+sealevel[i];
     896                                newbase[i]    = newthickness[i]*(-rho_ice/rho_water)+sealevel[i];
    885897                        }
    886898                        else if(hydroadjustment==IncrementalEnum){
    887899                                newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i])+sealevel[i]; //surface = oldsurface + (1-di) * dH
    888                                 newbase[i]     = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base               = oldbed + di * dH
     900                                newbase[i]    = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base               = oldbed + di * dH
    889901                        }
    890902                        else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet");
    … …  
    893905
    894906        /*Add input to the element: */
    895         element->AddBasalInput(ThicknessEnum,newthickness,P1Enum);
    896907        element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
    897908        element->AddBasalInput(BaseEnum,newbase,P1Enum);
    … …  
    911922        xDelete<IssmDouble>(sealevel);
    912923        xDelete<IssmDouble>(bed);
    913 
    914924        xDelete<int>(doflist);
    915925        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r24089 r24091  
    2020/* SEMIC prototype {{{*/
    2121extern "C" void run_semic_(double *sf_in, double *rf_in, double *swd_in, double *lwd_in, double *wind_in, double *sp_in, double *rhoa_in,
    22                                 double *qq_in, double *tt_in, double *tsurf_out, double *smb_out, double *saccu_out, double *smelt_out);
     22                        double *qq_in, double *tt_in, double *tsurf_out, double *smb_out, double *saccu_out, double *smelt_out);
    2323#endif
    2424// _HAVE_SEMIC_
    … …  
    161161
    162162        /* Fetch number of nodes and allocate output*/
    163    int numnodes = this->GetNumberOfNodes();
    164    IssmDouble* newD = xNew<IssmDouble>(numnodes);
     163        int numnodes = this->GetNumberOfNodes();
     164        IssmDouble* newD = xNew<IssmDouble>(numnodes);
    165165
    166166        /* Retrieve domain-dependent inputs */
    167167        Input* n_input=this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
    168    Input* damage_input = NULL;
    169    Input* B_input = NULL;
     168        Input* damage_input = NULL;
     169        Input* B_input = NULL;
    170170        int domaintype;
    171    parameters->FindParam(&domaintype,DomainTypeEnum);
     171        parameters->FindParam(&domaintype,DomainTypeEnum);
    172172        if(domaintype==Domain2DhorizontalEnum){
    173       damage_input = this->GetInput(DamageDbarEnum);  _assert_(damage_input);
    174       B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
    175    }
    176    else{
    177       damage_input = this->GetInput(DamageDEnum);   _assert_(damage_input);
    178       B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
    179    }
     173                damage_input = this->GetInput(DamageDbarEnum);  _assert_(damage_input);
     174                B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
     175        }
     176        else{
     177                damage_input = this->GetInput(DamageDEnum);   _assert_(damage_input);
     178                B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
     179        }
    180180
    181181        /* Start looping on the number of nodes: */
    … …  
    198198
    199199                B_input->GetInputValue(&B,gauss);
    200       n_input->GetInputValue(&n,gauss);
    201       damage_input->GetInputValue(&D,gauss);
     200                n_input->GetInputValue(&n,gauss);
     201                damage_input->GetInputValue(&D,gauss);
    202202
    203203                /* Compute threshold strain rate from threshold stress */
    … …  
    268268
    269269                if(dim==2){
    270                          /* epsilon=[exx,eyy,exy];*/
     270                        /* epsilon=[exx,eyy,exy];*/
    271271                        eps_xx[iv]=epsilon[0];
    272272                        eps_yy[iv]=epsilon[1];
    … …  
    705705        }
    706706        switch(this->ObjectEnum()){
    707         case TriaEnum:
    708                 this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb[0],P1Enum));
    709                 this->inputs->AddInput(new TriaInput(SmbRunoffEnum,&runoff[0],P1Enum));
    710                 break;
    711         case PentaEnum:
    712                 this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&smb[0],P1Enum));
    713                 this->inputs->AddInput(new PentaInput(SmbRunoffEnum,&runoff[0],P1Enum));
    714                 this->InputExtrude(SmbMassBalanceEnum,-1);
    715                 this->InputExtrude(SmbRunoffEnum,-1);
    716                 break;
    717         case TetraEnum:
    718                 this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&smb[0],P1Enum));
    719                 this->inputs->AddInput(new TetraInput(SmbRunoffEnum,&runoff[0],P1Enum));
    720                 this->InputExtrude(SmbMassBalanceEnum,-1);
    721                 this->InputExtrude(SmbRunoffEnum,-1);
    722                 break;
    723         default: _error_("Not implemented yet");
     707                case TriaEnum:
     708                        this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb[0],P1Enum));
     709                        this->inputs->AddInput(new TriaInput(SmbRunoffEnum,&runoff[0],P1Enum));
     710                        break;
     711                case PentaEnum:
     712                        this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&smb[0],P1Enum));
     713                        this->inputs->AddInput(new PentaInput(SmbRunoffEnum,&runoff[0],P1Enum));
     714                        this->InputExtrude(SmbMassBalanceEnum,-1);
     715                        this->InputExtrude(SmbRunoffEnum,-1);
     716                        break;
     717                case TetraEnum:
     718                        this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&smb[0],P1Enum));
     719                        this->inputs->AddInput(new TetraInput(SmbRunoffEnum,&runoff[0],P1Enum));
     720                        this->InputExtrude(SmbMassBalanceEnum,-1);
     721                        this->InputExtrude(SmbRunoffEnum,-1);
     722                        break;
     723                default: _error_("Not implemented yet");
    724724        }
    725725        /* this->AddInput(SmbMassBalanceEnum,smb,P1Enum); */
    … …  
    858858        IssmDouble eps_eff;
    859859
    860          if(dim==2){
    861                  /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
    862                  this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
    863                  eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
    864          }
    865          else{
    866                  /* eps_eff^2 = 1/2 exx^2*/
    867                  this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
    868                  eps_eff = sqrt(epsilon1d*epsilon1d/2.);
    869          }
     860        if(dim==2){
     861                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
     862                this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
     863                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
     864        }
     865        else{
     866                /* eps_eff^2 = 1/2 exx^2*/
     867                this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
     868                eps_eff = sqrt(epsilon1d*epsilon1d/2.);
     869        }
    870870
    871871        /*Get viscosity*/
    … …  
    13191319}/*}}}*/
    13201320Node*      Element::GetNode(int nodeindex){/*{{{*/
    1321    _assert_(nodeindex>=0);
    1322    _assert_(nodeindex<this->GetNumberOfNodes(this->element_type));
    1323    return this->nodes[nodeindex];
     1321        _assert_(nodeindex>=0);
     1322        _assert_(nodeindex<this->GetNumberOfNodes(this->element_type));
     1323        return this->nodes[nodeindex];
    13241324}/*}}}*/
    13251325int        Element::GetNodeIndex(Node* node){/*{{{*/
    … …  
    15111511
    15121512        switch(type){
    1513         case VertexPIdEnum:
    1514                 doflist = xNew<int>(NUM_VERTICES);
    1515                 values = xNew<IssmDouble>(NUM_VERTICES);
    1516                 /*Fill in values*/
    1517                 this->GetVerticesPidList(doflist);
    1518                 this->GetInputListOnVerticesAtTime(values,input_enum,time);
    1519                 vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
    1520                 break;
    1521         case VertexSIdEnum:
    1522                 doflist = xNew<int>(NUM_VERTICES);
    1523                 values = xNew<IssmDouble>(NUM_VERTICES);
    1524                 /*Fill in values*/
    1525                 this->GetVerticesSidList(doflist);
    1526                 this->GetInputListOnVerticesAtTime(values,input_enum,time);
    1527                 vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
    1528                 break;
    1529         default:
    1530                 _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
     1513                case VertexPIdEnum:
     1514                        doflist = xNew<int>(NUM_VERTICES);
     1515                        values = xNew<IssmDouble>(NUM_VERTICES);
     1516                        /*Fill in values*/
     1517                        this->GetVerticesPidList(doflist);
     1518                        this->GetInputListOnVerticesAtTime(values,input_enum,time);
     1519                        vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
     1520                        break;
     1521                case VertexSIdEnum:
     1522                        doflist = xNew<int>(NUM_VERTICES);
     1523                        values = xNew<IssmDouble>(NUM_VERTICES);
     1524                        /*Fill in values*/
     1525                        this->GetVerticesSidList(doflist);
     1526                        this->GetInputListOnVerticesAtTime(values,input_enum,time);
     1527                        vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
     1528                        break;
     1529                default:
     1530                        _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    15311531        }
    15321532
    … …  
    17411741void       Element::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/
    17421742
    1743     /*Intermediaries*/
    1744     int i,t;
    1745 
    1746     /*Branch on type of vector: nodal or elementary: */
    1747     if(vector_type==1){ //nodal vector
    1748 
    1749         const int NUM_VERTICES = this->GetNumberOfVertices();
    1750 
    1751         int        *vertexids   = xNew<int>(NUM_VERTICES);
    1752         IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
    1753 
    1754         /*Recover vertices ids needed to initialize inputs*/
    1755         _assert_(iomodel->elements);
    1756         for(i=0;i<NUM_VERTICES;i++){
    1757             vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    1758         }
    1759 
    1760         /*Are we in transient or static? */
    1761                   if(M==1){
    1762                           values[0]=vector[0];
    1763                           this->AddInput(vector_enum,values,P0Enum);
    1764                   }
    1765                   else if(M==iomodel->numberofvertices){
    1766             for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
    1767             this->AddInput(vector_enum,values,P1Enum);
    1768         }
    1769         else if(M==iomodel->numberofvertices+1){
    1770             /*create transient input: */
    1771             IssmDouble* times = xNew<IssmDouble>(N);
    1772             for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1773             TransientInput* transientinput=new TransientInput(vector_enum,times,N);
    1774             for(t=0;t<N;t++){
    1775                 for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    1776                 switch(this->ObjectEnum()){
    1777                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
    1778                     case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break;
    1779                     case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break;
    1780                     default: _error_("Not implemented yet");
    1781                 }
    1782             }
    1783             this->inputs->AddInput(transientinput);
    1784             xDelete<IssmDouble>(times);
    1785         }
    1786         else if(M==iomodel->numberofelements){
    1787 
    1788                           /*This is a Patch!*/
    1789                           xDelete<IssmDouble>(values);
    1790                           values = xNew<IssmDouble>(N);
    1791                           for(int j=0;j<N;j++) values[j]=vector[this->Sid()*N+j];
    1792 
    1793                           if     (N==this->GetNumberOfNodes(P1Enum)   ) this->AddInput(vector_enum,values,P1Enum);
    1794                           else if(N==this->GetNumberOfNodes(P0Enum)   ) this->AddInput(vector_enum,values,P0Enum);
    1795                           else if(N==this->GetNumberOfNodes(P1xP2Enum)) this->AddInput(vector_enum,values,P1xP2Enum);
    1796                           else if(N==this->GetNumberOfNodes(P1xP3Enum)) this->AddInput(vector_enum,values,P1xP3Enum);
    1797                           else _error_("Patch interpolation not supported yet");
    1798 
    1799                   }
    1800                   else{
    1801                           _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    1802                   }
    1803 
    1804         xDelete<IssmDouble>(values);
    1805         xDelete<int>(vertexids);
    1806     }
    1807     else if(vector_type==2){ //element vector
    1808 
    1809         IssmDouble value;
    1810 
    1811         /*Are we in transient or static? */
    1812         if(M==iomodel->numberofelements){
    1813             if (code==5){ //boolean
    1814                 this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->Sid()])));
    1815             }
    1816             else if (code==6){ //integer
    1817                 this->inputs->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->Sid()])));
    1818             }
    1819             else if (code==7){ //IssmDouble
    1820                 this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]));
    1821             }
    1822             else _error_("could not recognize nature of vector from code " << code);
    1823         }
    1824         else if(M==iomodel->numberofelements+1){
    1825             /*create transient input: */
    1826             IssmDouble* times = xNew<IssmDouble>(N);
    1827             for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1828             TransientInput* transientinput=new TransientInput(vector_enum,times,N);
    1829             TriaInput* bof=NULL;
    1830             for(t=0;t<N;t++){
    1831                 value=vector[N*this->Sid()+t];
    1832                 switch(this->ObjectEnum()){
    1833                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break;
    1834                     case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break;
    1835                     case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break;
    1836                     default: _error_("Not implemented yet");
    1837                 }
    1838             }
    1839             this->inputs->AddInput(transientinput);
    1840             xDelete<IssmDouble>(times);
    1841         }
    1842         else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    1843     }
    1844     else if(vector_type==3){ //element vector
    1845 
    1846         /*For right now we are static */
    1847         if(M==iomodel->numberofelements){
    1848             /*create transient input: */
    1849             IssmDouble* layers = xNewZeroInit<IssmDouble>(N);;
    1850             for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
    1851             DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N);
    1852             this->inputs->AddInput(arrayinput);
    1853             xDelete<IssmDouble>(layers);
    1854         }
    1855         else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    1856     }
    1857     else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
     1743        /*Intermediaries*/
     1744        int i,t;
     1745
     1746        /*Branch on type of vector: nodal or elementary: */
     1747        if(vector_type==1){ //nodal vector
     1748
     1749                const int NUM_VERTICES = this->GetNumberOfVertices();
     1750
     1751                int        *vertexids   = xNew<int>(NUM_VERTICES);
     1752                IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
     1753
     1754                /*Recover vertices ids needed to initialize inputs*/
     1755                _assert_(iomodel->elements);
     1756                for(i=0;i<NUM_VERTICES;i++){
     1757                        vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
     1758                }
     1759
     1760                /*Are we in transient or static? */
     1761                if(M==1){
     1762                        values[0]=vector[0];
     1763                        this->AddInput(vector_enum,values,P0Enum);
     1764                }
     1765                else if(M==iomodel->numberofvertices){
     1766                        for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
     1767                        this->AddInput(vector_enum,values,P1Enum);
     1768                }
     1769                else if(M==iomodel->numberofvertices+1){
     1770                        /*create transient input: */
     1771                        IssmDouble* times = xNew<IssmDouble>(N);
     1772                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1773                        TransientInput* transientinput=new TransientInput(vector_enum,times,N);
     1774                        for(t=0;t<N;t++){
     1775                                for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
     1776                                switch(this->ObjectEnum()){
     1777                                        case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
     1778                                        case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break;
     1779                                        case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break;
     1780                                        default: _error_("Not implemented yet");
     1781                                }
     1782                        }
     1783                        this->inputs->AddInput(transientinput);
     1784                        xDelete<IssmDouble>(times);
     1785                }
     1786                else if(M==iomodel->numberofelements){
     1787
     1788                        /*This is a Patch!*/
     1789                        xDelete<IssmDouble>(values);
     1790                        values = xNew<IssmDouble>(N);
     1791                        for(int j=0;j<N;j++) values[j]=vector[this->Sid()*N+j];
     1792
     1793                        if     (N==this->GetNumberOfNodes(P1Enum)   ) this->AddInput(vector_enum,values,P1Enum);
     1794                        else if(N==this->GetNumberOfNodes(P0Enum)   ) this->AddInput(vector_enum,values,P0Enum);
     1795                        else if(N==this->GetNumberOfNodes(P1xP2Enum)) this->AddInput(vector_enum,values,P1xP2Enum);
     1796                        else if(N==this->GetNumberOfNodes(P1xP3Enum)) this->AddInput(vector_enum,values,P1xP3Enum);
     1797                        else _error_("Patch interpolation not supported yet");
     1798
     1799                }
     1800                else{
     1801                        _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     1802                }
     1803
     1804                xDelete<IssmDouble>(values);
     1805                xDelete<int>(vertexids);
     1806        }
     1807        else if(vector_type==2){ //element vector
     1808
     1809                IssmDouble value;
     1810
     1811                /*Are we in transient or static? */
     1812                if(M==iomodel->numberofelements){
     1813                        if (code==5){ //boolean
     1814                                this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->Sid()])));
     1815                        }
     1816                        else if (code==6){ //integer
     1817                                this->inputs->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->Sid()])));
     1818                        }
     1819                        else if (code==7){ //IssmDouble
     1820                                this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]));
     1821                        }
     1822                        else _error_("could not recognize nature of vector from code " << code);
     1823                }
     1824                else if(M==iomodel->numberofelements+1){
     1825                        /*create transient input: */
     1826                        IssmDouble* times = xNew<IssmDouble>(N);
     1827                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1828                        TransientInput* transientinput=new TransientInput(vector_enum,times,N);
     1829                        TriaInput* bof=NULL;
     1830                        for(t=0;t<N;t++){
     1831                                value=vector[N*this->Sid()+t];
     1832                                switch(this->ObjectEnum()){
     1833                                        case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break;
     1834                                        case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break;
     1835                                        case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break;
     1836                                        default: _error_("Not implemented yet");
     1837                                }
     1838                        }
     1839                        this->inputs->AddInput(transientinput);
     1840                        xDelete<IssmDouble>(times);
     1841                }
     1842                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     1843        }
     1844        else if(vector_type==3){ //element vector
     1845
     1846                /*For right now we are static */
     1847                if(M==iomodel->numberofelements){
     1848                        /*create transient input: */
     1849                        IssmDouble* layers = xNewZeroInit<IssmDouble>(N);;
     1850                        for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
     1851                        DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N);
     1852                        this->inputs->AddInput(arrayinput);
     1853                        xDelete<IssmDouble>(layers);
     1854                }
     1855                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     1856        }
     1857        else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
    18581858}
    18591859/*}}}*/
    … …  
    18941894
    18951895        else if(M==iomodel->numberofvertices+1){
    1896             /*create transient input: */
    1897             IssmDouble* times = xNew<IssmDouble>(N);
    1898             for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1899                                 /*Create the three transient inputs for the control input*/
    1900             TransientInput* values_input=new TransientInput(input_enum,times,N);
    1901                                 TransientInput* mins_input = new TransientInput(ControlInputMinsEnum,times,N);
    1902                                 TransientInput* maxs_input = new TransientInput(ControlInputMaxsEnum,times,N);
    1903                            TransientInput* grad_input = new TransientInput(ControlInputGradEnum);
    1904                                 for(int t=0;t<N;t++){
    1905                 for(int i=0;i<NUM_VERTICES;i++){
    1906                                                 values[i]=vector[N*(vertexids[i]-1)+t];
    1907                                                 values_min[i] = min_vector[N*(vertexids[i]-1)+t];
    1908                                                 values_max[i] = max_vector[N*(vertexids[i]-1)+t];
    1909                                          }
    1910                                         switch(this->ObjectEnum()){
    1911                     case TriaEnum:
    1912                                                                         values_input->AddTimeInput(new TriaInput(input_enum,values,P1Enum));
    1913                                                                         mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum));
    1914                                                                         maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum));
    1915                                                                 break;
    1916                     case PentaEnum:
    1917                                                                         values_input->AddTimeInput(new PentaInput(input_enum,values,P1Enum));
    1918                                                                         mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum));
    1919                                                                         maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum));
    1920                                                                         break;
    1921                     case TetraEnum:
    1922                                                                         values_input->AddTimeInput(new TetraInput(input_enum,values,P1Enum));
    1923                                                                         mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum));
    1924                                                                         maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum));
    1925                                                                         break;
    1926                     default: _error_("Not implemented yet");
    1927                 }
    1928             }
    1929             this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id));
    1930             xDelete<IssmDouble>(times);
    1931         }
    1932                   else _error_("not currently supported type of M and N attempted");
     1896                /*create transient input: */
     1897                IssmDouble* times = xNew<IssmDouble>(N);
     1898                for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1899                /*Create the three transient inputs for the control input*/
     1900                TransientInput* values_input=new TransientInput(input_enum,times,N);
     1901                TransientInput* mins_input = new TransientInput(ControlInputMinsEnum,times,N);
     1902                TransientInput* maxs_input = new TransientInput(ControlInputMaxsEnum,times,N);
     1903                TransientInput* grad_input = new TransientInput(ControlInputGradEnum);
     1904                for(int t=0;t<N;t++){
     1905                        for(int i=0;i<NUM_VERTICES;i++){
     1906                                values[i]=vector[N*(vertexids[i]-1)+t];
     1907                                values_min[i] = min_vector[N*(vertexids[i]-1)+t];
     1908                                values_max[i] = max_vector[N*(vertexids[i]-1)+t];
     1909                        }
     1910                        switch(this->ObjectEnum()){
     1911                                case TriaEnum:
     1912                                        values_input->AddTimeInput(new TriaInput(input_enum,values,P1Enum));
     1913                                        mins_input->AddTimeInput(new TriaInput(ControlInputMinsEnum,values_min,P1Enum));
     1914                                        maxs_input->AddTimeInput(new TriaInput(ControlInputMaxsEnum,values_max,P1Enum));
     1915                                        break;
     1916                                case PentaEnum:
     1917                                        values_input->AddTimeInput(new PentaInput(input_enum,values,P1Enum));
     1918                                        mins_input->AddTimeInput(new PentaInput(ControlInputMinsEnum,values_min,P1Enum));
     1919                                        maxs_input->AddTimeInput(new PentaInput(ControlInputMaxsEnum,values_max,P1Enum));
     1920                                        break;
     1921                                case TetraEnum:
     1922                                        values_input->AddTimeInput(new TetraInput(input_enum,values,P1Enum));
     1923                                        mins_input->AddTimeInput(new TetraInput(ControlInputMinsEnum,values_min,P1Enum));
     1924                                        maxs_input->AddTimeInput(new TetraInput(ControlInputMaxsEnum,values_max,P1Enum));
     1925                                        break;
     1926                                default: _error_("Not implemented yet");
     1927                        }
     1928                }
     1929                this->inputs->AddInput(new ControlInput(input_enum,TransientInputEnum,values_input,mins_input,maxs_input,grad_input,P1Enum,id));
     1930                xDelete<IssmDouble>(times);
     1931        }
     1932        else _error_("not currently supported type of M and N attempted");
    19331933
    19341934        /*clean up*/
    … …  
    19381938/*}}}*/
    19391939void       Element::DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int input_enum,int code,int input_id){/*{{{*/
    1940     /*enum_type: the name of the DatasetInput (eg Outputdefinition1)
    1941           * vector: information being stored (eg observations)
    1942           * vector_type: is if by element or by vertex
    1943           * input_enum: is the name of the vector being stored
    1944           * code: what type of data is in the vector (booleans, ints, doubles)
    1945           */
    1946 
    1947          /*Intermediaries*/
    1948     int                                 i,t;
    1949          DatasetInput*          datasetinput = NULL;
    1950 
    1951          /*Get input if it already exists*/
    1952          Input*  tempinput = GetInput(enum_type);
    1953          if(tempinput){
    1954                  /*Cast it to a Datasetinput*/
    1955                  if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
    1956                  datasetinput = (DatasetInput*)tempinput;
    1957          }
    1958          else{
    1959                  datasetinput=new DatasetInput(enum_type);
    1960            this->inputs->AddInput(datasetinput);
    1961         }
    1962 
    1963     /*Branch on type of vector: nodal or elementary: */
    1964     if(vector_type==1){ //nodal vector
    1965 
    1966         const int NUM_VERTICES = this->GetNumberOfVertices();
    1967 
    1968         int        *vertexids   = xNew<int>(NUM_VERTICES);
    1969         IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
    1970 
    1971         /*Recover vertices ids needed to initialize inputs*/
    1972         _assert_(iomodel->elements);
    1973         for(i=0;i<NUM_VERTICES;i++){
    1974             vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    1975         }
    1976 
    1977         /*Are we in transient or static? */
    1978                   if(M==1){
    1979                           values[0]=vector[0];
     1940        /*enum_type: the name of the DatasetInput (eg Outputdefinition1)
     1941         * vector: information being stored (eg observations)
     1942         * vector_type: is if by element or by vertex
     1943         * input_enum: is the name of the vector being stored
     1944         * code: what type of data is in the vector (booleans, ints, doubles)
     1945         */
     1946
     1947        /*Intermediaries*/
     1948        int                                     i,t;
     1949        DatasetInput*           datasetinput = NULL;
     1950
     1951        /*Get input if it already exists*/
     1952        Input*  tempinput = GetInput(enum_type);
     1953        if(tempinput){
     1954                /*Cast it to a Datasetinput*/
     1955                if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
     1956                datasetinput = (DatasetInput*)tempinput;
     1957        }
     1958        else{
     1959                datasetinput=new DatasetInput(enum_type);
     1960                this->inputs->AddInput(datasetinput);
     1961        }
     1962
     1963        /*Branch on type of vector: nodal or elementary: */
     1964        if(vector_type==1){ //nodal vector
     1965
     1966                const int NUM_VERTICES = this->GetNumberOfVertices();
     1967
     1968                int        *vertexids   = xNew<int>(NUM_VERTICES);
     1969                IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
     1970
     1971                /*Recover vertices ids needed to initialize inputs*/
     1972                _assert_(iomodel->elements);
     1973                for(i=0;i<NUM_VERTICES;i++){
     1974                        vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
     1975                }
     1976
     1977                /*Are we in transient or static? */
     1978                if(M==1){
     1979                        values[0]=vector[0];
     1980                        switch(this->ObjectEnum()){
     1981                                case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break;
     1982                                case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break;
     1983                                case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break;
     1984                                default: _error_("Not implemented yet");
     1985                        }
     1986                }
     1987                else if(M==iomodel->numberofvertices){
     1988                        for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
     1989                        switch(this->ObjectEnum()){
     1990                                case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
     1991                                case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break;
     1992                                case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break;
     1993                                default: _error_("Not implemented yet");
     1994                        }  }
     1995                else if(M==iomodel->numberofvertices+1){
     1996                        /*create transient input: */
     1997                        IssmDouble* times = xNew<IssmDouble>(N);
     1998                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1999                        TransientInput* transientinput=new TransientInput(input_enum,times,N);
     2000                        for(t=0;t<N;t++){
     2001                                for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    19802002                                switch(this->ObjectEnum()){
    1981                     case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break;
    1982                     case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break;
    1983                     case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break;
    1984                     default: _error_("Not implemented yet");
     2003                                        case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break;
     2004                                        case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); break;
     2005                                        case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); break;
     2006                                        default: _error_("Not implemented yet");
    19852007                                }
    1986                   }
    1987                   else if(M==iomodel->numberofvertices){
    1988             for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
     2008                        }
     2009                        datasetinput->AddInput(transientinput,input_id);
     2010                        xDelete<IssmDouble>(times);
     2011                }
     2012                else if(M==iomodel->numberofelements){
     2013
     2014                        /*This is a Patch!*/
     2015                        xDelete<IssmDouble>(values);
     2016                        values = xNew<IssmDouble>(N);
     2017                        for(int j=0;j<N;j++) values[j]=vector[this->Sid()*N+j];
     2018
     2019                        if     (N==this->GetNumberOfNodes(P1Enum)   ){
    19892020                                switch(this->ObjectEnum()){
    1990                     case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
    1991                     case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break;
    1992                     case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break;
    1993                     default: _error_("Not implemented yet");
    1994                                 }  }
    1995         else if(M==iomodel->numberofvertices+1){
    1996             /*create transient input: */
    1997             IssmDouble* times = xNew<IssmDouble>(N);
    1998             for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1999             TransientInput* transientinput=new TransientInput(input_enum,times,N);
    2000             for(t=0;t<N;t++){
    2001                 for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    2002                 switch(this->ObjectEnum()){
    2003                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break;
    2004                     case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); break;
    2005                     case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); break;
    2006                     default: _error_("Not implemented yet");
    2007                 }
    2008             }
    2009             datasetinput->AddInput(transientinput,input_id);
    2010             xDelete<IssmDouble>(times);
    2011         }
    2012         else if(M==iomodel->numberofelements){
    2013 
    2014                           /*This is a Patch!*/
    2015                           xDelete<IssmDouble>(values);
    2016                           values = xNew<IssmDouble>(N);
    2017                           for(int j=0;j<N;j++) values[j]=vector[this->Sid()*N+j];
    2018 
    2019                           if     (N==this->GetNumberOfNodes(P1Enum)   ){
    2020                                   switch(this->ObjectEnum()){
    2021                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
    2022                                           case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break;
    2023                                           case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break;
    2024                                           default: _error_("Not implemented yet");
    2025                                   }
    2026                           }
    2027                           else if(N==this->GetNumberOfNodes(P0Enum)   ){
    2028                                   switch(this->ObjectEnum()){
    2029                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break;
    2030                                           case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break;
    2031                                           case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break;
    2032                                           default: _error_("Not implemented yet");
    2033                                   }
    2034                           }
    2035                           else if(N==this->GetNumberOfNodes(P1xP2Enum)){
    2036                                   switch(this->ObjectEnum()){
    2037                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1xP2Enum),input_id); break;
    2038                                           case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP2Enum),input_id); break;
    2039                                           case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP2Enum),input_id); break;
    2040                                           default: _error_("Not implemented yet");
    2041                                   }
    2042                           }
    2043                           else if(N==this->GetNumberOfNodes(P1xP3Enum)) {
    2044                                  switch(this->ObjectEnum()){
    2045                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1xP3Enum),input_id); break;
    2046                                           case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP3Enum),input_id); break;
    2047                                           case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP3Enum),input_id); break;
    2048                                           default: _error_("Not implemented yet");
    2049                                   }
    2050                           }
    2051                           else _error_("Patch interpolation not supported yet");
    2052 
    2053                   }
    2054                   else{
    2055                           _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
    2056                   }
    2057 
    2058         xDelete<IssmDouble>(values);
    2059         xDelete<int>(vertexids);
    2060     }
    2061     else if(vector_type==2){ //element vector
    2062 
    2063         IssmDouble value;
    2064 
    2065         /*Are we in transient or static? */
    2066         if(M==iomodel->numberofelements){
    2067             if (code==5){ //boolean
    2068                 datasetinput->AddInput(new BoolInput(input_enum,reCast<bool>(vector[this->Sid()])),input_id);
    2069             }
    2070             else if (code==6){ //integer
    2071                 datasetinput->AddInput(new IntInput(input_enum,reCast<int>(vector[this->Sid()])),input_id);
    2072             }
    2073             else if (code==7){ //IssmDouble
    2074                 datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id);
    2075             }
    2076             else _error_("could not recognize nature of vector from code " << code);
    2077         }
    2078         else if(M==iomodel->numberofelements+1){
    2079             /*create transient input: */
    2080             IssmDouble* times = xNew<IssmDouble>(N);
    2081             for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    2082             TransientInput* transientinput=new TransientInput(input_enum,times,N);
    2083             TriaInput* bof=NULL;
    2084             for(t=0;t<N;t++){
    2085                 value=vector[N*this->Sid()+t];
    2086                 switch(this->ObjectEnum()){
    2087                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break;
    2088                     case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break;
    2089                     case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break;
    2090                     default: _error_("Not implemented yet");
    2091                 }
    2092             }
    2093             datasetinput->AddInput(transientinput,input_id);
    2094             xDelete<IssmDouble>(times);
    2095         }
    2096         else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
    2097     }
    2098     else if(vector_type==3){ //element vector
    2099 
    2100         /*For right now we are static */
    2101         if(M==iomodel->numberofelements){
    2102             /*create transient input: */
    2103             IssmDouble* layers = xNewZeroInit<IssmDouble>(N);;
    2104             for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
    2105             DoubleArrayInput* arrayinput=new DoubleArrayInput(input_enum,layers,N);
    2106             datasetinput->AddInput(arrayinput,input_id);
    2107             xDelete<IssmDouble>(layers);
    2108         }
    2109         else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
    2110     }
    2111     else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
     2021                                        case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
     2022                                        case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break;
     2023                                        case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break;
     2024                                        default: _error_("Not implemented yet");
     2025                                }
     2026                        }
     2027                        else if(N==this->GetNumberOfNodes(P0Enum)   ){
     2028                                switch(this->ObjectEnum()){
     2029                                        case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break;
     2030                                        case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break;
     2031                                        case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break;
     2032                                        default: _error_("Not implemented yet");
     2033                                }
     2034                        }
     2035                        else if(N==this->GetNumberOfNodes(P1xP2Enum)){
     2036                                switch(this->ObjectEnum()){
     2037                                        case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1xP2Enum),input_id); break;
     2038                                        case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP2Enum),input_id); break;
     2039                                        case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP2Enum),input_id); break;
     2040                                        default: _error_("Not implemented yet");
     2041                                }
     2042                        }
     2043                        else if(N==this->GetNumberOfNodes(P1xP3Enum)) {
     2044                                switch(this->ObjectEnum()){
     2045                                        case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1xP3Enum),input_id); break;
     2046                                        case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP3Enum),input_id); break;
     2047                                        case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP3Enum),input_id); break;
     2048                                        default: _error_("Not implemented yet");
     2049                                }
     2050                        }
     2051                        else _error_("Patch interpolation not supported yet");
     2052
     2053                }
     2054                else{
     2055                        _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
     2056                }
     2057
     2058                xDelete<IssmDouble>(values);
     2059                xDelete<int>(vertexids);
     2060        }
     2061        else if(vector_type==2){ //element vector
     2062
     2063                IssmDouble value;
     2064
     2065                /*Are we in transient or static? */
     2066                if(M==iomodel->numberofelements){
     2067                        if (code==5){ //boolean
     2068                                datasetinput->AddInput(new BoolInput(input_enum,reCast<bool>(vector[this->Sid()])),input_id);
     2069                        }
     2070                        else if (code==6){ //integer
     2071                                datasetinput->AddInput(new IntInput(input_enum,reCast<int>(vector[this->Sid()])),input_id);
     2072                        }
     2073                        else if (code==7){ //IssmDouble
     2074                                datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id);
     2075                        }
     2076                        else _error_("could not recognize nature of vector from code " << code);
     2077                }
     2078                else if(M==iomodel->numberofelements+1){
     2079                        /*create transient input: */
     2080                        IssmDouble* times = xNew<IssmDouble>(N);
     2081                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     2082                        TransientInput* transientinput=new TransientInput(input_enum,times,N);
     2083                        TriaInput* bof=NULL;
     2084                        for(t=0;t<N;t++){
     2085                                value=vector[N*this->Sid()+t];
     2086                                switch(this->ObjectEnum()){
     2087                                        case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break;
     2088                                        case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break;
     2089                                        case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break;
     2090                                        default: _error_("Not implemented yet");
     2091                                }
     2092                        }
     2093                        datasetinput->AddInput(transientinput,input_id);
     2094                        xDelete<IssmDouble>(times);
     2095                }
     2096                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
     2097        }
     2098        else if(vector_type==3){ //element vector
     2099
     2100                /*For right now we are static */
     2101                if(M==iomodel->numberofelements){
     2102                        /*create transient input: */
     2103                        IssmDouble* layers = xNewZeroInit<IssmDouble>(N);;
     2104                        for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
     2105                        DoubleArrayInput* arrayinput=new DoubleArrayInput(input_enum,layers,N);
     2106                        datasetinput->AddInput(arrayinput,input_id);
     2107                        xDelete<IssmDouble>(layers);
     2108                }
     2109                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
     2110        }
     2111        else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
    21122112}
    21132113/*}}}*/
    … …  
    32523252        /*Some intputs need to be computed, even if they are already in inputs, they might not be up to date!*/
    32533253        switch(output_enum){
    3254         case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break;
    3255         case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break;
     3254                case ViscousHeatingEnum: this->ViscousHeatingCreateInput(); break;
     3255                case StressMaxPrincipalEnum: this->StressMaxPrincipalCreateInput(); break;
    32563256                case StressTensorxxEnum:
    32573257                case StressTensorxyEnum:
    … …  
    32663266                case StrainRateyzEnum:
    32673267                case StrainRatezzEnum:
    3268           case StrainRateeffectiveEnum: this->ComputeStrainRate(); break;
     3268                case StrainRateeffectiveEnum: this->ComputeStrainRate(); break;
    32693269                case DeviatoricStressxxEnum:
    32703270                case DeviatoricStressxyEnum:
    … …  
    32873287                case CalvingrateyEnum:
    32883288                case CalvingCalvingrateEnum:
    3289                         this->StrainRateparallel();
    3290                         this->StrainRateperpendicular();
    3291                         int calvinglaw;
    3292                         this->FindParam(&calvinglaw,CalvingLawEnum);
    3293                         switch(calvinglaw){
    3294                                 case DefaultCalvingEnum:
    3295                                         //do nothing
    3296                                         break;
    3297                                 case CalvingLevermannEnum:
    3298                                         this->CalvingRateLevermann();
    3299                                         break;
    3300                                 case CalvingVonmisesEnum:
    3301                                 case CalvingDev2Enum:
    3302                                         this->CalvingRateVonmises();
    3303                                         break;
    3304                                 case CalvingCrevasseDepthEnum:
    3305                                         this->CalvingCrevasseDepth();
    3306                                         break;
    3307                                 default:
    3308                                         _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
    3309                         }
    3310                         break;
     3289                                                                                                  this->StrainRateparallel();
     3290                                                                                                  this->StrainRateperpendicular();
     3291                                                                                                  int calvinglaw;
     3292                                                                                                  this->FindParam(&calvinglaw,CalvingLawEnum);
     3293                                                                                                  switch(calvinglaw){
     3294                                                                                                          case DefaultCalvingEnum:
     3295                                                                                                                  //do nothing
     3296                                                                                                                  break;
     3297                                                                                                          case CalvingLevermannEnum:
     3298                                                                                                                  this->CalvingRateLevermann();
     3299                                                                                                                  break;
     3300                                                                                                          case CalvingVonmisesEnum:
     3301                                                                                                          case CalvingDev2Enum:
     3302                                                                                                                  this->CalvingRateVonmises();
     3303                                                                                                                  break;
     3304                                                                                                          case CalvingCrevasseDepthEnum:
     3305                                                                                                                  this->CalvingCrevasseDepth();
     3306                                                                                                                  break;
     3307                                                                                                          default:
     3308                                                                                                                  _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     3309                                                                                                  }
     3310                                                                                                  break;
    33113311                case CalvingFluxLevelsetEnum: this->CalvingFluxLevelset(); break;
    33123312                case CalvingMeltingFluxLevelsetEnum: this->CalvingMeltingFluxLevelset(); break;
    … …  
    33513351        switch(input->GetResultInterpolation()){
    33523352                case P0Enum:{
    3353                         IssmDouble  value;
    3354                         bool        bvalue;
    3355                         Input*      input = this->GetInput(output_enum); _assert_(input);
    3356                         switch(input->ObjectEnum()){
    3357                                 case DoubleInputEnum:
    3358                                         input->GetInputValue(&value);
    3359                                         break;
    3360                                 case BoolInputEnum:
    3361                                         input->GetInputValue(&bvalue);
    3362                                         value=reCast<IssmDouble>(bvalue);
    3363                                         break;
    3364                                 default:
    3365                                         Gauss* gauss = this->NewGauss();
    3366                                         input->GetInputValue(&value,gauss);
    3367                                         delete gauss;
    3368                         }
    3369                         vector->SetValue(this->Sid(),value,INS_VAL);
    3370                         break;
    3371                 }
     3353                                                        IssmDouble  value;
     3354                                                        bool        bvalue;
     3355                                                        Input*      input = this->GetInput(output_enum); _assert_(input);
     3356                                                        switch(input->ObjectEnum()){
     3357                                                                case DoubleInputEnum:
     3358                                                                        input->GetInputValue(&value);
     3359                                                                        break;
     3360                                                                case BoolInputEnum:
     3361                                                                        input->GetInputValue(&bvalue);
     3362                                                                        value=reCast<IssmDouble>(bvalue);
     3363                                                                        break;
     3364                                                                default:
     3365                                                                        Gauss* gauss = this->NewGauss();
     3366                                                                        input->GetInputValue(&value,gauss);
     3367                                                                        delete gauss;
     3368                                                        }
     3369                                                        vector->SetValue(this->Sid(),value,INS_VAL);
     3370                                                        break;
     3371                                                }
    33723372                case P1Enum:{
    3373                         const int NUM_VERTICES = this->GetNumberOfVertices();
    3374 
    3375                         IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
    3376                         int        *connectivity= xNew<int>(NUM_VERTICES);
    3377                         int        *sidlist     = xNew<int>(NUM_VERTICES);
    3378 
    3379                         this->GetVerticesSidList(sidlist);
    3380                         this->GetVerticesConnectivityList(connectivity);
    3381                         this->GetInputListOnVertices(values,output_enum);
    3382                         for(int i=0;i<NUM_VERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
    3383 
    3384                         vector->SetValues(NUM_VERTICES,sidlist,values,ADD_VAL);
    3385 
    3386                         xDelete<IssmDouble>(values);
    3387                         xDelete<int>(connectivity);
    3388                         xDelete<int>(sidlist);
    3389                         break;
    3390                 }
     3373                                                        const int NUM_VERTICES = this->GetNumberOfVertices();
     3374
     3375                                                        IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
     3376                                                        int        *connectivity= xNew<int>(NUM_VERTICES);
     3377                                                        int        *sidlist     = xNew<int>(NUM_VERTICES);
     3378
     3379                                                        this->GetVerticesSidList(sidlist);
     3380                                                        this->GetVerticesConnectivityList(connectivity);
     3381                                                        this->GetInputListOnVertices(values,output_enum);
     3382                                                        for(int i=0;i<NUM_VERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
     3383
     3384                                                        vector->SetValues(NUM_VERTICES,sidlist,values,ADD_VAL);
     3385
     3386                                                        xDelete<IssmDouble>(values);
     3387                                                        xDelete<int>(connectivity);
     3388                                                        xDelete<int>(sidlist);
     3389                                                        break;
     3390                                                }
    33913391                default:
    33923392                                         _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet");

  • issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp

    r24089 r24091  
    243243
    244244        /*Build B for this segment*/
    245         B[0] = +basis[index1];
    246         B[1] = +basis[index2];
    247         B[2] = -basis[index1];
    248         B[3] = -basis[index2];
     245        switch(finiteelement){
     246                case P0DGEnum:
     247                        B[0] = +basis[0];
     248                        B[1] = -basis[0];
     249                        break;
     250                case P1DGEnum:
     251                        B[0] = +basis[index1];
     252                        B[1] = +basis[index2];
     253                        B[2] = -basis[index1];
     254                        B[3] = -basis[index2];
     255                        break;
     256                default:
     257                        _error_("not supported yet");
     258        }
    249259
    250260        /*Clean-up*/
    … …  
    268278
    269279        /*Build B'*/
    270         Bprime[0] = basis[index1];
    271         Bprime[1] = basis[index2];
    272         Bprime[2] = basis[index1];
    273         Bprime[3] = basis[index2];
     280        /*Build B for this segment*/
     281        switch(finiteelement){
     282                case P0DGEnum:
     283                        Bprime[0] = basis[0];
     284                        Bprime[1] = basis[0];
     285                        break;
     286                case P1DGEnum:
     287                        Bprime[0] = basis[index1];
     288                        Bprime[1] = basis[index2];
     289                        Bprime[2] = basis[index1];
     290                        Bprime[3] = basis[index2];
     291                        break;
     292                default:
     293                        _error_("not supported yet");
     294        }
    274295
    275296        /*Clean-up*/
    … …  
    305326
    306327        switch(finiteelement){
     328                case P0Enum: case P0DGEnum:
     329                        basis[0]=triabasis[0];
     330                        xDelete<IssmDouble>(triabasis);
     331                        return;
    307332                case P1Enum: case P1DGEnum:
    308333                        basis[0]=triabasis[index1];

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r23998 r24091  
    9696int  TriaInput::GetResultInterpolation(void){/*{{{*/
    9797
    98         if(this->interpolation_type==P0Enum){
     98        if(this->interpolation_type==P0Enum || this->interpolation_type==P0DGEnum){
    9999                return P0Enum;
    100100        }

  • issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp

    r24090 r24091  
    6262   /*FIXME: hardcode element degree for now*/
    6363   this->flux_degree= P1DGEnum;
    64    //printf("-------------- file: Numericalflux.cpp line: %i\n",__LINE__);
    6564   //this->flux_degree= P0DGEnum;
    6665

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp

    r24088 r24091  
    512512        }
    513513        /*}}}*/
    514         if(finite_element==P1DGEnum && analysis!=UzawaPressureAnalysisEnum){/*Special case for DG...{{{*/
    515                 int node_list[4];
    516                 if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet");
    517                 CreateEdges(iomodel);
    518                 CreateFaces(iomodel);
    519                 for(int i=0;i<iomodel->numberoffaces;i++){
    520                         int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
    521                         int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
    522                         if(e2!=-2){
    523                                 if(epart[e1]!=epart[e2]){
    524                                         int i1=iomodel->faces[4*i+0];
    525                                         int i2=iomodel->faces[4*i+1];
    526                                         int pos=-1;
    527                                         for(int j=0;j<3;j++) if(iomodel->elements[3*e2+j]==i1) pos=j;
    528                                         if(     pos==0){ node_list[0] = e2*3+0; node_list[1] = e2*3+2;}
    529                                         else if(pos==1){ node_list[0] = e2*3+1; node_list[1] = e2*3+0;}
    530                                         else if(pos==2){ node_list[0] = e2*3+2; node_list[1] = e2*3+1;}
    531                                         else _error_("not supposed to happen");
    532                                         pos=-1;
    533                                         for(int j=0;j<3;j++) if(iomodel->elements[3*e1+j]==i1) pos=j;
    534                                         if(     pos==0){ node_list[2] = e1*3+0; node_list[3] = e1*3+1;}
    535                                         else if(pos==1){ node_list[2] = e1*3+1; node_list[3] = e1*3+2;}
    536                                         else if(pos==2){ node_list[2] = e1*3+2; node_list[3] = e1*3+0;}
    537                                         else _error_("not supposed to happen");
    538                                         for(int j=0;j<4;j++){
    539                                                 int  nid = node_list[j];
    540                                                 AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e1]);
    541                                                 AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e2]);
     514        if((finite_element==P0DGEnum || finite_element==P1DGEnum) && analysis!=UzawaPressureAnalysisEnum){/*Special case for DG...{{{*/
     515                if(finite_element==P1DGEnum){
     516                        int node_list[4];
     517                        if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet");
     518                        CreateEdges(iomodel);
     519                        CreateFaces(iomodel);
     520                        for(int i=0;i<iomodel->numberoffaces;i++){
     521                                int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     522                                int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
     523                                if(e2!=-2){
     524                                        if(epart[e1]!=epart[e2]){
     525                                                int i1=iomodel->faces[4*i+0];
     526                                                int i2=iomodel->faces[4*i+1];
     527                                                int pos=-1;
     528                                                for(int j=0;j<3;j++) if(iomodel->elements[3*e2+j]==i1) pos=j;
     529                                                if(     pos==0){ node_list[0] = e2*3+0; node_list[1] = e2*3+2;}
     530                                                else if(pos==1){ node_list[0] = e2*3+1; node_list[1] = e2*3+0;}
     531                                                else if(pos==2){ node_list[0] = e2*3+2; node_list[1] = e2*3+1;}
     532                                                else _error_("not supposed to happen");
     533                                                pos=-1;
     534                                                for(int j=0;j<3;j++) if(iomodel->elements[3*e1+j]==i1) pos=j;
     535                                                if(     pos==0){ node_list[2] = e1*3+0; node_list[3] = e1*3+1;}
     536                                                else if(pos==1){ node_list[2] = e1*3+1; node_list[3] = e1*3+2;}
     537                                                else if(pos==2){ node_list[2] = e1*3+2; node_list[3] = e1*3+0;}
     538                                                else _error_("not supposed to happen");
     539                                                for(int j=0;j<4;j++){
     540                                                        int  nid = node_list[j];
     541                                                        AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e1]);
     542                                                        AddNodeToRank(nodes_ranks,nodes_proc_count,nid,epart[e2]);
     543                                                }
    542544                                        }
    543545                                }
    544546                        }
     547                }
     548                else if(finite_element==P0DGEnum){
     549                        int node_list[2];
     550                        if(iomodel->domaintype!=Domain2DhorizontalEnum) _error_("not implemented yet");
     551                        CreateEdges(iomodel);
     552                        CreateFaces(iomodel);
     553                        for(int i=0;i<iomodel->numberoffaces;i++){
     554                                int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     555                                int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
     556                                if(e2!=-2){
     557                                        if(epart[e1]!=epart[e2]){
     558                                                AddNodeToRank(nodes_ranks,nodes_proc_count,e2,epart[e1]);
     559                                                //AddNodeToRank(nodes_ranks,nodes_proc_count,e1,epart[e2]);
     560                                        }
     561                                }
     562                        }
     563                }
     564                else{
     565                        _error_("not supported");
    545566                }
    546567        }/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.