Changeset 15695


Ignore:
Timestamp:
08/02/13 15:37:37 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: extending P1+ and P1+condensed to 3d

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r15694 r15695  
    391391}
    392392/*}}}*/
    393 /*FUNCTION Penta::CreateKMatrix {{{*/
     393/*FUNCTION Penta::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs) {{{*/
    394394void  Penta::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
    395 
    396         /*retrieve parameters: */
    397         ElementMatrix* Ke=NULL;
    398         int analysis_type;
    399         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    400 
    401         /*Checks in debugging {{{*/
    402         _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    403         /*}}}*/
    404395
    405396        /*Skip if water element*/
    406397        if(IsOnWater()) return;
    407398
    408         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    409         switch(analysis_type){
    410                 #ifdef _HAVE_DIAGNOSTIC_
    411                 case DiagnosticHorizAnalysisEnum:
    412                         Ke=CreateKMatrixDiagnosticHoriz();
    413                         break;
    414                 case AdjointHorizAnalysisEnum:
    415                         Ke=CreateKMatrixAdjointHoriz();
    416                         break;
    417                 case DiagnosticSIAAnalysisEnum:
    418                         Ke=CreateKMatrixDiagnosticSIA();
    419                         break;
    420                 case DiagnosticVertAnalysisEnum:
    421                         Ke=CreateKMatrixDiagnosticVert();
    422                         break;
    423                 #endif
    424                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    425                         Ke=CreateBasalMassMatrix();
    426                         break;
    427                 case PrognosticAnalysisEnum:
    428                         Ke=CreateKMatrixPrognostic();
    429                         break;
    430                 #ifdef _HAVE_BALANCED_
    431                 case BalancethicknessAnalysisEnum:
    432                         Ke=CreateKMatrixBalancethickness();
    433                         break;
    434                 #endif
    435                 #ifdef _HAVE_THERMAL_
    436                 case ThermalAnalysisEnum:
    437                         Ke=CreateKMatrixThermal();
    438                         break;
    439                 case EnthalpyAnalysisEnum:
    440                         Ke=CreateKMatrixEnthalpy();
    441                         break;
    442                 case MeltingAnalysisEnum:
    443                         Ke=CreateKMatrixMelting();
    444                         break;
    445                 #endif
    446                 #ifdef _HAVE_HYDROLOGY_
    447           case HydrologyDCInefficientAnalysisEnum:
    448                         Ke=CreateKMatrixHydrologyDCInefficient();
    449                         break;
    450           case HydrologyDCEfficientAnalysisEnum:
    451                         Ke=CreateKMatrixHydrologyDCEfficient();
    452                         break;
    453                 #endif
    454                 default:
    455                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    456         }
     399        /*Create element stiffness matrix*/
     400        ElementMatrix* Ke=CreateKMatrix();
    457401
    458402        if(Ke){
    … …  
    462406                        Ke->StaticCondensation(3,&indices[0]);
    463407                }
     408                else if(this->element_type==P1bubblecondensedEnum){
     409                        int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     410                        int offset = 0;
     411                        for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     412                        int* indices=xNew<int>(size);
     413                        for(int i=0;i<size;i++) indices[i] = offset+i;
     414                        Ke->StaticCondensation(size,indices);
     415                        xDelete<int>(indices);
     416                }
    464417
    465418                /*Add to global matrix*/
    … …  
    467420                delete Ke;
    468421        }
     422}
     423/*}}}*/
     424/*FUNCTION Penta::CreateKMatrix(void){{{*/
     425ElementMatrix* Penta::CreateKMatrix(void){
     426
     427        /*retrieve parameters: */
     428        ElementMatrix* Ke=NULL;
     429        int analysis_type;
     430        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     431
     432        /*Checks in debugging {{{*/
     433        _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
     434        /*}}}*/
     435
     436        /*Skip if water element*/
     437        if(IsOnWater()) return NULL;
     438
     439        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     440        switch(analysis_type){
     441                #ifdef _HAVE_DIAGNOSTIC_
     442                case DiagnosticHorizAnalysisEnum:
     443                        return CreateKMatrixDiagnosticHoriz();
     444                        break;
     445                case AdjointHorizAnalysisEnum:
     446                        return CreateKMatrixAdjointHoriz();
     447                        break;
     448                case DiagnosticSIAAnalysisEnum:
     449                        return CreateKMatrixDiagnosticSIA();
     450                        break;
     451                case DiagnosticVertAnalysisEnum:
     452                        return CreateKMatrixDiagnosticVert();
     453                        break;
     454                #endif
     455                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
     456                        return CreateBasalMassMatrix();
     457                        break;
     458                case PrognosticAnalysisEnum:
     459                        return CreateKMatrixPrognostic();
     460                        break;
     461                #ifdef _HAVE_BALANCED_
     462                case BalancethicknessAnalysisEnum:
     463                        return CreateKMatrixBalancethickness();
     464                        break;
     465                #endif
     466                #ifdef _HAVE_THERMAL_
     467                case ThermalAnalysisEnum:
     468                        return CreateKMatrixThermal();
     469                        break;
     470                case EnthalpyAnalysisEnum:
     471                        return CreateKMatrixEnthalpy();
     472                        break;
     473                case MeltingAnalysisEnum:
     474                        return CreateKMatrixMelting();
     475                        break;
     476                #endif
     477                #ifdef _HAVE_HYDROLOGY_
     478          case HydrologyDCInefficientAnalysisEnum:
     479                        return CreateKMatrixHydrologyDCInefficient();
     480                        break;
     481          case HydrologyDCEfficientAnalysisEnum:
     482                        return CreateKMatrixHydrologyDCEfficient();
     483                        break;
     484                #endif
     485                default:
     486                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
     487        }
     488
     489        /*Make compiler happy*/
     490        return NULL;
    469491}
    470492/*}}}*/
    … …  
    531553}
    532554/*}}}*/
    533 /*FUNCTION Penta::CreatePVector {{{*/
     555/*FUNCTION Penta::CreatePVector(Vector<IssmDouble>* pf) {{{*/
    534556void  Penta::CreatePVector(Vector<IssmDouble>* pf){
    535557
    536         /*retrive parameters: */
    537         ElementVector* pe=NULL;
     558        /*Skip if water element*/
     559        if(IsOnWater()) return;
     560
     561        /*Create element load vector*/
     562        ElementVector* pe = CreatePVector();
     563
     564        if(pe){
     565                /*StaticCondensation if requested*/
     566                if(this->element_type==MINIcondensedEnum){
     567                        int indices[3]={18,19,20};
     568
     569                        this->element_type=MINIEnum;
     570                        ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
     571                        this->element_type=MINIcondensedEnum;
     572
     573                        pe->StaticCondensation(Ke,3,&indices[0]);
     574                        delete Ke;
     575                }
     576                else if(this->element_type==P1bubblecondensedEnum){
     577                        int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     578                        int offset = 0;
     579                        for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     580                        int* indices=xNew<int>(size);
     581                        for(int i=0;i<size;i++) indices[i] = offset+i;
     582
     583                        this->element_type=P1bubbleEnum;
     584                        ElementMatrix* Ke = CreateKMatrix();
     585                        this->element_type=P1bubblecondensedEnum;
     586                        pe->StaticCondensation(Ke,size,indices);
     587                        xDelete<int>(indices);
     588                        delete Ke;
     589                }
     590
     591                /*Add to global Vector*/
     592                pe->AddToGlobal(pf);
     593                delete pe;
     594        }
     595}
     596/*}}}*/
     597/*FUNCTION Penta::CreatePVector(void) {{{*/
     598ElementVector* Penta::CreatePVector(void){
     599
     600        /*retrieve parameters: */
    538601        int analysis_type;
    539602        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    544607
    545608        /*Skip if water element*/
    546         if(IsOnWater()) return;
     609        if(IsOnWater()) return NULL;
    547610
    548611        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    … …  
    550613                #ifdef _HAVE_DIAGNOSTIC_
    551614                case DiagnosticHorizAnalysisEnum:
    552                         pe=CreatePVectorDiagnosticHoriz();
     615                        return CreatePVectorDiagnosticHoriz();
    553616                        break;
    554617                case DiagnosticSIAAnalysisEnum:
    555                         pe=CreatePVectorDiagnosticSIA();
     618                        return CreatePVectorDiagnosticSIA();
    556619                        break;
    557620                case DiagnosticVertAnalysisEnum:
    558                         pe=CreatePVectorDiagnosticVert();
     621                        return CreatePVectorDiagnosticVert();
    559622                        break;
    560623                #endif
    561624                #ifdef _HAVE_CONTROL_
    562625                case AdjointHorizAnalysisEnum:
    563                         pe=CreatePVectorAdjointHoriz();
     626                        return CreatePVectorAdjointHoriz();
    564627                        break;
    565628                #endif
    566629                #ifdef _HAVE_THERMAL_
    567630                case ThermalAnalysisEnum:
    568                         pe=CreatePVectorThermal();
     631                        return CreatePVectorThermal();
    569632                        break;
    570633                case EnthalpyAnalysisEnum:
    571                         pe=CreatePVectorEnthalpy();
     634                        return CreatePVectorEnthalpy();
    572635                        break;
    573636                case MeltingAnalysisEnum:
    574                         pe=CreatePVectorMelting();
     637                        return CreatePVectorMelting();
    575638                        break;
    576639                #endif
    577640                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    578                         pe=CreatePVectorSlope();
     641                        return CreatePVectorSlope();
    579642                        break;
    580643                case PrognosticAnalysisEnum:
    581                         pe=CreatePVectorPrognostic();
     644                        return CreatePVectorPrognostic();
    582645                        break;
    583646                #ifdef _HAVE_BALANCED_
    584647                case BalancethicknessAnalysisEnum:
    585                         pe=CreatePVectorBalancethickness();
     648                        return CreatePVectorBalancethickness();
    586649                        break;
    587650                #endif
    588651                #ifdef _HAVE_HYDROLOGY_
    589652          case HydrologyDCInefficientAnalysisEnum:
    590                         pe=CreatePVectorHydrologyDCInefficient();
     653                        return CreatePVectorHydrologyDCInefficient();
    591654                        break;
    592655          case HydrologyDCEfficientAnalysisEnum:
    593                         pe=CreatePVectorHydrologyDCEfficient();
     656                        return CreatePVectorHydrologyDCEfficient();
    594657                        break;
    595658                #endif
    … …  
    598661        }
    599662
    600         if(pe){
    601 
    602                 /*StaticCondensation if requested*/
    603                 if(this->element_type==MINIcondensedEnum){
    604                         int indices[3]={18,19,20};
    605 
    606                         this->element_type=MINIEnum;
    607                         ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
    608                         this->element_type=MINIcondensedEnum;
    609 
    610                         pe->StaticCondensation(Ke,3,&indices[0]);
    611                         delete Ke;
    612                 }
    613 
    614                 /*Add to global Vector*/
    615                 pe->AddToGlobal(pf);
    616                 delete pe;
    617         }
     663        /*make compiler happy*/
     664        return NULL;
    618665}
    619666/*}}}*/
    … …  
    83288375ElementVector* Penta::CreatePVectorDiagnosticHODrivingStress(void){
    83298376
    8330         /*Constants*/
    8331         const int    numdof=NDOF2*NUMVERTICES;
    8332 
    83338377        /*Intermediaries*/
    83348378        int         i,j;
    … …  
    83378381        IssmDouble  driving_stress_baseline,thickness;
    83388382        IssmDouble  xyz_list[NUMVERTICES][3];
    8339         IssmDouble  basis[6];
    83408383        GaussPenta  *gauss=NULL;
    83418384
    8342         /*Initialize Element vector*/
    8343         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
     8385        /*Fetch number of nodes and dof for this finite element*/
     8386        int numnodes = this->NumberofNodes();
     8387        int numdof   = numnodes*NDOF2;
     8388
     8389        /*Initialize vector*/
     8390        ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
     8391        IssmDouble*    basis = xNewZeroInit<IssmDouble>(numnodes);
    83448392
    83458393        /*Retrieve all inputs and parameters*/
    … …  
    83558403
    83568404                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8357                 GetNodalFunctionsP1(basis, gauss);
     8405                GetNodalFunctions(basis,gauss);
    83588406
    83598407                thickness_input->GetInputValue(&thickness, gauss);
    … …  
    83628410                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
    83638411
    8364                 for(i=0;i<NUMVERTICES;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
     8412                for(i=0;i<numnodes;i++){
     8413                        pe->values[i*NDOF2+0]+= -driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
     8414                        pe->values[i*NDOF2+1]+= -driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
     8415                }
    83658416        }
    83668417
    83678418        /*Transform coordinate system*/
    8368         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);
     8419        TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
    83698420
    83708421        /*Clean up and return*/
    83718422        delete gauss;
     8423        xDelete<IssmDouble>(basis);
    83728424        return pe;
    83738425}

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r15654 r15695  
    179179                void      BedNormal(IssmDouble* bed_normal, IssmDouble xyz_list[3][3]);
    180180                ElementMatrix* CreateBasalMassMatrix(void);
     181                ElementMatrix* CreateKMatrix(void);
    181182                ElementMatrix* CreateKMatrixPrognostic(void);
     183                ElementVector* CreatePVector(void);
    182184                ElementVector* CreatePVectorPrognostic(void);
    183185                ElementVector* CreatePVectorSlope(void);

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

    r15694 r15695  
    172172}
    173173/*}}}*/
    174 /*FUNCTION Tria::CreateKMatrix {{{*/
     174/*FUNCTION Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs) {{{*/
    175175void  Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
    176 
    177         /*retreive parameters: */
    178         ElementMatrix* Ke=NULL;
    179         int analysis_type;
    180         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    181 
    182         /*Checks in debugging mode{{{*/
    183         _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);
    184         /*}}}*/
    185176
    186177        /*Skip if water element*/
    187178        if(IsOnWater()) return;
    188179
    189         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    190         switch(analysis_type){
    191                 #ifdef _HAVE_DIAGNOSTIC_
    192                 case DiagnosticHorizAnalysisEnum:
    193                         Ke=CreateKMatrixDiagnosticSSA();
    194                         break;
    195                 case DiagnosticSIAAnalysisEnum:
    196                         Ke=CreateKMatrixDiagnosticSIA();
    197                         break;
    198                  #endif
    199                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    200                         Ke=CreateMassMatrix();
    201                         break;
    202                 case PrognosticAnalysisEnum:
    203                         Ke=CreateKMatrixPrognostic();
    204                         break;
    205                 #ifdef _HAVE_HYDROLOGY_
    206                 case HydrologyShreveAnalysisEnum:
    207                         Ke=CreateKMatrixHydrologyShreve();
    208                         break;
    209                 case HydrologyDCInefficientAnalysisEnum:
    210                         Ke=CreateKMatrixHydrologyDCInefficient();
    211                         break;
    212                 case HydrologyDCEfficientAnalysisEnum:
    213                         Ke=CreateKMatrixHydrologyDCEfficient();
    214                         break;
    215                 #endif
    216                 #ifdef _HAVE_BALANCED_
    217                 case BalancethicknessAnalysisEnum:
    218                         Ke=CreateKMatrixBalancethickness();
    219                         break;
    220                 #endif
    221                 #ifdef _HAVE_CONTROL_
    222                 case AdjointBalancethicknessAnalysisEnum:
    223                         Ke=CreateKMatrixAdjointBalancethickness();
    224                         break;
    225                 case AdjointHorizAnalysisEnum:
    226                         Ke=CreateKMatrixAdjointSSA();
    227                         break;
    228                 #endif
    229                 default:
    230                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    231         }
     180        /*Create element stiffness matrix*/
     181        ElementMatrix* Ke=CreateKMatrix();
    232182
    233183        if(Ke){
    … …  
    239189                        int* indices=xNew<int>(size);
    240190                        for(int i=0;i<size;i++) indices[i] = offset+i;
    241                         printarray(indices,1,size);
    242191                        Ke->StaticCondensation(size,indices);
    243192                        xDelete<int>(indices);
    … …  
    248197                delete Ke;
    249198        }
     199}
     200/*}}}*/
     201/*FUNCTION Tria::CreateKMatrix(void) {{{*/
     202ElementMatrix* Tria::CreateKMatrix(void){
     203
     204        /*retreive parameters: */
     205        int analysis_type;
     206        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     207
     208        /*Checks in debugging mode{{{*/
     209        _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);
     210        /*}}}*/
     211
     212        /*Skip if water element*/
     213        if(IsOnWater()) return NULL;
     214
     215        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     216        switch(analysis_type){
     217                #ifdef _HAVE_DIAGNOSTIC_
     218                case DiagnosticHorizAnalysisEnum:
     219                        return CreateKMatrixDiagnosticSSA();
     220                        break;
     221                case DiagnosticSIAAnalysisEnum:
     222                        return CreateKMatrixDiagnosticSIA();
     223                        break;
     224                 #endif
     225                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
     226                        return CreateMassMatrix();
     227                        break;
     228                case PrognosticAnalysisEnum:
     229                        return CreateKMatrixPrognostic();
     230                        break;
     231                #ifdef _HAVE_HYDROLOGY_
     232                case HydrologyShreveAnalysisEnum:
     233                        return CreateKMatrixHydrologyShreve();
     234                        break;
     235                case HydrologyDCInefficientAnalysisEnum:
     236                        return CreateKMatrixHydrologyDCInefficient();
     237                        break;
     238                case HydrologyDCEfficientAnalysisEnum:
     239                        return CreateKMatrixHydrologyDCEfficient();
     240                        break;
     241                #endif
     242                #ifdef _HAVE_BALANCED_
     243                case BalancethicknessAnalysisEnum:
     244                        return CreateKMatrixBalancethickness();
     245                        break;
     246                #endif
     247                #ifdef _HAVE_CONTROL_
     248                case AdjointBalancethicknessAnalysisEnum:
     249                        return CreateKMatrixAdjointBalancethickness();
     250                        break;
     251                case AdjointHorizAnalysisEnum:
     252                        return CreateKMatrixAdjointSSA();
     253                        break;
     254                #endif
     255                default:
     256                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
     257        }
     258
     259        /*Make compiler happy*/
     260        return NULL;
    250261}
    251262/*}}}*/
    … …  
    296307}
    297308/*}}}*/
    298 /*FUNCTION Tria::CreatePVector {{{*/
     309/*FUNCTION Tria::CreatePVector(Vector<IssmDouble>* pf) {{{*/
    299310void  Tria::CreatePVector(Vector<IssmDouble>* pf){
    300 
    301         /*retrive parameters: */
    302         ElementVector* pe=NULL;
    303         int analysis_type;
    304         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    305 
    306         /*asserts: {{{*/
    307         /*if debugging mode, check that all pointers exist*/
    308         _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);
    309         /*}}}*/
    310311
    311312        /*Skip if water element*/
    312313        if(IsOnWater()) return;
    313314
    314         /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    315         switch(analysis_type){
    316                 #ifdef _HAVE_DIAGNOSTIC_
    317                 case DiagnosticHorizAnalysisEnum:
    318                         pe=CreatePVectorDiagnosticSSA();
    319                         break;
    320                 case DiagnosticSIAAnalysisEnum:
    321                         pe=CreatePVectorDiagnosticSIA();
    322                         break;
    323                 #endif
    324                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
    325                         pe=CreatePVectorSlope();
    326                         break;
    327                 case PrognosticAnalysisEnum:
    328                         pe=CreatePVectorPrognostic();
    329                         break;
    330                 #ifdef _HAVE_HYDROLOGY_
    331                 case HydrologyShreveAnalysisEnum:
    332                         pe=CreatePVectorHydrologyShreve();
    333                         break;
    334                 case HydrologyDCInefficientAnalysisEnum:
    335                         pe=CreatePVectorHydrologyDCInefficient();
    336                         break;
    337                 case HydrologyDCEfficientAnalysisEnum:
    338                         pe=CreatePVectorHydrologyDCEfficient();
    339                         break;
    340                 #endif
    341                 #ifdef _HAVE_BALANCED_
    342                 case BalancethicknessAnalysisEnum:
    343                         pe=CreatePVectorBalancethickness();
    344                         break;
    345                 #endif
    346                 #ifdef _HAVE_CONTROL_
    347                 case AdjointBalancethicknessAnalysisEnum:
    348                         pe=CreatePVectorAdjointBalancethickness();
    349                         break;
    350                 case AdjointHorizAnalysisEnum:
    351                         pe=CreatePVectorAdjointHoriz();
    352                         break;
    353                 #endif
    354                 default:
    355                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    356         }
     315        /*Create element load vector*/
     316        ElementVector* pe = CreatePVector();
    357317
    358318        /*Add to global Vector*/
    … …  
    367327
    368328                        this->element_type=P1bubbleEnum;
    369                         ElementMatrix* Ke = CreateKMatrixDiagnosticSSA();
     329                        ElementMatrix* Ke = CreateKMatrix();
    370330                        this->element_type=P1bubblecondensedEnum;
    371331
    … …  
    378338                delete pe;
    379339        }
     340}
     341/*}}}*/
     342/*FUNCTION Tria::CreatePVector(void){{{*/
     343ElementVector* Tria::CreatePVector(void){
     344
     345        /*retrive parameters: */
     346        ElementVector* pe=NULL;
     347        int analysis_type;
     348        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     349
     350        /*asserts: {{{*/
     351        /*if debugging mode, check that all pointers exist*/
     352        _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);
     353        /*}}}*/
     354
     355        /*Skip if water element*/
     356        if(IsOnWater()) return NULL;
     357
     358        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
     359        switch(analysis_type){
     360#ifdef _HAVE_DIAGNOSTIC_
     361                case DiagnosticHorizAnalysisEnum:
     362                        return CreatePVectorDiagnosticSSA();
     363                        break;
     364                case DiagnosticSIAAnalysisEnum:
     365                        return CreatePVectorDiagnosticSIA();
     366                        break;
     367#endif
     368                case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
     369                        return CreatePVectorSlope();
     370                        break;
     371                case PrognosticAnalysisEnum:
     372                        return CreatePVectorPrognostic();
     373                        break;
     374#ifdef _HAVE_HYDROLOGY_
     375                case HydrologyShreveAnalysisEnum:
     376                        return CreatePVectorHydrologyShreve();
     377                        break;
     378                case HydrologyDCInefficientAnalysisEnum:
     379                        return CreatePVectorHydrologyDCInefficient();
     380                        break;
     381                case HydrologyDCEfficientAnalysisEnum:
     382                        return CreatePVectorHydrologyDCEfficient();
     383                        break;
     384#endif
     385#ifdef _HAVE_BALANCED_
     386                case BalancethicknessAnalysisEnum:
     387                        return CreatePVectorBalancethickness();
     388                        break;
     389#endif
     390#ifdef _HAVE_CONTROL_
     391                case AdjointBalancethicknessAnalysisEnum:
     392                        return CreatePVectorAdjointBalancethickness();
     393                        break;
     394                case AdjointHorizAnalysisEnum:
     395                        return CreatePVectorAdjointHoriz();
     396                        break;
     397#endif
     398                default:
     399                        _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
     400        }
     401
     402        /*make compiler happy*/
     403        return NULL;
    380404}
    381405/*}}}*/
    … …  
    31533177
    31543178                /*Build load vector: */
    3155                 for (i=0;i<numnodes;i++){
    3156                         for(j=0;j<NDOF2;j++){
    3157                                 pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
    3158                         }
     3179                for(i=0;i<numnodes;i++){
     3180                        pe->values[i*NDOF2+0]+=-driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
     3181                        pe->values[i*NDOF2+1]+=-driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
    31593182                }
    31603183        }

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r15613 r15695  
    179179                /*}}}*/
    180180                /*Tria specific routines:{{{*/
     181                ElementMatrix* CreateKMatrix(void);
    181182                ElementMatrix* CreateKMatrixBalancethickness(void);
    182183                ElementMatrix* CreateKMatrixBalancethickness_DG(void);
    … …  
    187188                ElementMatrix* CreateKMatrixPrognostic_DG(void);
    188189                ElementMatrix* CreateMassMatrix(void);
     190                ElementVector* CreatePVector(void);
    189191                ElementVector* CreatePVectorBalancethickness(void);
    190192                ElementVector* CreatePVectorBalancethickness_DG(void);
Note: See TracChangeset for help on using the changeset viewer.