Changeset 19127


Ignore:
Timestamp:
02/19/15 16:23:50 (10 years ago)
Author:
Mathieu Morlighem
Message:

BUG: LevelsetAnalysis.cpp was not updated properly

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r18301 r19127  
    1111#include "../solutionsequences/solutionsequences.h"
    1212
    13 int LevelsetAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
     13void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
     14        return;
     15}
     16/*}}}*/
     17void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
     18        return;
     19}/*}}}*/
     20void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
     21        int finiteelement=P1Enum;
     22        if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
     23        ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
     24        iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
     25}
     26/*}}}*/
     27int  LevelsetAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
    1428        return 1;
    1529}
    1630/*}}}*/
    17 void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    18         return;
    19 }
    20 /*}}}*/
    2131void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    22         int  finiteelement;
    2332
    2433        /*Finite element type*/
    25         finiteelement = P1Enum;
     34        int finiteelement = P1Enum;
    2635
    2736        /*Update elements: */
    … …  
    3847        iomodel->FetchDataToInput(elements,VxEnum);
    3948        iomodel->FetchDataToInput(elements,VyEnum);
    40         iomodel->FetchDataToInput(elements,MasstransportCalvingrateEnum);
    41 }
    42 /*}}}*/
    43 void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    44         int finiteelement=P1Enum;
    45         if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
    46         ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
    47         iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
    48 }
    49 /*}}}*/
    50 void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
     49
     50        /*Get calving parameters*/
     51        bool iscalving;
     52        int  calvinglaw;
     53        iomodel->Constant(&iscalving,TransientIscalvingEnum);
     54        if(iscalving){
     55                iomodel->Constant(&calvinglaw,CalvingLawEnum);
     56                iomodel->Constant(&iscalving,TransientIscalvingEnum);
     57                switch(calvinglaw){
     58                        case DefaultCalvingEnum:
     59                                iomodel->FetchDataToInput(elements,CalvingCalvingrateEnum);
     60                                iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
     61                                break;
     62                        case CalvingLevermannEnum:
     63                                iomodel->FetchDataToInput(elements,CalvinglevermannCoeffEnum);
     64                                iomodel->FetchDataToInput(elements,CalvinglevermannMeltingrateEnum);
     65                                break;
     66                        case CalvingPiEnum:
     67                                iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
     68                                iomodel->FetchDataToInput(elements,CalvingpiMeltingrateEnum);
     69                                break;
     70                        case CalvingDevEnum:
     71                                iomodel->FetchDataToInput(elements,CalvingpiCoeffEnum);
     72                                iomodel->FetchDataToInput(elements,CalvingMeltingrateEnum);
     73                                break;
     74                        default:
     75                                _error_("Calving law "<<EnumToStringx(calvinglaw)<<" not supported yet");
     76                }
     77        }
     78}
     79/*}}}*/
     80void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     81        parameters->AddObject(iomodel->CopyConstantObject(LevelsetStabilizationEnum));
    5182        return;
    5283}
    5384/*}}}*/
    54 void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    55         return;
    56 }/*}}}*/
    5785
    5886/*Finite element Analysis*/
    59 void LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
     87void           LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/
    6088
    6189        #if !defined(_DEVELOPMENT_)
    … …  
    6492
    6593        /*parameters: */
     94        int  stabilization;
    6695        bool save_results;
    6796        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     97        femmodel->parameters->FindParam(&stabilization,LevelsetStabilizationEnum);
    6898
    6999        /*activate formulation: */
    … …  
    71101
    72102        if(VerboseSolution()) _printf0_("call computational core:\n");
    73         solutionsequence_linear(femmodel);
     103        if(stabilization==4){
     104                solutionsequence_fct(femmodel);
     105        }
     106        else{
     107                solutionsequence_linear(femmodel);
     108        }
    74109
    75110        if(save_results){
    76111                if(VerboseSolution()) _printf0_("   saving results\n");
    77                 int outputs = MaskIceLevelsetEnum;
    78                 femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
     112                int outputs[2] = {MaskIceLevelsetEnum, CalvingCalvingrateEnum};
     113                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
    79114        }
    80115}/*}}}*/
    … …  
    97132        _error_("not implemented yet");
    98133}/*}}}*/
    99 void LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
     134void           LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     135        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     136         * For node i, Bi can be expressed in the actual coordinate system
     137         * by:
     138         *       Bi=[ N ]
     139         *          [ N ]
     140         * where N is the finiteelement function for node i.
     141         *
     142         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
     143         */
     144
     145        /*Fetch number of nodes for this finite element*/
     146        int numnodes = element->GetNumberOfNodes();
     147
     148        /*Get nodal functions*/
     149        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     150        element->NodalFunctions(basis,gauss);
     151
     152        /*Build B: */
     153        for(int i=0;i<numnodes;i++){
     154                B[numnodes*0+i] = basis[i];
     155                B[numnodes*1+i] = basis[i];
     156        }
     157
     158        /*Clean-up*/
     159        xDelete<IssmDouble>(basis);
     160}/*}}}*/
     161void           LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     162        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     163         * For node i, Bi' can be expressed in the actual coordinate system
     164         * by:
     165         *       Bi_prime=[ dN/dx ]
     166         *                [ dN/dy ]
     167         * where N is the finiteelement function for node i.
     168         *
     169         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     170         */
     171
     172        /*Fetch number of nodes for this finite element*/
     173        int numnodes = element->GetNumberOfNodes();
     174
     175        /*Get nodal functions derivatives*/
     176        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     177        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     178
     179        /*Build B': */
     180        for(int i=0;i<numnodes;i++){
     181                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
     182                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
     183        }
     184
     185        /*Clean-up*/
     186        xDelete<IssmDouble>(dbasis);
     187
     188}/*}}}*/
     189IssmDouble     LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
     190        // returns distance d of point q to straight going through points s0, s1
     191        // d=|a x b|/|b|
     192        // with a=q-s0, b=s1-s0
     193       
     194        /* Intermediaries */
     195        const int dim=2;
     196        int i;
     197        IssmDouble a[dim], b[dim];
     198        IssmDouble norm_b;
     199
     200        for(i=0;i<dim;i++){
     201                a[i]=q[i]-s0[i];
     202                b[i]=s1[i]-s0[i];
     203        }
     204       
     205        norm_b=0.;
     206        for(i=0;i<dim;i++)
     207                norm_b+=b[i]*b[i];
     208        norm_b=sqrt(norm_b);
     209        _assert_(norm_b>0.);
     210
     211        return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
     212}/*}}}*/
     213void           LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    100214        _error_("not implemented yet");
    101215}/*}}}*/
    102 void LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
     216void           LevelsetAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
    103217        _error_("Not implemented yet");
    104218}/*}}}*/
    105 void LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
     219void           LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    106220
    107221        int domaintype;
    … …  
    117231        }
    118232}/*}}}*/
    119 void LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    120         /*Default, do nothing*/
    121         return;
    122 }/*}}}*/
    123 void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    124         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    125          * For node i, Bi can be expressed in the actual coordinate system
    126          * by:
    127          *       Bi=[ N ]
    128          *          [ N ]
    129          * where N is the finiteelement function for node i.
    130          *
    131          * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
    132          */
    133 
    134         /*Fetch number of nodes for this finite element*/
    135         int numnodes = element->GetNumberOfNodes();
    136 
    137         /*Get nodal functions*/
    138         IssmDouble* basis=xNew<IssmDouble>(numnodes);
    139         element->NodalFunctions(basis,gauss);
    140 
    141         /*Build B: */
    142         for(int i=0;i<numnodes;i++){
    143                 B[numnodes*0+i] = basis[i];
    144                 B[numnodes*1+i] = basis[i];
    145         }
    146 
    147         /*Clean-up*/
    148         xDelete<IssmDouble>(basis);
    149 }/*}}}*/
    150 void LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    151         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    152          * For node i, Bi' can be expressed in the actual coordinate system
    153          * by:
    154          *       Bi_prime=[ dN/dx ]
    155          *                [ dN/dy ]
    156          * where N is the finiteelement function for node i.
    157          *
    158          * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
    159          */
    160 
    161         /*Fetch number of nodes for this finite element*/
    162         int numnodes = element->GetNumberOfNodes();
    163 
    164         /*Get nodal functions derivatives*/
    165         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
    166         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    167 
    168         /*Build B': */
    169         for(int i=0;i<numnodes;i++){
    170                 Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
    171                 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
    172         }
    173 
    174         /*Clean-up*/
    175         xDelete<IssmDouble>(dbasis);
    176 
    177 }/*}}}*/
    178 void LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
     233void           LevelsetAnalysis::SetDistanceOnIntersectedElements(FemModel* femmodel){/*{{{*/
    179234
    180235        /* Intermediaries */
    … …  
    187242
    188243        Vector<IssmDouble>* vec_dist_zerolevelset = NULL;
    189         GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexEnum);
     244        GetVectorFromInputsx(&vec_dist_zerolevelset, femmodel, MaskIceLevelsetEnum, VertexPIdEnum);
    190245       
    191246        /* set NaN on elements intersected by zero levelset */
    192247        for(i=0;i<femmodel->elements->Size();i++){
    193                 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     248                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    194249                if(element->IsZeroLevelset(MaskIceLevelsetEnum))
    195250                        for(k=0;k<element->GetNumberOfVertices();k++)
    … …  
    199254        /* set distance on elements intersected by zero levelset */
    200255        for(i=0;i<femmodel->elements->Size();i++){
    201                 element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
     256                element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    202257                if(element->IsZeroLevelset(MaskIceLevelsetEnum)){
    203258                        SetDistanceToZeroLevelsetElement(vec_dist_zerolevelset, element);
    … …  
    232287        delete dist_zerolevelset;
    233288}/*}}}*/
    234 void LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
     289void           LevelsetAnalysis::SetDistanceToZeroLevelsetElement(Vector<IssmDouble>* vec_signed_dist, Element* element){/*{{{*/
    235290
    236291        if(!element->IsZeroLevelset(MaskIceLevelsetEnum))
    … …  
    285340        xDelete<IssmDouble>(xyz_list_zero);
    286341}/*}}}*/
    287 IssmDouble LevelsetAnalysis::GetDistanceToStraight(IssmDouble* q, IssmDouble* s0, IssmDouble* s1){/*{{{*/
    288         // returns distance d of point q to straight going through points s0, s1
    289         // d=|a x b|/|b|
    290         // with a=q-s0, b=s1-s0
    291        
    292         /* Intermediaries */
    293         const int dim=2;
    294         int i;
    295         IssmDouble a[dim], b[dim];
    296         IssmDouble norm_b;
    297 
    298         for(i=0;i<dim;i++){
    299                 a[i]=q[i]-s0[i];
    300                 b[i]=s1[i]-s0[i];
    301         }
    302        
    303         norm_b=0.;
    304         for(i=0;i<dim;i++)
    305                 norm_b+=b[i]*b[i];
    306         norm_b=sqrt(norm_b);
    307         _assert_(norm_b>0.);
    308 
    309         return fabs(a[0]*b[1]-a[1]*b[0])/norm_b;
    310 }/*}}}*/
     342void           LevelsetAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
     343        /*Default, do nothing*/
     344        return;
     345}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.