Ignore:
Timestamp:
11/01/19 12:01:57 (5 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 24310

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        55*.obj
        66*.exe
         7*.mod
        78appscan.*
        89ouncemake*
        … …  
        1415probe.results
        1516stXXXX*
        16 *.deps
        17 *.dirstamp
         17.deps
         18.dirstamp
        1819.libs
        1920issm
        … …  
        2122issm_slr
        2223issm_ocean
        23 lnb_param.mod
        24 lovenb_sub.mod
        25 model.mod
        26 util.mod
         24issm_dakota

  • issm/trunk/src/c/classes/Loads/Moulin.cpp

    r23189 r24313  
    2323        this->parameters=NULL;
    2424        this->hnode=NULL;
     25        this->hvertex=NULL;
    2526        this->node=NULL;
    2627        this->helement=NULL;
    2728        this->element=NULL;
    28         this->hmatpar=NULL;
    29         this->matpar=NULL;
    30 }
    31 /*}}}*/
    32 Moulin::Moulin(int id, int index, IoModel* iomodel, int in_analysis_type){ //i is the element index/*{{{*/
     29}
     30/*}}}*/
     31Moulin::Moulin(int id, int index, IoModel* iomodel){ /*{{{*/
    3332
    3433        int pengrid_node_id;
    35         int pengrid_matpar_id;
    3634        int pengrid_element_id;
    3735
    … …  
    4341        /*id: */
    4442        this->id=id;
    45         this->analysis_type=in_analysis_type;
    4643
    4744        /*hooks: */
    48         pengrid_node_id=iomodel->nodecounter+index+1;
     45        pengrid_node_id=index+1;
    4946        pengrid_element_id=iomodel->singlenodetoelementconnectivity[index];
    5047        _assert_(pengrid_element_id);
    51         pengrid_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
    5248
    5349        this->hnode=new Hook(&pengrid_node_id,1);
     50        this->hvertex=new Hook(&pengrid_node_id,1);
    5451        this->helement=new Hook(&pengrid_element_id,1);
    55         this->hmatpar=new Hook(&pengrid_matpar_id,1);
    5652
    5753        //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
    5854        this->parameters=NULL;
    5955        this->node=NULL;
     56        this->vertex=NULL;
    6057        this->element=NULL;
    61         this->matpar=NULL;
    6258}
    6359/*}}}*/
    6460Moulin::~Moulin(){/*{{{*/
    6561        delete hnode;
     62        delete hvertex;
    6663        delete helement;
    67         delete hmatpar;
    6864        return;
    6965}
    … …  
    7975        /*copy fields: */
    8076        pengrid->id=this->id;
    81         pengrid->analysis_type=this->analysis_type;
    8277
    8378        /*point parameters: */
    … …  
    8681        /*now deal with hooks and objects: */
    8782        pengrid->hnode=(Hook*)this->hnode->copy();
    88         pengrid->hmatpar=(Hook*)this->hmatpar->copy();
     83        pengrid->hvertex=(Hook*)this->hvertex->copy();
    8984        pengrid->helement=(Hook*)this->helement->copy();
    9085
    9186        /*corresponding fields*/
    9287        pengrid->node  =(Node*)pengrid->hnode->delivers();
    93         pengrid->matpar =(Matpar*)pengrid->hmatpar->delivers();
     88        pengrid->vertex=(Vertex*)pengrid->hvertex->delivers();
    9489        pengrid->element=(Element*)pengrid->helement->delivers();
    9590
    … …  
    10196        _printf_("Moulin:\n");
    10297        _printf_("   id: " << id << "\n");
    103         _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
    10498        hnode->DeepEcho();
     99        hvertex->DeepEcho();
    105100        helement->DeepEcho();
    106         hmatpar->DeepEcho();
    107101        _printf_("   parameters\n");
    108102        parameters->DeepEcho();
    … …  
    122116        MARSHALLING_ENUM(MoulinEnum);
    123117        MARSHALLING(id);
    124         MARSHALLING(analysis_type);
    125118
    126119        if(marshall_direction==MARSHALLING_BACKWARD){
    127120                this->hnode      = new Hook();
     121                this->hvertex      = new Hook();
    128122                this->helement   = new Hook();
    129                 this->hmatpar    = new Hook();
    130123        }
    131124
    132125        this->hnode->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
     126        this->hvertex->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    133127        this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    134         this->hmatpar->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    135128
    136129        /*corresponding fields*/
    137130        node   =(Node*)this->hnode->delivers();
    138         matpar =(Matpar*)this->hmatpar->delivers();
     131        vertex =(Vertex*)this->hvertex->delivers();
    139132        element=(Element*)this->helement->delivers();
    140133}
    … …  
    152145         * datasets, using internal ids and offsets hidden in hooks: */
    153146        hnode->configure(nodesin);
     147        hvertex->configure(verticesin);
    154148        helement->configure(elementsin);
    155         hmatpar->configure(materialsin);
    156149
    157150        /*Get corresponding fields*/
    158151        node=(Node*)hnode->delivers();
     152        vertex=(Vertex*)hvertex->delivers();
    159153        element=(Element*)helement->delivers();
    160         matpar=(Matpar*)hmatpar->delivers();
    161154
    162155        /*point parameters to real dataset: */
    … …  
    166159void  Moulin::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){/*{{{*/
    167160
    168         /*No loads applied, do nothing: */
    169         return;
     161        /*recover some parameters*/
     162        ElementMatrix* Ke=NULL;
     163        int analysis_type;
     164        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     165
     166        switch(analysis_type){
     167                case HydrologyGlaDSAnalysisEnum:
     168                        Ke = this->CreateKMatrixHydrologyGlaDS();
     169                        break;
     170                case HydrologyShaktiAnalysisEnum:
     171                        /*do nothing: */
     172                        return;
     173                case HydrologyDCInefficientAnalysisEnum:
     174                        /*do nothing: */
     175                        return;
     176                case HydrologyDCEfficientAnalysisEnum:
     177                        /*do nothing: */
     178                        return;
     179                default:
     180                        _error_("Don't know why we should be here");
     181
     182        }
     183        /*Add to global matrix*/
     184        if(Ke){
     185                Ke->AddToGlobal(Kff,Kfs);
     186                delete Ke;
     187        }
    170188
    171189}
    … …  
    178196
    179197        switch(analysis_type){
    180 
    181         case HydrologyShaktiAnalysisEnum:
    182                 pe = this->CreatePVectorHydrologyShakti();
    183                 break;
    184         case HydrologyDCInefficientAnalysisEnum:
    185                 pe = CreatePVectorHydrologyDCInefficient();
    186                 break;
    187         case HydrologyDCEfficientAnalysisEnum:
    188                 pe = CreatePVectorHydrologyDCEfficient();
    189                 break;
    190         default:
    191                 _error_("Don't know why we should be here");
    192                 /*No loads applied, do nothing: */
    193                 return;
     198                case HydrologyGlaDSAnalysisEnum:
     199                        pe = this->CreatePVectorHydrologyGlaDS();
     200                        break;
     201                case HydrologyShaktiAnalysisEnum:
     202                        pe = this->CreatePVectorHydrologyShakti();
     203                        break;
     204                case HydrologyDCInefficientAnalysisEnum:
     205                        pe = CreatePVectorHydrologyDCInefficient();
     206                        break;
     207                case HydrologyDCEfficientAnalysisEnum:
     208                        pe = CreatePVectorHydrologyDCEfficient();
     209                        break;
     210                default:
     211                        _error_("Don't know why we should be here");
     212                        /*No loads applied, do nothing: */
     213                        return;
    194214        }
    195215        if(pe){
    … …  
    221241}
    222242/*}}}*/
    223 bool  Moulin::InAnalysis(int in_analysis_type){/*{{{*/
    224         if (in_analysis_type==this->analysis_type)return true;
    225         else return false;
    226 }
    227 /*}}}*/
    228243bool  Moulin::IsPenalty(void){/*{{{*/
    229         return true;
     244        return false;
    230245}
    231246/*}}}*/
    … …  
    245260        this->node=NULL;
    246261        this->element=NULL;
    247         this->matpar=NULL;
    248262        this->parameters=NULL;
    249263
    250264        /*Get Element type*/
    251265        this->hnode->reset();
     266        this->hvertex->reset();
    252267        this->helement->reset();
    253         this->hmatpar->reset();
    254268
    255269}
    … …  
    277291                switch(set2_enum){
    278292                        case FsetEnum:
    279                                 if(node->indexing.fsize){
     293                                if(node->fsize){
    280294                                        if(this->node->IsClone())
    281295                                         o_nz += 1;
    … …  
    285299                                break;
    286300                        case GsetEnum:
    287                                 if(node->indexing.gsize){
     301                                if(node->gsize){
    288302                                        if(this->node->IsClone())
    289303                                         o_nz += 1;
    … …  
    293307                                break;
    294308                        case SsetEnum:
    295                                 if(node->indexing.ssize){
     309                                if(node->ssize){
    296310                                        if(this->node->IsClone())
    297311                                         o_nz += 1;
    … …  
    337351/*}}}*/
    338352
     353ElementMatrix* Moulin::CreateKMatrixHydrologyGlaDS(void){/*{{{*/
     354
     355        /*If this node is not the master node (belongs to another partition of the
     356         * mesh), don't add the moulin input a second time*/
     357        if(node->IsClone()) return NULL;
     358
     359        /*Initialize Element matrix*/
     360        ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters);
     361
     362        /*Get all inputs and parameters*/
     363        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
     364        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     365        IssmDouble g         = element->FindParam(ConstantsGEnum);
     366        IssmDouble Am        = 0.; //For now...
     367
     368        /*Load vector*/
     369        Ke->values[0] = - Am/(rho_water*g)/dt;
     370
     371        /*Clean up and return*/
     372        return Ke;
     373}
     374/*}}}*/
     375ElementVector* Moulin::CreatePVectorHydrologyGlaDS(void){/*{{{*/
     376
     377        /*If this node is not the master node (belongs to another partition of the
     378         * mesh), don't add the moulin input a second time*/
     379        if(node->IsClone()) return NULL;
     380
     381        /*Initialize Element vector*/
     382        ElementVector* pe=new ElementVector(&node,1,this->parameters);
     383
     384        /*Get all inputs and parameters*/
     385        IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
     386        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
     387        IssmDouble g         = element->FindParam(ConstantsGEnum);
     388        IssmDouble Am        = 0.; //For now...
     389
     390        /*Get hydraulic potential*/
     391        IssmDouble phi_old,moulin_load;
     392        element->GetInputValue(&phi_old,node,HydraulicPotentialOldEnum);
     393        element->GetInputValue(&moulin_load,node,HydrologyMoulinInputEnum);
     394
     395        pe->values[0] = moulin_load -Am/(rho_water*g) * phi_old/dt;
     396
     397        /*Clean up and return*/
     398        return pe;
     399}
     400/*}}}*/
    339401ElementVector* Moulin::CreatePVectorHydrologyShakti(void){/*{{{*/
    340402
Note: See TracChangeset for help on using the changeset viewer.