Changeset 10440


Ignore:
Timestamp:
11/02/11 16:51:01 (13 years ago)
Author:
Mathieu Morlighem
Message:

Temporary commit for spc in Stokes under grounded ice, still not perfect

Location:
issm/trunk/src
Files:
7 added
1 deleted
16 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/DataSet.cpp

    r10300 r10440  
    503503
    504504        /*Check index in debugging mode*/
     505        _assert_(this!=NULL);
    505506        _assert_(offset<this->Size());
    506507
    … …  
    572573/*FUNCTION DataSet::Size{{{1*/
    573574int  DataSet::Size(void){
     575        _assert_(this!=NULL);
    574576
    575577        return objects.size();

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h

    r10287 r10440  
    484484        SpcDynamicEnum,
    485485        IceVolumeEnum,
     486        DiagnosticStokesBasalPenaltyEnum,
    486487        MaximumNumberOfEnums
    487488};

  • issm/trunk/src/c/Makefile.am

    r10407 r10440  
    330330                                        ./modules/ResetConstraintsx/ResetConstraintsx.h\
    331331                                        ./modules/ResetConstraintsx/ResetConstraintsx.cpp\
     332                                        ./modules/ResetCoordinateSystemx/ResetCoordinateSystemx.h\
     333                                        ./modules/ResetCoordinateSystemx/ResetCoordinateSystemx.cpp\
    332334                                        ./modules/Solverx/Solverx.cpp\
    333335                                        ./modules/Solverx/Solverx.h\

  • issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r10287 r10440  
    428428                case SpcDynamicEnum : return "SpcDynamic";
    429429                case IceVolumeEnum : return "IceVolume";
     430                case DiagnosticStokesBasalPenaltyEnum : return "DiagnosticStokesBasalPenalty";
    430431                default : return "unknown";
    431432

  • issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r9902 r10440  
    6565        parameters->AddObject(iomodel->CopyConstantObject(DiagnosticRiftPenaltyThresholdEnum));
    6666        parameters->AddObject(iomodel->CopyConstantObject(DiagnosticStokesreconditioningEnum));
     67        parameters->AddObject(iomodel->CopyConstantObject(DiagnosticStokesBasalPenaltyEnum));
    6768        parameters->AddObject(iomodel->CopyConstantObject(DiagnosticShelfDampeningEnum));
    6869        parameters->AddObject(iomodel->CopyConstantObject(DiagnosticViscosityOvershootEnum));

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

    r10126 r10440  
    2020        double  rho_ice;
    2121        double  stokesreconditioning;
    22         bool    isstokes,ismacayealpattyn;
     22        bool    isstokes,ismacayealpattyn,stokes_basal_penalty;
    2323   bool    spcpresent=false;
    2424        int Mx,Nx;
    2525        int My,Ny;
    2626        int Mz,Nz;
    27         double *spcvx           = NULL;
    28         double *spcvy           = NULL;
    29         double *spcvz           = NULL;
    30         double *nodeonmacayeal  = NULL;
    31         double *nodeonpattyn    = NULL;
    32         double *nodeonstokes    = NULL;
    33         double *vertices_type   = NULL;
    34         double *surface         = NULL;
    35         double *z               = NULL;
     27        double *spcvx          = NULL;
     28        double *spcvy          = NULL;
     29        double *spcvz          = NULL;
     30        double *nodeonmacayeal = NULL;
     31        double *nodeonpattyn   = NULL;
     32        double *nodeonstokes   = NULL;
     33        double *nodeonbed      = NULL;
     34        double *nodeonicesheet = NULL;
     35        double *vertices_type  = NULL;
     36        double *surface        = NULL;
     37        double *z              = NULL;
    3638        double *timesx=NULL;
    3739        double *timesy=NULL;
    … …  
    4042
    4143        /*Output*/
    42         Constraints* constraints = NULL;
    43         SpcStatic*    spcstatic  = NULL;
    44         int     node1,node2;
    45         int    dim;
    46         int    numberofvertices;
     44        Constraints *constraints      = NULL;
     45        SpcStatic   *spcstatic        = NULL;
     46        int          node1,node2;
     47        int          dim;
     48        int          numberofvertices;
    4749
    4850        /*Fetch parameters: */
    … …  
    5557        iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
    5658        iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
     59        iomodel->Constant(&stokes_basal_penalty,DiagnosticStokesBasalPenaltyEnum);
    5760
    5861        /*Recover pointer: */
    … …  
    7578        if(dim==3)iomodel->FetchData(&nodeonpattyn,NULL,NULL,FlowequationBorderpattynEnum);
    7679        if(dim==3)iomodel->FetchData(&nodeonstokes,NULL,NULL,FlowequationBorderstokesEnum);
     80        if(dim==3)iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     81        if(dim==3)iomodel->FetchData(&nodeonicesheet,NULL,NULL,MaskVertexongroundediceEnum);
    7782        iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
    7883        iomodel->FetchData(&surface,NULL,NULL,SurfaceEnum);
    … …  
    291296                                                }
    292297                                                xfree((void**)&values);
     298                                        }
     299
     300                                        /*Constraint at the bedrock interface (v.n ~ vz ~ 0)*/
     301                                        if(!stokes_basal_penalty){
     302                                                if(nodeonbed[i] && nodeonicesheet[i] && nodeonstokes[i]){
     303                                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,3,0.,DiagnosticHorizAnalysisEnum));
     304                                                        count++;
     305                                                }
    293306                                        }
    294307                                }
    … …  
    309322        xfree((void**)&nodeonpattyn);
    310323        xfree((void**)&nodeonstokes);
     324        xfree((void**)&nodeonicesheet);
     325        xfree((void**)&nodeonbed);
    311326        xfree((void**)&vertices_type);
    312327        xfree((void**)&surface);

  • issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

    r9733 r10440  
    2727        int dim;
    2828        int numberofvertices;
    29         bool ismacayealpattyn,isstokes;
     29        bool ismacayealpattyn,isstokes,stokes_basal_penalty;
    3030        int  numpenalties,numberofpressureloads,numrifts,numriftsegments;
    3131        double *pressureload   = NULL;
    … …  
    4545        iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
    4646        iomodel->Constant(&numrifts,RiftsNumriftsEnum);
     47        iomodel->Constant(&stokes_basal_penalty,DiagnosticStokesBasalPenaltyEnum);
    4748
    4849        /*Recover pointer: */
    … …  
    123124        xfree((void**)&pressureload);
    124125
    125         /*create penalties for nodes on the base of icesheet. We must have wb=ub*db/dx+vb*db/dy */
    126         iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
    127         iomodel->FetchData(&nodeonicesheet,NULL,NULL,MaskVertexongroundediceEnum);
    128         iomodel->FetchData(&nodeonstokes,NULL,NULL,FlowequationBorderstokesEnum);
    129         iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
    130         iomodel->FetchData(1,MeshElementsEnum);
    131         CreateSingleNodeToElementConnectivity(iomodel);
    132        
    133         for (i=0;i<numberofvertices;i++){
    134                 if(iomodel->my_vertices[i]==1 && iomodel->singlenodetoelementconnectivity[i]!=0){
     126        if(stokes_basal_penalty){
     127                /*create penalties for nodes on the base of icesheet. We must have wb=ub*db/dx+vb*db/dy */
     128                iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);
     129                iomodel->FetchData(&nodeonicesheet,NULL,NULL,MaskVertexongroundediceEnum);
     130                iomodel->FetchData(&nodeonstokes,NULL,NULL,FlowequationBorderstokesEnum);
     131                iomodel->FetchData(&vertices_type,NULL,NULL,FlowequationVertexEquationEnum);
     132                iomodel->FetchData(1,MeshElementsEnum);
     133                CreateSingleNodeToElementConnectivity(iomodel);
    135134
    136                         //if ((nodeonbed[i]) && (nodeonstokes[i])){
    137                         if ((nodeonbed[i]) && (nodeonicesheet[i]) && (nodeonstokes[i])){
    138                                
    139                                 loads->AddObject(new Pengrid(iomodel->loadcounter+count+1,i,iomodel,DiagnosticHorizAnalysisEnum));
    140                                 count++;
     135                for (i=0;i<numberofvertices;i++){
     136                        if(iomodel->my_vertices[i]==1 && iomodel->singlenodetoelementconnectivity[i]!=0){
     137
     138                                //if ((nodeonbed[i]) && (nodeonstokes[i])){
     139                                if ((nodeonbed[i]) && (nodeonicesheet[i]) && (nodeonstokes[i])){
     140
     141                                        loads->AddObject(new Pengrid(iomodel->loadcounter+count+1,i,iomodel,DiagnosticHorizAnalysisEnum));
     142                                        count++;
     143                                }
    141144                        }
    142145                }
     146                xfree((void**)&nodeonbed);
     147                xfree((void**)&nodeonstokes);
     148                xfree((void**)&nodeonicesheet);
     149                xfree((void**)&vertices_type);
     150                iomodel->DeleteData(1,MeshElementsEnum);
    143151        }
    144         xfree((void**)&nodeonbed);
    145         xfree((void**)&nodeonstokes);
    146         xfree((void**)&nodeonicesheet);
    147         xfree((void**)&vertices_type);
    148         iomodel->DeleteData(1,MeshElementsEnum);
    149152
    150153        /*Create Penpair for penalties: */

  • issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r10287 r10440  
    426426        else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
    427427        else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
     428        else if (strcmp(name,"DiagnosticStokesBasalPenalty")==0) return DiagnosticStokesBasalPenaltyEnum;
    428429        else _error_("Enum %s not found",name);
    429430

  • issm/trunk/src/c/modules/modules.h

    r10292 r10440  
    9696#include "./RequestedOutputsx/RequestedOutputsx.h"
    9797#include "./ResetConstraintsx/ResetConstraintsx.h"
     98#include "./ResetCoordinateSystemx/ResetCoordinateSystemx.h"
    9899#include "./Responsex/Responsex.h"
    99100#include "./RheologyBbarx/RheologyBbarx.h"

  • issm/trunk/src/c/objects/Elements/Element.h

    r10377 r10440  
    6868                virtual void   PotentialSheetUngrounding(Vec potential_sheet_ungrounding)=0;
    6969                virtual int    UpdatePotentialSheetUngrounding(double* potential_sheet_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf)=0;
     70                virtual void   ResetCoordinateSystem()=0;
    7071
    7172                #ifdef _HAVE_RESPONSES_

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r10434 r10440  
    23502350        }
    23512351
     2352}
     2353/*}}}*/
     2354/*FUNCTION Penta::ResetCoordinateSystem{{{1*/
     2355void  Penta::ResetCoordinateSystem(void){
     2356
     2357        int    approximation;
     2358        bool   stokes_basal_penalty;
     2359        double slopex[NUMVERTICES];
     2360        double slopey[NUMVERTICES];
     2361        double xz_plane[6];
     2362
     2363        /*Do not change CS if stokes_basal_penalty is activated*/
     2364        parameters->FindParam(&stokes_basal_penalty,DiagnosticStokesBasalPenaltyEnum);
     2365        if(stokes_basal_penalty) return;
     2366
     2367        /*For Stokes only: we want the CS to be tangential to the bedrock*/
     2368        inputs->GetInputValue(&approximation,ApproximationEnum);
     2369        if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return;
     2370
     2371        /*Get slope on each node*/
     2372        GetInputListOnVertices(&slopex[0],BedSlopeXEnum);
     2373        GetInputListOnVertices(&slopey[0],BedSlopeYEnum);
     2374
     2375        /*Loop over basal nodes (first 3) and update their CS*/
     2376        for(int i=0;i<NUMVERTICES2D;i++){
     2377
     2378                /*New X axis             New Z axis*/
     2379                xz_plane[0]=1.;          xz_plane[3]=-slopex[i]; 
     2380                xz_plane[1]=0.;          xz_plane[4]=-slopey[i]; 
     2381                xz_plane[2]=slopex[i];   xz_plane[5]=1.;         
     2382
     2383                XZvectorsToCoordinateSystem(&this->nodes[i]->coord_system[0][0],&xz_plane[0]);
     2384        }
    23522385}
    23532386/*}}}*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r10407 r10440  
    110110                void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
    111111                void   ProcessResultsUnits(void);
     112                void   ResetCoordinateSystem(void);
    112113                double SurfaceArea(void);
    113114                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);

  • issm/trunk/src/c/objects/Elements/Tria.h

    r10407 r10440  
    109109                void   PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes);
    110110                void   ProcessResultsUnits(void);
     111                void   ResetCoordinateSystem(void){_error_("not implemented yet");};
    111112                double SurfaceArea(void);
    112113                void   Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type);

  • issm/trunk/src/c/solutions/diagnostic_core.cpp

    r9880 r10440  
    4949        /*Compute slopes: */
    5050        if(ishutter) surfaceslope_core(femmodel);
    51         if(isstokes) bedslope_core(femmodel);
     51        if(isstokes){
     52                bedslope_core(femmodel);
     53                femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum);
     54                ResetCoordinateSystemx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     55        }
    5256
    5357        if(ishutter){

  • issm/trunk/src/m/classes/diagnostic.m

    r10398 r10440  
    1616                icefront                 = modelfield('default',NaN,'marshall',true,'preprocess','marshallicefront','format','DoubleMat','mattype',3);
    1717                maxiter                  = modelfield('default',0,'marshall',true,'format','Integer');
     18                stokes_basal_penalty     = modelfield('default',0,'marshall',true,'format','Boolean');
    1819                shelf_dampening          = modelfield('default',0,'marshall',true,'format','Integer');
    1920                vertex_pairing           = modelfield('default',NaN,'marshall',true,'format','DoubleMat','mattype',3);
    … …  
    7576                         %parameter is often used.
    7677                         obj.rift_penalty_lock=10;
     78
     79                         %penalty scheme for basal Dirichlet condition of Full-Stokes
     80                         %activated by default
     81                         obj.stokes_basal_penalty=1;
    7782
    7883                end % }}}
    … …  
    145150                        disp(sprintf('\n      %s','Other:'));
    146151                        fielddisplay(obj,'shelf_dampening','use dampening for floating ice ? Only for Stokes model');
     152                        fielddisplay(obj,'stokesreconditioning','multiplier for incompressibility equation. Only for Stokes model');
     153                        fielddisplay(obj,'stokes_basal_penalty','penalty based basal condition for grounded ice. Only for Stokes model');
    147154                        fielddisplay(obj,'referential','local referential');
    148155                        fielddisplay(obj,'requested_outputs','additional outputs requested');

  • issm/trunk/src/mex/Makefile.am

    r10299 r10440  
    6666                                Reducevectorgtof\
    6767                                Response\
     68                                ResetCoordinateSystem\
    6869                                Scotch\
    6970                                Solver\
    … …  
    302303                                                                        Response/Response.h
    303304
     305ResetCoordinateSystem_SOURCES = ResetCoordinateSystem/ResetCoordinateSystem.cpp\
     306                                                 ResetCoordinateSystem/ResetCoordinateSystem.h
     307
    304308Scotch_SOURCES = Scotch/Scotch.cpp\
    305309                          Scotch/Scotch.h
Note: See TracChangeset for help on using the changeset viewer.