Changeset 17208


Ignore:
Timestamp:
02/04/14 14:16:16 (11 years ago)
Author:
jbondzio
Message:

CHG: setting of BCs for LSM, minor bugs

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

Legend:

Unmodified
Added
Removed

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

    r17159 r17208  
    6767        bool save_results;
    6868
    69         /* extrapolate */
    70         Analysis* analysis = new ExtrapolationAnalysis();
    71         femmodel->parameters->SetParam(VxEnum,ExtrapolationVariableEnum);
    72         analysis->Core(femmodel);
    73         delete analysis;
    74 
    7569        /*activate formulation: */
    7670        femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
    7771
     72        if(VerboseSolution()) _printf0_("call computational core:\n");
     73        solutionsequence_linear(femmodel);
     74
    7875        /*recover parameters: */
    7976        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    80 
    81         if(VerboseSolution()) _printf0_("call computational core:\n");
    82         solutionsequence_linear(femmodel);
    83 
    8477        if(save_results){
    8578                if(VerboseSolution()) _printf0_("   saving results\n");
    … …  
    324317}/*}}}*/
    325318
     319/* Update of constraints */
     320void LevelsetAnalysis::UpdateNoIceConstraints(FemModel* femmodel){/*{{{*/
     321
     322        IssmDouble* mask_ice = GetMaskOfIce(femmodel->elements, femmodel->nodes);
     323
     324        for (int i=0;i<femmodel->nodes->Size();i++){
     325                Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
     326                if(node->InAnalysis(LevelsetAnalysisEnum)){
     327                        if(mask_ice[node->Sid()]==1.){//FIXME: what should be done with actual spcs to ice model?
     328//                              node->DofInFSet(0); /*remove spc*/
     329                        }
     330                        else{
     331                                IssmDouble defval=0.;
     332                                node->ApplyConstraint(1,defval); /*apply spc*/
     333                        }
     334                }
     335        }
     336}/*}}}*/
     337IssmDouble* LevelsetAnalysis::GetMaskOfIce(Elements* elements, Nodes* nodes){/*{{{*/
     338
     339        int                 i;
     340        IssmDouble*         mask_ice      = NULL;
     341        Vector<IssmDouble>* vec_mask_ice  = NULL;
     342        Element*            element       = NULL;
     343
     344        /*Initialize vector with number of vertices*/
     345        IssmDouble numnodes=nodes->NumberOfNodes(LevelsetAnalysisEnum);
     346        vec_mask_ice=new Vector<IssmDouble>(numnodes); //nodes at ice front that have ice at next time step
     347        for(i=0;i<numnodes;i++)
     348                vec_mask_ice[i]=0.;
     349        /*Fill vector vertices that have no contact to ice: */
     350        for(i=0;i<elements->Size();i++){
     351                element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
     352                SetMaskOfIceElement(vec_mask_ice, element);
     353        }
     354
     355        /*Assemble vector and serialize */
     356        vec_mask_ice->Assemble();
     357        mask_ice=vec_mask_ice->ToMPISerial();
     358        delete vec_mask_ice;
     359        return mask_ice;
     360
     361}/*}}}*/
     362void LevelsetAnalysis::SetMaskOfIceElement(Vector<IssmDouble>* vec_mask_ice, Element* element){/*{{{*/
     363
     364        /* Intermediaries */
     365        int numnodes = element->GetNumberOfNodes();
     366       
     367        if(element->IsIceInElement()){
     368                for(int i = 0;i<numnodes;i++){
     369                        vec_mask_ice->SetValue(element->nodes[i]->Sid(),1.,INS_VAL);
     370                }
     371        }
     372}/*}}}*/

  • issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.h

    r17098 r17208  
    3131        void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    3232
     33        /* Updating constraints */
     34        void UpdateNoIceConstraints(FemModel* femmodel);
     35        IssmDouble* GetMaskOfIce(Elements* elements, Nodes* nodes);
     36        void SetMaskOfIceElement(Vector<IssmDouble>* vec_mask_ice, Element* element);
    3337};
    3438#endif

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r17193 r17208  
    13681368
    13691369        /*If no front, return NULL*/
    1370         if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
     1370        if(!element->IsIcefront()) return NULL;
    13711371
    13721372        /*Intermediaries*/
    … …  
    13911391        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    13921392        element->GetVerticesCoordinates(&xyz_list);
    1393         element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     1393
     1394//      element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     1395        element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     1396
    13941397        element->NormalSection(&normal[0],xyz_list_front);
     1398//      element->GetNormalFromLSF(&normal[0]);
    13951399
    13961400        /*Start looping on Gaussian points*/
    … …  
    17811785
    17821786        /*If no front, return NULL*/
    1783         if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
     1787        if(!element->IsIcefront()) return NULL;
    17841788
    17851789        /*Intermediaries*/
    … …  
    22042208
    22052209        /*If no front, return NULL*/
    2206         if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
     2210        if(!element->IsIcefront()) return NULL;
    22072211
    22082212        /*Intermediaries*/
    … …  
    29702974
    29712975        /*If no front, return NULL*/
    2972         if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
     2976        if(!element->IsIcefront()) return NULL;
    29732977
    29742978        /*Intermediaries*/

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

    r17194 r17208  
    220220                virtual int    PressureInterpolation()=0;
    221221                virtual bool   IsZeroLevelset(int levelset_enum)=0;
     222                virtual bool   IsIcefront(void)=0;
    222223                virtual void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
    223224                virtual void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0;

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

    r17194 r17208  
    114114                int    PressureInterpolation();
    115115                bool   IsZeroLevelset(int levelset_enum);
     116                bool   IsIcefront(void){_error_("not implemented yet");};
    116117                void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    117118                void   GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};

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

    r17194 r17208  
    159159                void        ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
    160160                bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
     161                bool            IsIcefront(void){_error_("not implemented yet");};
    161162                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
    162163                void            GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented yet");};

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

    r17199 r17208  
    909909                }
    910910        }
    911 
     911        _assert_(nrfrontnodes==2);
     912       
     913        /* arrange order of frontnodes such that they are oriented counterclockwise */
     914        if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=2){
     915                int index=indicesfront[0];
     916                indicesfront[0]=indicesfront[1];
     917                indicesfront[1]=index;
     918        }
    912919        IssmDouble* xyz_front = xNew<IssmDouble>(3*nrfrontnodes);
    913920        /* Return nodes */
    914921        for(i=0;i<nrfrontnodes;i++){
    915922                for(dir=0;dir<3;dir++){
    916                         xyz_front[3*i+dir]=xyz_list[indicesfront[i]+dir];
    917                 }
    918         }
     923                        xyz_front[3*i+dir]=xyz_list[3*indicesfront[i]+dir];
     924                }
     925        }
     926
     927//      for(i=0;i<nrfrontnodes;i++){
     928//              _printf0_("coords frontnode " << i << " :[");
     929//              for(dir=0;dir<3;dir++){
     930//                      xyz_front[3*i+dir]=xyz_list[indicesfront[i]+dir];
     931//                      _printf0_(xyz_front[3*i+dir] << "; ");
     932//              }
     933//              _printf0_("]\n");
     934//      }
    919935
    920936        *pxyz_front=xyz_front;
    … …  
    922938        xDelete<int>(indicesfront);
    923939}/*}}}*/
    924 void  Tria::GetNormalFromLSF(IssmDouble *pnormal){/*{{{*/
     940void  Tria::GetNormalFromLSF(IssmDouble *normal){/*{{{*/
    925941
    926942        /* Intermediaries */
    … …  
    929945        IssmDouble* xyz_list = NULL;
    930946        IssmDouble  dlevelset[dim], norm_dlevelset;
    931         IssmDouble  normal[dim]={0.};
    932947
    933948        /*Retrieve all inputs and parameters*/
    … …  
    936951       
    937952        counter=0;
     953        for(i=0;i<dim;i++) normal[i]=0.;
     954
    938955        Gauss* gauss = this->NewGauss(2);
    939956        for(int ig=gauss->begin();ig<gauss->end();ig++){
    … …  
    947964        }
    948965        _assert_(counter>0);
    949         for(i=0;i<dim;i++) normal[i]/counter;
    950        
    951         pnormal=&normal[0];
     966        for(i=0;i<dim;i++) normal[i]/=IssmDouble(counter);
    952967
    953968        delete gauss;
    … …  
    25862601        GetInputListOnVertices(&ls[0],levelset_enum);
    25872602
    2588         /*If the level set is awlays <0, there is no ice front here*/
     2603        /* If levelset function changes sign, or is zero on at least two vertices, there is a zero level set here */
    25892604        iszerols= false;
    25902605        if(IsIceInElement()){
    … …  
    25972612}
    25982613/*}}}*/
     2614/*FUNCTION Tria::IsIcefront{{{*/
     2615bool Tria::IsIcefront(void){
     2616
     2617        bool isicefront;
     2618        int i,nrice;
     2619        IssmDouble ls[NUMVERTICES];
     2620
     2621        /*Retrieve all inputs and parameters*/
     2622        GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
     2623
     2624        /* If only one vertex has ice, there is an ice front here */
     2625        isicefront= false;
     2626        if(IsIceInElement()){
     2627                nrice=0;       
     2628                for(i=0;i<NUMVERTICES;i++)
     2629                        if(ls[i]<0.) nrice++;
     2630                if(nrice==1) isicefront= true;
     2631        }
     2632        return isicefront;
     2633}
     2634/*}}}*/
     2635
    25992636
    26002637#ifdef _HAVE_RESPONSES_

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

    r17194 r17208  
    132132                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    133133                void        GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
    134             void        GetNormalFromLSF(IssmDouble *pnormal);
     134                void        GetNormalFromLSF(IssmDouble *normal);
    135135                bool        IsZeroLevelset(int levelset_enum);
     136                bool        IsIcefront(void);
    136137
    137138                #ifdef _HAVE_RESPONSES_

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r17139 r17208  
    162162                if(islevelset){
    163163                        if(VerboseSolution()) _printf0_("   computing movement of ice boundaries\n");
    164                         analysis = new LevelsetAnalysis();
    165                         analysis->Core(femmodel);
    166                         delete analysis;
     164
     165                        /* extrapolate */
     166                        Analysis* extanalysis = new ExtrapolationAnalysis();
     167                        int vars[2] = {VxEnum, VyEnum};
     168                        for(int iv=0; i<2;i++){
     169                                femmodel->parameters->SetParam(vars[i],ExtrapolationVariableEnum);
     170                                extanalysis->Core(femmodel);
     171                        }
     172                        delete extanalysis;
     173
     174                        LevelsetAnalysis* lsanalysis = new LevelsetAnalysis();
     175                        /* solve level-set equation */
     176                        lsanalysis->Core(femmodel);
     177                        /* update spc boundary conditions for ice model */
     178//                      lsanalysis->UpdateNoIceConstraints(femmodel);
     179                        delete lsanalysis;
    167180                }
    168181
Note: See TracChangeset for help on using the changeset viewer.