Ignore:
Timestamp:
04/01/20 21:54:40 (5 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 24684

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

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

    r24313 r24686  
    111111        return 1;
    112112}/*}}}*/
    113 void MasstransportAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
     113void MasstransportAnalysis::UpdateElements(Elements* elements,Inputs2* inputs2,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    114114
    115115        int    stabilization,finiteelement;
    … …  
    140140                if(iomodel->my_elements[i]){
    141141                        Element* element=(Element*)elements->GetObjectByOffset(counter);
    142                         element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
     142                        element->Update(inputs2,i,iomodel,analysis_counter,analysis_type,finiteelement);
    143143                        counter++;
    144144                }
    145145        }
    146146
    147         iomodel->FetchDataToInput(elements,"md.geometry.thickness",ThicknessEnum);
    148         iomodel->FetchDataToInput(elements,"md.geometry.surface",SurfaceEnum);
    149         iomodel->FetchDataToInput(elements,"md.geometry.base",BaseEnum);
    150         iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
    151         iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    152         iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
    153         iomodel->FetchDataToInput(elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
    154         iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
    155         iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
    156         if(isgroundingline)     iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum);
     147        iomodel->FetchDataToInput(inputs2,elements,"md.geometry.thickness",ThicknessEnum);
     148        iomodel->FetchDataToInput(inputs2,elements,"md.geometry.surface",SurfaceEnum);
     149        iomodel->FetchDataToInput(inputs2,elements,"md.geometry.base",BaseEnum);
     150        iomodel->FetchDataToInput(inputs2,elements,"md.slr.sealevel",SealevelEnum,0);
     151        iomodel->FetchDataToInput(inputs2,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
     152        iomodel->FetchDataToInput(inputs2,elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
     153        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum);
     154        iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vx",VxEnum);
     155        iomodel->FetchDataToInput(inputs2,elements,"md.initialization.vy",VyEnum);
     156        if(isgroundingline)     iomodel->FetchDataToInput(inputs2,elements,"md.geometry.bed",BedEnum);
    157157        /*Initialize cumdeltalthickness input*/
    158         InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum);
     158        InputUpdateFromConstantx(inputs2,elements,0.,SealevelriseCumDeltathicknessEnum);
    159159        /*Initialize ThicknessResidual input*/
    160         InputUpdateFromConstantx(elements,0.,ThicknessResidualEnum);
     160        InputUpdateFromConstantx(inputs2,elements,0.,ThicknessResidualEnum);
    161161
    162162        /*Get what we need for ocean-induced basal melting*/
    … …  
    165165        switch(basalforcing_model){
    166166                case FloatingMeltRateEnum:
    167                         iomodel->FetchDataToInput(elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
     167                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.floatingice_melting_rate",BasalforcingsFloatingiceMeltingRateEnum);
    168168                        break;
    169169                case LinearFloatingMeltRateEnum:
     170                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.perturbation_melting_rate",BasalforcingsPerturbationMeltingRateEnum,0.);
    170171                        break;
    171172                case MismipFloatingMeltRateEnum:
    … …  
    174175                        break;
    175176                case SpatialLinearFloatingMeltRateEnum:
    176                         iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
    177                         iomodel->FetchDataToInput(elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
    178                         iomodel->FetchDataToInput(elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
     177                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_melting_rate",BasalforcingsDeepwaterMeltingRateEnum);
     178                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.deepwater_elevation",BasalforcingsDeepwaterElevationEnum);
     179                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.upperwater_elevation",BasalforcingsUpperwaterElevationEnum);
    179180                        break;
    180181                case BasalforcingsPicoEnum:
    181                         iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
     182                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsPicoBasinIdEnum);
     183                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.overturning_coeff",BasalforcingsPicoOverturningCoeffEnum);
    182184                        break;
    183185                case BasalforcingsIsmip6Enum:{
    184                         iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
    185                         iomodel->FetchDataToInput(elements,"md.basalforcings.melt_anomaly",BasalforcingsIsmip6MeltAnomalyEnum,0.);
     186                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
     187                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.melt_anomaly",BasalforcingsIsmip6MeltAnomalyEnum,0.);
    186188                        IssmDouble** array3d = NULL; int* Ms = NULL; int* Ns = NULL; int K;
    187189                        iomodel->FetchData(&array3d,&Ms,&Ns,&K,"md.basalforcings.tf");
    … …  
    191193                                if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
    192194                                for(int kk=0;kk<K;kk++){
    193                                         element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array3d[kk],iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmip6TfEnum,7,kk);
     195                                        element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array3d[kk],inputs2,iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmip6TfEnum,7,kk);
    194196                                }
    195197                        }
    … …  
    200202                        break;
    201203                case BeckmannGoosseFloatingMeltRateEnum:
    202                         iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
    203                         iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
     204                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
     205                        iomodel->FetchDataToInput(inputs2,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
    204206                        break;
    205207                default:
    … …  
    208210
    209211        if(!issmb){
    210                 iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum);
     212                iomodel->FetchDataToInput(inputs2,elements,"md.smb.mass_balance",SmbMassBalanceEnum);
    211213        }
    212214        if(stabilization==3){
    213                 iomodel->FetchDataToInput(elements,"md.masstransport.spcthickness",MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
     215                iomodel->FetchDataToInput(inputs2,elements,"md.masstransport.spcthickness",MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
    214216        }
    215217        if(stabilization==4){
    216                 iomodel->FetchDataToInput(elements,"md.masstransport.spcthickness",MasstransportSpcthicknessEnum); //for FCT, we need the spc in the element (penlaties)
     218                iomodel->FetchDataToInput(inputs2,elements,"md.masstransport.spcthickness",MasstransportSpcthicknessEnum); //for FCT, we need the spc in the element (penlaties)
    217219        }
    218220
    219221        if(iomodel->domaintype!=Domain2DhorizontalEnum){
    220                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
    221                 iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
     222                iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
     223                iomodel->FetchDataToInput(inputs2,elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
    222224        }
    223225
    … …  
    291293        IssmDouble xi,tau;
    292294        IssmDouble dvx[2],dvy[2];
     295        IssmDouble D[4];
    293296        IssmDouble* xyz_list = NULL;
    294297
    … …  
    309312        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    310313        IssmDouble*             dbasis = xNew<IssmDouble>(dim*numnodes);
    311         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    312         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
    313         IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
    314314
    315315        /*Retrieve all inputs and parameters*/
    … …  
    318318        element->FindParam(&domaintype,DomainTypeEnum);
    319319        element->FindParam(&stabilization,MasstransportStabilizationEnum);
    320         Input* vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    321         Input* vyaverage_input=NULL;
     320        Input2* vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
     321        Input2* vyaverage_input=NULL;
    322322        if(dim==2){
    323                 vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     323                vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
    324324        }
    325325
    … …  
    333333                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    334334                element->NodalFunctions(basis,gauss);
    335 
     335                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     336
     337                /*Transient term*/
     338                D_scalar=gauss->weight*Jdet;
     339                for(int i=0;i<numnodes;i++){
     340                        for(int j=0;j<numnodes;j++){
     341                                Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
     342                        }
     343                }
     344
     345                /*Advection terms*/
    336346                vxaverage_input->GetInputValue(&vx,gauss);
    337347                vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     348                D_scalar=dt*gauss->weight*Jdet;
    338349                if(dim==2){
    339350                        vyaverage_input->GetInputValue(&vy,gauss);
    340351                        vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
    341                 }
    342 
    343                 D_scalar=gauss->weight*Jdet;
    344                 TripleMultiply(basis,1,numnodes,1,
    345                                         &D_scalar,1,1,0,
    346                                         basis,1,numnodes,0,
    347                                         &Ke->values[0],1);
    348 
    349                 GetB(B,element,dim,xyz_list,gauss);
    350                 GetBprime(Bprime,element,dim,xyz_list,gauss);
    351 
    352                 dvxdx=dvx[0];
    353                 if(dim==2) dvydy=dvy[1];
    354                 D_scalar=dt*gauss->weight*Jdet;
    355 
    356                 D[0*dim+0]=D_scalar*dvxdx;
    357                 if(dim==2) D[1*dim+1]=D_scalar*dvydy;
    358 
    359                 TripleMultiply(B,dim,numnodes,1,
    360                                         D,dim,dim,0,
    361                                         B,dim,numnodes,0,
    362                                         &Ke->values[0],1);
    363 
    364                 D[0*dim+0]=D_scalar*vx;
    365                 if(dim==2) D[1*dim+1]=D_scalar*vy;
    366 
    367                 TripleMultiply(B,dim,numnodes,1,
    368                                         D,dim,dim,0,
    369                                         Bprime,dim,numnodes,0,
    370                                         &Ke->values[0],1);
    371 
     352                        dvxdx=dvx[0];
     353                        dvydy=dvy[1];
     354                        for(int i=0;i<numnodes;i++){
     355                                for(int j=0;j<numnodes;j++){
     356                                        /*\phi_i \phi_j \nabla\cdot v*/
     357                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j]*(dvxdx+dvydy);
     358                                        /*\phi_i v\cdot\nabla\phi_j*/
     359                                        Ke->values[i*numnodes+j] += D_scalar*basis[i]*(vx*dbasis[0*numnodes+j] + vy*dbasis[1*numnodes+j]);
     360                                }
     361                        }
     362                }
     363                else{
     364                        dvxdx=dvx[0];
     365                        for(int i=0;i<numnodes;i++){
     366                                for(int j=0;j<numnodes;j++){
     367                                        Ke->values[i*numnodes+j] += D_scalar*dvxdx*dbasis[0*numnodes+i]*dbasis[0*numnodes+j];
     368                                        Ke->values[i*numnodes+j] += D_scalar*vx*dbasis[0*numnodes+j]*basis[i];
     369                                }
     370                        }
     371                }
     372
     373                for(int i=0;i<4;i++) D[i]=0.;
    372374                switch(stabilization){
    373375                        case 0:
    … …  
    383385                        case 2:
    384386                                /*Streamline upwinding*/
    385                                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    386387                                vxaverage_input->GetInputAverage(&vx);
    387388                                if(dim==1){
    … …  
    397398                                /*SUPG*/
    398399                                if(dim!=2) _error_("Stabilization "<<stabilization<<" not supported yet for dim != 2");
    399                                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    400400                                vxaverage_input->GetInputAverage(&vx);
    401401                                vyaverage_input->GetInputAverage(&vy);
    … …  
    418418                        }
    419419
    420                         TripleMultiply(Bprime,dim,numnodes,1,
    421                                                 D,dim,dim,0,
    422                                                 Bprime,dim,numnodes,0,
    423                                                 &Ke->values[0],1);
     420                        if(dim==2){
     421                                for(int i=0;i<numnodes;i++){
     422                                        for(int j=0;j<numnodes;j++){
     423                                                Ke->values[i*numnodes+j] += (
     424                                                                        dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
     425                                                                        dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
     426                                                                        );
     427                                        }
     428                                }
     429                        }
     430                        else{
     431                                for(int i=0;i<numnodes;i++){
     432                                        for(int j=0;j<numnodes;j++){
     433                                                Ke->values[i*numnodes+j] += D_scalar*h/(2.*vel)*dbasis[0*numnodes+i] *D[0]*dbasis[0*numnodes+j];
     434                                        }
     435                                }
     436                        }
    424437                }
    425438                if(stabilization==2){
    426439                        /*Streamline upwind*/
     440                        _assert_(dim==2);
    427441                        for(int i=0;i<numnodes;i++){
    428442                                for(int j=0;j<numnodes;j++){
    … …  
    477491        xDelete<IssmDouble>(basis);
    478492        xDelete<IssmDouble>(dbasis);
    479         xDelete<IssmDouble>(B);
    480         xDelete<IssmDouble>(Bprime);
    481         xDelete<IssmDouble>(D);
    482493        delete gauss;
    483494        return Ke;
    … …  
    499510        ElementMatrix* Ke     = element->NewElementMatrix();
    500511        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    501         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    502         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    503         IssmDouble     D[2][2];
     512        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
    504513
    505514        /*Retrieve all inputs and parameters*/
    … …  
    507516        element->FindParam(&dt,TimesteppingTimeStepEnum);
    508517        element->FindParam(&domaintype,DomainTypeEnum);
    509         Input* vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    510         Input* vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     518        Input2* vxaverage_input=element->GetInput2(VxAverageEnum); _assert_(vxaverage_input);
     519        Input2* vyaverage_input=element->GetInput2(VyAverageEnum); _assert_(vyaverage_input);
    511520
    512521        /* Start  looping on the number of gaussian points: */
    … …  
    517526                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    518527                element->NodalFunctions(basis,gauss);
     528                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    519529
    520530                vxaverage_input->GetInputValue(&vx,gauss);
    … …  
    522532
    523533                D_scalar=gauss->weight*Jdet;
    524                 TripleMultiply(basis,1,numnodes,1,
    525                                         &D_scalar,1,1,0,
    526                                         basis,1,numnodes,0,
    527                                         &Ke->values[0],1);
    528 
    529                 /*WARNING: B and Bprime are inverted compared to CG*/
    530                 GetB(Bprime,element,2,xyz_list,gauss);
    531                 GetBprime(B,element,2,xyz_list,gauss);
    532 
     534                for(int i=0;i<numnodes;i++){
     535                        for(int j=0;j<numnodes;j++){
     536                                Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
     537                        }
     538                }
     539
     540                /*WARNING: basis and dbasis are inverted compared to CG*/
    533541                D_scalar = - dt*gauss->weight*Jdet;
    534                 D[0][0]  = D_scalar*vx;
    535                 D[0][1]  = 0.;
    536                 D[1][0]  = 0.;
    537                 D[1][1]  = D_scalar*vy;
    538                 TripleMultiply(B,2,numnodes,1,
    539                                         &D[0][0],2,2,0,
    540                                         Bprime,2,numnodes,0,
    541                                         &Ke->values[0],1);
    542 
     542                for(int i=0;i<numnodes;i++){
     543                        for(int j=0;j<numnodes;j++){
     544                                Ke->values[i*numnodes+j] += D_scalar*(vx*dbasis[0*numnodes+i]*basis[j] + vy*dbasis[1*numnodes+i]*basis[j]);
     545                        }
     546                }
    543547        }
    544548
    … …  
    546550        xDelete<IssmDouble>(xyz_list);
    547551        xDelete<IssmDouble>(basis);
    548         xDelete<IssmDouble>(B);
    549         xDelete<IssmDouble>(Bprime);
     552        xDelete<IssmDouble>(dbasis);
    550553        delete gauss;
    551554        return Ke;
    … …  
    617620        element->FindParam(&dt,TimesteppingTimeStepEnum);
    618621        element->FindParam(&stabilization,MasstransportStabilizationEnum);
    619         Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
    620         Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
    621         Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
    622         Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
    623         Input* thickness_input  = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
    624         Input* vxaverage_input  = element->GetInput(VxAverageEnum);                                                                             _assert_(vxaverage_input);
    625         Input* vyaverage_input  = element->GetInput(VyAverageEnum);                                                                             _assert_(vyaverage_input);
     622        Input2* gmb_input        = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
     623        Input2* fmb_input        = element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
     624        Input2* gllevelset_input = element->GetInput2(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
     625        Input2* ms_input         = element->GetInput2(SmbMassBalanceEnum);                       _assert_(ms_input);
     626        Input2* thickness_input  = element->GetInput2(ThicknessEnum);                            _assert_(thickness_input);
     627        Input2* vxaverage_input  = element->GetInput2(VxAverageEnum);                                                                           _assert_(vxaverage_input);
     628        Input2* vyaverage_input  = element->GetInput2(VyAverageEnum);                                                                           _assert_(vyaverage_input);
     629
     630//      if(element->Id()==9){
     631//              gmb_input->Echo();
     632//              _error_("S");
     633//      }
    626634
    627635        h=element->CharacteristicLength();
    … …  
    726734        element->FindParam(&dt,TimesteppingTimeStepEnum);
    727735        element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
    728         Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
    729         Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
    730         Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                      _assert_(ms_input);
    731         Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
    732         Input* thickness_input  = element->GetInput(ThicknessEnum);                           _assert_(thickness_input);
     736        Input2* gmb_input        = element->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(gmb_input);
     737        Input2* fmb_input        = element->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(fmb_input);
     738        Input2* ms_input         = element->GetInput2(SmbMassBalanceEnum);                      _assert_(ms_input);
     739        Input2* gllevelset_input = element->GetInput2(MaskGroundediceLevelsetEnum);             _assert_(gllevelset_input);
     740        Input2* thickness_input  = element->GetInput2(ThicknessEnum);                           _assert_(thickness_input);
    733741
    734742   /*Recover portion of element that is grounded*/
    … …  
    783791        return pe;
    784792}/*}}}*/
    785 void           MasstransportAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    786         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    787          * For node i, Bi can be expressed in the actual coordinate system
    788          * by:
    789          *       Bi=[ N ]
    790          *          [ N ]
    791          * where N is the finiteelement function for node i.
    792          *
    793          * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
    794          */
    795 
    796         /*Fetch number of nodes for this finite element*/
    797         int numnodes = element->GetNumberOfNodes();
    798 
    799         /*Get nodal functions*/
    800         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    801         element->NodalFunctions(basis,gauss);
    802 
    803         /*Build B: */
    804         for(int i=0;i<numnodes;i++){
    805                 for(int j=0;j<dim;j++){
    806                         B[numnodes*j+i] = basis[i];
    807                 }
    808         }
    809 
    810         /*Clean-up*/
    811         xDelete<IssmDouble>(basis);
    812 }/*}}}*/
    813 void           MasstransportAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    814         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    815          * For node i, Bi' can be expressed in the actual coordinate system
    816          * by:
    817          *       Bi_prime=[ dN/dx ]
    818          *                [ dN/dy ]
    819          * where N is the finiteelement function for node i.
    820          *
    821          * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
    822          */
    823 
    824         /*Fetch number of nodes for this finite element*/
    825         int numnodes = element->GetNumberOfNodes();
    826 
    827         /*Get nodal functions derivatives*/
    828         IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    829         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    830 
    831         /*Build B': */
    832         for(int i=0;i<numnodes;i++){
    833                 for(int j=0;j<dim;j++){
    834                         Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
    835                 }
    836         }
    837 
    838         /*Clean-up*/
    839         xDelete<IssmDouble>(dbasis);
    840 
    841 }/*}}}*/
    842793void           MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    843794        element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
    … …  
    872823                }
    873824        }
    874         element->AddBasalInput(ThicknessEnum,newthickness,element->GetElementType());
    875         element->AddBasalInput(ThicknessResidualEnum,thicknessresidual,element->GetElementType());
     825        element->AddBasalInput2(ThicknessEnum,newthickness,element->GetElementType());
     826        element->AddBasalInput2(ThicknessResidualEnum,thicknessresidual,element->GetElementType());
    876827       
    877828        xDelete<int>(doflist);
    … …  
    905856        basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum);
    906857        basalelement->GetInputListOnVertices(&oldthickness[0],ThicknessOldEnum);
    907         basalelement->GetInputListOnVertices(&oldbase[0],BaseEnum);
    908         basalelement->GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
     858        basalelement->GetInputListOnVertices(&oldbase[0],BaseOldEnum);
     859        basalelement->GetInputListOnVertices(&oldsurface[0],SurfaceOldEnum);
    909860        basalelement->GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
    910861        basalelement->GetInputListOnVertices(&sealevel[0],SealevelEnum);
    911         basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
     862        basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessOldEnum);
    912863
    913864        /*Do we do grounding line migration?*/
    … …  
    953904
    954905        /*Add input to the element: */
    955         element->AddBasalInput(SurfaceEnum,newsurface,P1Enum);
    956         element->AddBasalInput(BaseEnum,newbase,P1Enum);
    957         element->AddBasalInput(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum);
    958         element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
     906        element->AddBasalInput2(SurfaceEnum,newsurface,P1Enum);
     907        element->AddBasalInput2(BaseEnum,newbase,P1Enum);
     908        element->AddBasalInput2(SealevelriseCumDeltathicknessEnum,cumdeltathickness,P1Enum);
     909        element->AddBasalInput2(SealevelriseDeltathicknessEnum,deltathickness,P1Enum);
    959910
    960911        /*Free ressources:*/
    … …  
    984935
    985936        /*Intermediaries */
    986         IssmDouble Jdet;
     937        IssmDouble Jdet,D_scalar;
    987938        IssmDouble vx,vy;
    988939        IssmDouble* xyz_list = NULL;
    … …  
    994945        /*Initialize Element vector and other vectors*/
    995946        ElementMatrix* Ke     = element->NewElementMatrix();
    996         IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    997         IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     947        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     948        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    998949        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
    999950
    1000951        /*Retrieve all inputs and parameters*/
    1001952        element->GetVerticesCoordinates(&xyz_list);
    1002         Input* vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
    1003         Input* vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input);
     953        Input2* vxaverage_input=element->GetInput2(VxEnum); _assert_(vxaverage_input);
     954        Input2* vyaverage_input=element->GetInput2(VyEnum); _assert_(vyaverage_input);
    1004955
    1005956        /* Start  looping on the number of gaussian points: */
    … …  
    1009960
    1010961                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    1011                 GetB(B,element,dim,xyz_list,gauss);
    1012                 GetBprime(Bprime,element,dim,xyz_list,gauss);
     962                element->NodalFunctions(basis,gauss);
     963                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    1013964                vxaverage_input->GetInputValue(&vx,gauss);
    1014965                vyaverage_input->GetInputValue(&vy,gauss);
    1015966
    1016                 D[0*dim+0] = -gauss->weight*vx*Jdet;
    1017                 D[1*dim+1] = -gauss->weight*vy*Jdet;
    1018 
    1019                 TripleMultiply(B,dim,numnodes,1,
    1020                                         D,dim,dim,0,
    1021                                         Bprime,dim,numnodes,0,
    1022                                         &Ke->values[0],1);
    1023 
     967                D_scalar = gauss->weight*Jdet;
     968                for(int i=0;i<numnodes;i++){
     969                        for(int j=0;j<numnodes;j++){
     970                                Ke->values[i*numnodes+j] += -D_scalar*(vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i]);
     971                        }
     972                }
    1024973        }
    1025974
    1026975        /*Clean up and return*/
    1027976        xDelete<IssmDouble>(xyz_list);
    1028         xDelete<IssmDouble>(B);
    1029         xDelete<IssmDouble>(Bprime);
    1030         xDelete<IssmDouble>(D);
     977        xDelete<IssmDouble>(basis);
     978        xDelete<IssmDouble>(dbasis);
    1031979        delete gauss;
    1032980        return Ke;
Note: See TracChangeset for help on using the changeset viewer.