Changeset 17098


Ignore:
Timestamp:
01/10/14 16:34:01 (11 years ago)
Author:
jbondzio
Message:

CHG: LevelsetAnalysis::CreateKMatrix and CreatePVector filled, required enums added

Location:
issm/trunk-jpl/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r17076 r17098  
    77
    88int LevelsetAnalysis::DofsPerNode(int** doflist,int meshtype,int approximation){/*{{{*/
    9         _error_("not implemented yet");
     9        return 1;
    1010}
    1111/*}}}*/
    1212void LevelsetAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    13         _error_("not implemented yet");
     13        // parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
    1414}
    1515/*}}}*/
    1616void LevelsetAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    17         _error_("not implemented yet");
     17        int    stabilization,finiteelement;
     18
     19        /*Finite element type*/
     20        finiteelement = P1Enum;
     21
     22        /*Update elements: */
     23        int counter=0;
     24        for(int i=0;i<iomodel->numberofelements;i++){
     25                if(iomodel->my_elements[i]){
     26                        Element* element=(Element*)elements->GetObjectByOffset(counter);
     27                        element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
     28                        counter++;
     29                }
     30        }
     31
     32        iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
     33        iomodel->FetchDataToInput(elements,VxEnum);
     34        iomodel->FetchDataToInput(elements,VyEnum);
     35       
    1836}
    1937/*}}}*/
    2038void LevelsetAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    21         _error_("not implemented yet");
     39        int finiteelement=P1Enum;
     40        ::CreateNodes(nodes,iomodel,LevelsetAnalysisEnum,finiteelement);
    2241}
    2342/*}}}*/
    2443void LevelsetAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    25         _error_("not implemented yet");
     44
     45        /*Intermediary*/
     46
     47        /*Fetch parameters: */
     48
    2649}
    2750/*}}}*/
    2851void LevelsetAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
    29         _error_("not implemented yet");
     52       
     53        /*Intermediary*/
     54
     55        /*Fetch parameters: */
     56
    3057}/*}}}*/
    3158
    … …  
    3562}/*}}}*/
    3663ElementVector* LevelsetAnalysis::CreateDVector(Element* element){/*{{{*/
    37         _error_("not implemented yet");
     64        /*Default, return NULL*/
     65        return NULL;
    3866}/*}}}*/
    3967ElementMatrix* LevelsetAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/
     68        /* Jacobian required for the Newton solver */
    4069        _error_("not implemented yet");
    4170}/*}}}*/
    4271ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/
    43         _error_("not implemented yet");
     72        /*Intermediaries */
     73
     74        int dim = 2; // solve for LSF in horizontal plane only
     75        IssmDouble kappa;
     76        IssmDouble Jdet, dt, D_scalar;
     77        IssmDouble u,v,um,vm,ub,vb,vx,vy;
     78        IssmDouble* xyz_list = NULL;
     79
     80        /*Fetch number of nodes and dof for this finite element*/
     81        int numnodes    = element->GetNumberOfNodes();
     82
     83        /*Initialize Element vector and other vectors*/
     84        ElementMatrix* Ke       = element->NewElementMatrix();
     85        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
     86        IssmDouble*    B        = xNew<IssmDouble>(dim*numnodes);
     87        IssmDouble*    Bprime   = xNew<IssmDouble>(dim*numnodes);
     88        IssmDouble     D[dim][dim];
     89
     90        /*Retrieve all inputs and parameters*/
     91        element->GetVerticesCoordinates(&xyz_list);
     92        element->FindParam(&dt,TimesteppingTimeStepEnum);
     93        Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
     94        Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
     95        Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
     96        Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
     97       
     98        /* Start  looping on the number of gaussian points: */
     99        Gauss* gauss=element->NewGauss(2);
     100        for(int ig=gauss->begin();ig<gauss->end();ig++){
     101                gauss->GaussPoint(ig);
     102
     103                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     104                D_scalar=gauss->weight*Jdet;
     105
     106                /* Transient */
     107                if(dt!=0.){
     108                        element->NodalFunctions(basis,gauss);
     109                        TripleMultiply(basis,numnodes,1,0,
     110                                                &D_scalar,1,1,0,
     111                                                basis,1,numnodes,0,
     112                                                &Ke->values[0],1);
     113                        D_scalar*=dt;
     114                }
     115
     116                /* Advection */
     117                GetB(B,element,xyz_list,gauss);
     118                GetBprime(Bprime,element,xyz_list,gauss);
     119                vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
     120                vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
     121                ub=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here
     122                vb=0.;
     123                D[0][0]=D_scalar*(vx+ub);
     124                D[0][1]=0.;
     125                D[1][0]=0.;
     126                D[1][1]=D_scalar*(vy+vb);
     127                TripleMultiply(B,dim,numnodes,1,
     128                                        &D[0][0],dim,dim,0,
     129                                        Bprime,dim,numnodes,0,
     130                                        &Ke->values[0],1);
     131
     132                /* Artificial Diffusion */
     133                kappa=0.; //FIXME: insert suitable value for kappa
     134                GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed?
     135                D[0][0]=D_scalar*kappa;
     136                D[1][1]=D_scalar*kappa;
     137                TripleMultiply(Bprime,dim,numnodes,1,
     138                                        &D[0][0],dim,dim,0,
     139                                        Bprime,dim,numnodes,0,
     140                                        &Ke->values[0],1);
     141        }
     142
     143        /*Clean up and return*/
     144        xDelete<IssmDouble>(xyz_list);
     145        xDelete<IssmDouble>(basis);
     146        xDelete<IssmDouble>(B);
     147        xDelete<IssmDouble>(Bprime);
     148        delete gauss;
     149        return Ke;
    44150}/*}}}*/
    45151ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
    46         _error_("not implemented yet");
     152
     153        /*Intermediaries */
     154        int i, ig;
     155        IssmDouble  Jdet,dt;
     156        IssmDouble  phi;
     157        IssmDouble* xyz_list = NULL;
     158       
     159        /*Fetch number of nodes and dof for this finite element*/
     160        int numnodes = element->GetNumberOfNodes();
     161
     162        /*Initialize Element vector*/
     163        ElementVector* pe    = element->NewElementVector();
     164        element->FindParam(&dt,TimesteppingTimeStepEnum);
     165       
     166        if(dt!=0.){
     167                /*Initialize basis vector*/
     168                IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     169
     170                /*Retrieve all inputs and parameters*/
     171                element->GetVerticesCoordinates(&xyz_list);
     172                Input* levelset_input     = element->GetInput(MaskIceLevelsetEnum);                    _assert_(levelset_input);
     173
     174                /* Start  looping on the number of gaussian points: */
     175                Gauss* gauss=element->NewGauss(2);
     176                for(ig=gauss->begin();ig<gauss->end();ig++){
     177                        gauss->GaussPoint(ig);
     178
     179                        element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     180                        element->NodalFunctions(basis,gauss);
     181                        levelset_input->GetInputValue(&phi,gauss);
     182                        for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*phi*basis[i];
     183                }
     184
     185                /*Clean up and return*/
     186                xDelete<IssmDouble>(xyz_list);
     187                xDelete<IssmDouble>(basis);
     188                delete gauss;
     189        }
     190        else{for(i=0;i<numnodes;i++) pe->values[i]=0.;}
     191        return pe;
    47192}/*}}}*/
    48193void LevelsetAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    … …  
    52197        _error_("not implemented yet");
    53198}/*}}}*/
     199void LevelsetAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     200        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
     201         * For node i, Bi can be expressed in the actual coordinate system
     202         * by:
     203         *       Bi=[ N ]
     204         *          [ N ]
     205         * where N is the finiteelement function for node i.
     206         *
     207         * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
     208         */
     209
     210        /*Fetch number of nodes for this finite element*/
     211        int numnodes = element->GetNumberOfNodes();
     212
     213        /*Get nodal functions*/
     214        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     215        element->NodalFunctions(basis,gauss);
     216
     217        /*Build B: */
     218        for(int i=0;i<numnodes;i++){
     219                B[numnodes*0+i] = basis[i];
     220                B[numnodes*1+i] = basis[i];
     221        }
     222
     223        /*Clean-up*/
     224        xDelete<IssmDouble>(basis);
     225}/*}}}*/
     226void LevelsetAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     227        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
     228         * For node i, Bi' can be expressed in the actual coordinate system
     229         * by:
     230         *       Bi_prime=[ dN/dx ]
     231         *                [ dN/dy ]
     232         * where N is the finiteelement function for node i.
     233         *
     234         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
     235         */
     236
     237        /*Fetch number of nodes for this finite element*/
     238        int numnodes = element->GetNumberOfNodes();
     239
     240        /*Get nodal functions derivatives*/
     241        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     242        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     243
     244        /*Build B': */
     245        for(int i=0;i<numnodes;i++){
     246                Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
     247                Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
     248        }
     249
     250        /*Clean-up*/
     251        xDelete<IssmDouble>(dbasis);
     252
     253}/*}}}*/
     254

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

    r17076 r17098  
    2828        void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2929        void InputUpdateFromSolution(IssmDouble* solution,Element* element);
     30        void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     31        void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     32
    3033};
    3134#endif

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r17085 r17098  
    273273        TransientIsgiaEnum,
    274274        TransientIsdamageEnum,
    275         TransientIslevelsetEnum,
    276275        TransientNumRequestedOutputsEnum,
    277276        TransientRequestedOutputsEnum,
    … …  
    321320        DamageEvolutionSolutionEnum,
    322321        DamageEvolutionAnalysisEnum,
    323         LevelsetAnalysisEnum,
    324322        StressbalanceAnalysisEnum,
    325323        StressbalanceSIAAnalysisEnum,
    … …  
    674672        LliboutryDuvalEnum,
    675673        /*}}}*/
     674        /*Levelset related enums (will be moved to appropriate place when finished){{{*/
     675        LevelsetAnalysisEnum,
     676        TransientIslevelsetEnum,
     677        /*}}}*/
    676678        MaximumNumberOfDefinitionsEnum
    677679};

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r17085 r17098  
    279279                case TransientIsgiaEnum : return "TransientIsgia";
    280280                case TransientIsdamageEnum : return "TransientIsdamage";
    281                 case TransientIslevelsetEnum : return "TransientIslevelset";
    282281                case TransientNumRequestedOutputsEnum : return "TransientNumRequestedOutputs";
    283282                case TransientRequestedOutputsEnum : return "TransientRequestedOutputs";
    … …  
    322321                case DamageEvolutionSolutionEnum : return "DamageEvolutionSolution";
    323322                case DamageEvolutionAnalysisEnum : return "DamageEvolutionAnalysis";
    324                 case LevelsetAnalysisEnum : return "LevelsetAnalysis";
    325323                case StressbalanceAnalysisEnum : return "StressbalanceAnalysis";
    326324                case StressbalanceSIAAnalysisEnum : return "StressbalanceSIAAnalysis";
    … …  
    634632                case ArrheniusEnum : return "Arrhenius";
    635633                case LliboutryDuvalEnum : return "LliboutryDuval";
     634                case LevelsetAnalysisEnum : return "LevelsetAnalysis";
     635                case TransientIslevelsetEnum : return "TransientIslevelset";
    636636                case MaximumNumberOfDefinitionsEnum : return "MaximumNumberOfDefinitions";
    637637                default : return "unknown";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r17085 r17098  
    285285              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
    286286              else if (strcmp(name,"TransientIsdamage")==0) return TransientIsdamageEnum;
    287               else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
    288287              else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum;
    289288              else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum;
    … …  
    328327              else if (strcmp(name,"DamageEvolutionSolution")==0) return DamageEvolutionSolutionEnum;
    329328              else if (strcmp(name,"DamageEvolutionAnalysis")==0) return DamageEvolutionAnalysisEnum;
    330               else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    331329              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    332330              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
    … …  
    383381              else if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum;
    384382              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
     383              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
     384              else if (strcmp(name,"Contour")==0) return ContourEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    389               else if (strcmp(name,"Contour")==0) return ContourEnum;
    390               else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
     388              if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
    391389              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    392390              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    … …  
    506504              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    507505              else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
     506              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
     507              else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
    512               else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
    513               else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
     511              if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
    514512              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
    515513              else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
    … …  
    629627              else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
    630628              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
     629              else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
     630              else if (strcmp(name,"XY")==0) return XYEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum;
    635               else if (strcmp(name,"XY")==0) return XYEnum;
    636               else if (strcmp(name,"XYZ")==0) return XYZEnum;
     634              if (strcmp(name,"XYZ")==0) return XYZEnum;
    637635              else if (strcmp(name,"Dense")==0) return DenseEnum;
    638636              else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    … …  
    649647              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    650648              else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum;
     649              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
     650              else if (strcmp(name,"TransientIslevelset")==0) return TransientIslevelsetEnum;
    651651              else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
    652652         else stage=7;

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r17086 r17098  
    271271def TransientIsgiaEnum(): return StringToEnum("TransientIsgia")[0]
    272272def TransientIsdamageEnum(): return StringToEnum("TransientIsdamage")[0]
    273 def TransientIslevelsetEnum(): return StringToEnum("TransientIslevelset")[0]
    274273def TransientNumRequestedOutputsEnum(): return StringToEnum("TransientNumRequestedOutputs")[0]
    275274def TransientRequestedOutputsEnum(): return StringToEnum("TransientRequestedOutputs")[0]
    … …  
    314313def DamageEvolutionSolutionEnum(): return StringToEnum("DamageEvolutionSolution")[0]
    315314def DamageEvolutionAnalysisEnum(): return StringToEnum("DamageEvolutionAnalysis")[0]
    316 def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0]
    317315def StressbalanceAnalysisEnum(): return StringToEnum("StressbalanceAnalysis")[0]
    318316def StressbalanceSIAAnalysisEnum(): return StringToEnum("StressbalanceSIAAnalysis")[0]
    … …  
    626624def ArrheniusEnum(): return StringToEnum("Arrhenius")[0]
    627625def LliboutryDuvalEnum(): return StringToEnum("LliboutryDuval")[0]
     626def LevelsetAnalysisEnum(): return StringToEnum("LevelsetAnalysis")[0]
     627def TransientIslevelsetEnum(): return StringToEnum("TransientIslevelset")[0]
    628628def MaximumNumberOfDefinitionsEnum(): return StringToEnum("MaximumNumberOfDefinitions")[0]
Note: See TracChangeset for help on using the changeset viewer.