Changeset 17429


Ignore:
Timestamp:
03/13/14 15:32:02 (11 years ago)
Author:
jbondzio
Message:

BUG: deactivate nodes on non-basal layer for 2d solving process of extrapolation and LSM

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

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp

    r17428 r17429  
    2929                }
    3030        }
    31        
    3231        if(iomodel->meshtype==Mesh3DEnum){
    3332                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     33                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
    3434        }
    3535}
    … …  
    8686
    8787        /*Intermediaries */
    88         int               meshtype,dim;
     88        int                meshtype,dim;
    8989        int        i,row,col,stabilization;
    9090        bool       extrapolatebydiffusion = true;
    … …  
    110110        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
    111111        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
    112         IssmDouble*             dlsf   = xNew<IssmDouble>(dim);
    113         IssmDouble*             normal = xNew<IssmDouble>(dim);
    114         IssmDouble*    D                 = xNew<IssmDouble>(dim*dim);
     112        IssmDouble*             D         = xNew<IssmDouble>(dim*dim);
     113        IssmDouble*             dlsf = xNew<IssmDouble>(dim);
     114        IssmDouble*             normal= xNew<IssmDouble>(dim);
    115115
    116116        /*Retrieve all inputs and parameters*/
    … …  
    162162                                for(col=0;col<dim;col++)
    163163                                        if(row==col)
    164                                                 D[row*dim+col] = D_scalar*normal[row];
     164                                                D[row*dim+col]=D_scalar*normal[row];
    165165                                        else
    166                                                 D[row*dim+col] = 0.;
     166                                                D[row*dim+col]=0.;
    167167                        TripleMultiply(B,dim,numnodes,1,
    168168                                                D,dim,dim,0,
    … …  
    181181                                        for(col=0;col<dim;col++)
    182182                                                if(row==col)
    183                                                         D[row*dim+col] = D_scalar*kappa;
     183                                                        D[row*dim+col]=D_scalar*kappa;
    184184                                                else
    185                                                         D[row*dim+col] = 0.;
     185                                                        D[row*dim+col]=0.;
    186186                                TripleMultiply(Bprime,dim,numnodes,1,
    187187                                                        D,dim,dim,0,
    … …  
    207207        xDelete<IssmDouble>(B);
    208208        xDelete<IssmDouble>(Bprime);
     209        xDelete<IssmDouble>(D);
    209210        xDelete<IssmDouble>(dlsf);
    210211        xDelete<IssmDouble>(normal);
     212        delete gauss;
    211213        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    212         delete gauss;
    213214        return Ke;
    214215
    … …  
    224225        /*Fetch number of nodes */
    225226        int numnodes = basalelement->GetNumberOfNodes();
    226         basalelement->FindParam(&meshtype,MeshTypeEnum);
    227227
    228228        /*Initialize Element vector*/
    … …  
    230230        for(i=0;i<numnodes;i++)
    231231                pe->values[i]=0.;
     232
     233        basalelement->FindParam(&meshtype,MeshTypeEnum);
    232234        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    233235        return pe;
    … …  
    326328                node=element->GetNode(in);
    327329                levelset_input->GetInputValue(&phi,gauss);
    328                 if(phi<=0.){
    329                         /* if ice, set dirichlet BC */
    330                         extvar_input->GetInputValue(&value,gauss);
    331                         node->ApplyConstraint(1,value);
    332                 }
    333                 else {
    334                         /* no ice, set no spc */
    335                         node->DofInFSet(0);
     330                if(node->IsActive()){
     331                        if(phi<=0.){
     332                                /* if ice, set dirichlet BC */
     333                                extvar_input->GetInputValue(&value,gauss);
     334                                node->ApplyConstraint(1,value);
     335                        }
     336                        else {
     337                                /* no ice, set no spc */
     338                                node->DofInFSet(0);
     339                        }
    336340                }
    337341        }

  • TabularUnified issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

    r17427 r17429  
    3737        if(iomodel->meshtype==Mesh3DEnum){
    3838                iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
     39                iomodel->FetchDataToInput(elements,MeshElementonsurfaceEnum);
    3940        }
    4041}
    … …  
    101102
    102103        /*Intermediaries */
    103         int  dim, meshtype; // solve for LSF in horizontal plane only
     104        int  dim, meshtype;
    104105        int i, row, col;
    105106        IssmDouble kappa;
    … …  
    136137        basalelement->GetVerticesCoordinates(&xyz_list);
    137138        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    138         Input* vx_input  = NULL;
    139         Input* vy_input  = NULL;
     139        Input* vx_input=NULL;
     140        Input* vy_input=NULL;
     141        if(meshtype==Mesh2DhorizontalEnum){
     142                vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     143                vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
     144        }
     145        else{
     146                if(dim==1){
     147                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     148                }
     149                if(dim==2){
     150                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
     151                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
     152                }
     153        }
     154
    140155        Input* lsf_slopex_input  = basalelement->GetInput(LevelsetfunctionSlopeXEnum); _assert_(lsf_slopex_input);
    141156        Input* lsf_slopey_input  = basalelement->GetInput(LevelsetfunctionSlopeYEnum); _assert_(lsf_slopey_input);
    142157        Input* calvingrate_input  = basalelement->GetInput(MasstransportCalvingrateEnum);     _assert_(calvingrate_input);
    143 
    144         if(meshtype==Mesh2DhorizontalEnum){
    145                 vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    146                 vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
    147         }
    148         else{
    149                 if(dim==1){
    150                         vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
    151                 }
    152                 if(dim==2){
    153                         vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
    154                         vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
    155                 }
    156         }
    157158       
    158159        /* Start  looping on the number of gaussian points: */
    … …  
    260261        xDelete<IssmDouble>(c);
    261262        xDelete<IssmDouble>(dlsf);
     263        delete gauss;
    262264        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    263         delete gauss;
    264265        return Ke;
    265266}/*}}}*/
    … …  
    306307                xDelete<IssmDouble>(xyz_list);
    307308                xDelete<IssmDouble>(basis);
     309                basalelement->FindParam(&meshtype,MeshTypeEnum);
     310                if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    308311                delete gauss;
    309312        }
    … …  
    311314                for(i=0;i<numnodes;i++)
    312315                        pe->values[i]=0.;
    313         if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    314316        return pe;
    315317}/*}}}*/

  • TabularUnified issm/trunk-jpl/src/c/classes/Node.cpp

    r17323 r17429  
    9797                                analysis_enum==HydrologyDCInefficientAnalysisEnum ||
    9898                                analysis_enum==DamageEvolutionAnalysisEnum ||
    99                                 analysis_enum==HydrologyDCEfficientAnalysisEnum
     99                                analysis_enum==HydrologyDCEfficientAnalysisEnum ||
     100                                analysis_enum==LevelsetAnalysisEnum ||
     101                                analysis_enum==ExtrapolationAnalysisEnum
    100102                                ){
    101103                if(iomodel->meshtype==Mesh3DEnum || iomodel->meshtype==Mesh2DverticalEnum){
Note: See TracChangeset for help on using the changeset viewer.