source: issm/oecreview/Archive/16554-17801/ISSM-16839-16840.diff@ 17802

Last change on this file since 17802 was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 18.4 KB

  • ../trunk-jpl/src/c/analyses/ThermalAnalysis.h

     
    2121
    2222                /*Finite element Analysis*/
    2323                ElementMatrix* CreateKMatrix(Element* element);
     24                ElementMatrix* CreateKMatrixVolume(Element* element);
     25                ElementMatrix* CreateKMatrixShelf(Element* element);
    2426                ElementVector* CreatePVector(Element* element);
    2527                ElementVector* CreatePVectorVolume(Element* element);
    2628                ElementVector* CreatePVectorSheet(Element* element);
    2729                ElementVector* CreatePVectorShelf(Element* element);
     30                void GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     31                void GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     32                void GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    2833                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    2934                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3035};

  • ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp

     
    112112
    113113/*Finite Element Analysis*/
    114114ElementMatrix* ThermalAnalysis::CreateKMatrix(Element* element){/*{{{*/
    115         _error_("not implemented yet");
     115
     116        /*compute all stiffness matrices for this element*/
     117        ElementMatrix* Ke1=CreateKMatrixVolume(element);
     118        ElementMatrix* Ke2=CreateKMatrixShelf(element);
     119        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     120
     121        /*clean-up and return*/
     122        delete Ke1;
     123        delete Ke2;
     124        return Ke;
    116125}/*}}}*/
     126ElementMatrix* ThermalAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
     127
     128        /*Intermediaries */
     129        int         stabilization;
     130        IssmDouble  Jdet,dt,u,v,w,um,vm,wm,vel;
     131        IssmDouble  h,hx,hy,hz,vx,vy,vz;
     132        IssmDouble  tau_parameter,diameter;
     133        IssmDouble  D_scalar;
     134        IssmDouble* xyz_list = NULL;
     135
     136        /*Fetch number of nodes and dof for this finite element*/
     137        int numnodes = element->GetNumberOfNodes();
     138
     139        /*Initialize Element vector and other vectors*/
     140        ElementMatrix* Ke     = element->NewElementMatrix();
     141        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
     142        IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
     143        IssmDouble*    B      = xNew<IssmDouble>(3*numnodes);
     144        IssmDouble*    Bprime = xNew<IssmDouble>(3*numnodes);
     145        IssmDouble     D[3][3]={0.};
     146        IssmDouble     K[3][3];
     147
     148        /*Retrieve all inputs and parameters*/
     149        element->GetVerticesCoordinates(&xyz_list);
     150        element->FindParam(&dt,TimesteppingTimeStepEnum);
     151        element->FindParam(&stabilization,ThermalStabilizationEnum);
     152        IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     153        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
     154        IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
     155        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     156        IssmDouble  thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum);
     157        IssmDouble  kappa = thermalconductivity/(rho_ice*heatcapacity);
     158        Input* vx_input  = element->GetInput(VxEnum);     _assert_(vx_input);
     159        Input* vy_input  = element->GetInput(VyEnum);     _assert_(vy_input);
     160        Input* vz_input  = element->GetInput(VzEnum);     _assert_(vz_input);
     161        Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input);
     162        Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input);
     163        Input* vzm_input = element->GetInput(VzMeshEnum); _assert_(vzm_input);
     164        if(stabilization==2) diameter=element->MinEdgeLength(xyz_list);
     165
     166        /* Start  looping on the number of gaussian points: */
     167        Gauss* gauss=element->NewGauss(2);
     168        for(int ig=gauss->begin();ig<gauss->end();ig++){
     169                gauss->GaussPoint(ig);
     170
     171                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     172                D_scalar=gauss->weight*Jdet*kappa;
     173                if(dt!=0.) D_scalar=D_scalar*dt;
     174
     175                /*Conduction: */
     176                GetBConduct(B,element,xyz_list,gauss);
     177                D[0][0]=D_scalar;
     178                D[1][1]=D_scalar;
     179                D[2][2]=D_scalar;
     180                TripleMultiply(B,3,numnodes,1,
     181                                        &D[0][0],3,3,0,
     182                                        B,3,numnodes,0,
     183                                        &Ke->values[0],1);
     184
     185                /*Advection: */
     186                GetBAdvec(B,element,xyz_list,gauss);
     187                GetBAdvecprime(Bprime,element,xyz_list,gauss);
     188                vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
     189                vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
     190                vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
     191                D[0][0]=D_scalar*vx;
     192                D[1][1]=D_scalar*vy;
     193                D[2][2]=D_scalar*vz;
     194                TripleMultiply(B,3,numnodes,1,
     195                                        &D[0][0],3,3,0,
     196                                        Bprime,3,numnodes,0,
     197                                        &Ke->values[0],1);
     198
     199                /*Transient: */
     200                if(dt!=0.){
     201                        element->NodalFunctions(basis,gauss);
     202
     203                        TripleMultiply(basis,numnodes,1,0,
     204                                                &D_scalar,1,1,0,
     205                                                basis,1,numnodes,0,
     206                                                &Ke->values[0],1);
     207                }
     208
     209                /*Artifficial diffusivity*/
     210                if(stabilization==1){
     211                        element->ElementSizes(&hx,&hy,&hz);
     212                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
     213                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
     214                        K[0][0]=h/(2.*vel)*fabs(vx*vx);  K[0][1]=h/(2.*vel)*fabs(vx*vy); K[0][2]=h/(2.*vel)*fabs(vx*vz);
     215                        K[1][0]=h/(2.*vel)*fabs(vy*vx);  K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
     216                        K[2][0]=h/(2.*vel)*fabs(vz*vx);  K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
     217                        for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
     218
     219                        GetBAdvecprime(Bprime,element,xyz_list,gauss);
     220
     221                        TripleMultiply(Bprime,3,numnodes,1,
     222                                                &K[0][0],3,3,0,
     223                                                Bprime,3,numnodes,0,
     224                                                &Ke->values[0],1);
     225                }
     226                else if(stabilization==2){
     227                        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     228                        tau_parameter=element->StabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
     229                        for(int i=0;i<numnodes;i++){
     230                                for(int j=0;j<numnodes;j++){
     231                                        Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*
     232                                          ((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i])*((u-um)*dbasis[0*3+j]+(v-vm)*dbasis[1*3+j]+(w-wm)*dbasis[2*3+j]);
     233                                }
     234                        }
     235                        if(dt!=0.){
     236                                for(int i=0;i<numnodes;i++){
     237                                        for(int j=0;j<numnodes;j++){
     238                                                Ke->values[i*numnodes+j]+=tau_parameter*D_scalar*basis[j]*((u-um)*dbasis[0*3+i]+(v-vm)*dbasis[1*3+i]+(w-wm)*dbasis[2*3+i]);
     239                                        }
     240                                }
     241                        }
     242                }
     243        }
     244
     245        /*Clean up and return*/
     246        delete gauss;
     247        return Ke;
     248}/*}}}*/
     249ElementMatrix* ThermalAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/
     250
     251        /*Initialize Element matrix and return if necessary*/
     252        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     253
     254        IssmDouble  dt,Jdet,D;
     255        IssmDouble *xyz_list_base = NULL;
     256
     257        /*Get basal element*/
     258        if(!element->IsOnBed() || !element->IsFloating()) return NULL;
     259
     260        /*Fetch number of nodes for this finite element*/
     261        int numnodes = element->GetNumberOfNodes();
     262
     263        /*Initialize vectors*/
     264        ElementMatrix* Ke    = element->NewElementMatrix();
     265        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     266
     267        /*Retrieve all inputs and parameters*/
     268        element->GetVerticesCoordinatesBase(&xyz_list_base);
     269        element->FindParam(&dt,TimesteppingTimeStepEnum);
     270        IssmDouble  gravity             = element->GetMaterialParameter(ConstantsGEnum);
     271        IssmDouble  rho_water           = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     272        IssmDouble  rho_ice             = element->GetMaterialParameter(MaterialsRhoIceEnum);
     273        IssmDouble  heatcapacity        = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     274        IssmDouble  mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
     275        IssmDouble  thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
     276
     277        /* Start  looping on the number of gaussian points: */
     278        Gauss* gauss=element->NewGaussBase(2);
     279        for(int ig=gauss->begin();ig<gauss->end();ig++){
     280                gauss->GaussPoint(ig);
     281
     282                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     283                element->NodalFunctions(basis,gauss);
     284
     285                D=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel/(heatcapacity*rho_ice);
     286                if(reCast<bool,IssmDouble>(dt)) D=dt*D;
     287                TripleMultiply(basis,numnodes,1,0,
     288                                        &D,1,1,0,
     289                                        basis,1,numnodes,0,
     290                                        &Ke->values[0],1);
     291
     292        }
     293
     294        /*Clean up and return*/
     295        delete gauss;
     296        xDelete<IssmDouble>(basis);
     297        xDelete<IssmDouble>(xyz_list_base);
     298        return Ke;
     299}/*}}}*/
    117300ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/
    118301
    119302        /*compute all load vectors for this element*/
    … …  
    262445void ThermalAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    263446        element->GetSolutionFromInputsOneDof(solution,TemperatureEnum);
    264447}/*}}}*/
     448void ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     449        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
     450         * For node i, Bi' can be expressed in the actual coordinate system
     451         * by:
     452         *       Bi_conduct=[ dh/dx ]
     453         *                  [ dh/dy ]
     454         *                  [ dh/dz ]
     455         * where h is the interpolation function for node i.
     456         *
     457         * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
     458         */
     459
     460        /*Fetch number of nodes for this finite element*/
     461        int numnodes = element->GetNumberOfNodes();
     462
     463        /*Get nodal functions derivatives*/
     464        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     465        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     466
     467        /*Build B: */
     468        for(int i=0;i<numnodes;i++){
     469                B[numnodes*0+i] = dbasis[0*numnodes+i];
     470                B[numnodes*1+i] = dbasis[1*numnodes+i];
     471                B[numnodes*2+i] = dbasis[2*numnodes+i];
     472        }
     473
     474        /*Clean-up*/
     475        xDelete<IssmDouble>(dbasis);
     476}/*}}}*/
     477void ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     478        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
     479         * For node i, Bi' can be expressed in the actual coordinate system
     480         * by:
     481         *       Bi_advec =[ h ]
     482         *                 [ h ]
     483         *                 [ h ]
     484         * where h is the interpolation function for node i.
     485         *
     486         * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
     487         */
     488
     489        /*Fetch number of nodes for this finite element*/
     490        int numnodes = element->GetNumberOfNodes();
     491
     492        /*Get nodal functions*/
     493        IssmDouble* basis=xNew<IssmDouble>(numnodes);
     494        element->NodalFunctions(basis,gauss);
     495
     496        /*Build B: */
     497        for(int i=0;i<numnodes;i++){
     498                B[numnodes*0+i] = basis[i];
     499                B[numnodes*1+i] = basis[i];
     500                B[numnodes*2+i] = basis[i];
     501        }
     502
     503        /*Clean-up*/
     504        xDelete<IssmDouble>(basis);
     505}/*}}}*/
     506void ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     507        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
     508         * For node i, Bi' can be expressed in the actual coordinate system
     509         * by:
     510         *       Biprime_advec=[ dh/dx ]
     511         *                     [ dh/dy ]
     512         *                     [ dh/dz ]
     513         * where h is the interpolation function for node i.
     514         *
     515         * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
     516         */
     517
     518        /*Fetch number of nodes for this finite element*/
     519        int numnodes = element->GetNumberOfNodes();
     520
     521        /*Get nodal functions derivatives*/
     522        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
     523        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     524
     525        /*Build B: */
     526        for(int i=0;i<numnodes;i++){
     527                B[numnodes*0+i] = dbasis[0*numnodes+i];
     528                B[numnodes*1+i] = dbasis[1*numnodes+i];
     529                B[numnodes*2+i] = dbasis[2*numnodes+i];
     530        }
     531
     532        /*Clean-up*/
     533        xDelete<IssmDouble>(dbasis);
     534}/*}}}*/
    265535void ThermalAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    266536
    267537        bool        converged;

  • ../trunk-jpl/src/c/analyses/MeltingAnalysis.cpp

     
    7272
    7373/*Finite Element Analysis*/
    7474ElementMatrix* MeltingAnalysis::CreateKMatrix(Element* element){/*{{{*/
    75         _error_("not implemented yet");
     75
     76        /*Get basal element*/
     77        if(!element->IsOnBed()) return NULL;
     78        Element* basalelement = element->SpawnBasalElement();
     79
     80        /*Intermediaries */
     81        IssmDouble  D,Jdet;
     82        IssmDouble *xyz_list  = NULL;
     83
     84        /*Fetch number of nodes and dof for this finite element*/
     85        int numnodes = basalelement->GetNumberOfNodes();
     86
     87        /*Initialize Element vector*/
     88        ElementMatrix* Ke    = basalelement->NewElementMatrix(NoneApproximationEnum);
     89        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     90
     91        /*Retrieve all inputs and parameters*/
     92        basalelement->GetVerticesCoordinates(&xyz_list);
     93        IssmDouble latentheat   = element->GetMaterialParameter(MaterialsLatentheatEnum);
     94        IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum);
     95
     96        /* Start  looping on the number of gaussian points: */
     97        Gauss* gauss=basalelement->NewGauss(2);
     98        for(int ig=gauss->begin();ig<gauss->end();ig++){
     99                gauss->GaussPoint(ig);
     100
     101                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     102                basalelement->NodalFunctions(basis,gauss);
     103                D=latentheat/heatcapacity*gauss->weight*Jdet;
     104
     105                TripleMultiply(basis,1,numnodes,1,
     106                                        &D,1,1,0,
     107                                        basis,1,numnodes,0,
     108                                        &Ke->values[0],1);
     109        }
     110
     111        /*Clean up and return*/
     112        xDelete<IssmDouble>(xyz_list);
     113        xDelete<IssmDouble>(basis);
     114        delete gauss;
     115        basalelement->DeleteMaterials(); delete basalelement;
     116        return Ke;
    76117}/*}}}*/
    77118ElementVector* MeltingAnalysis::CreatePVector(Element* element){/*{{{*/
    78119        return NULL;

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

     
    4848                virtual ElementVector* CreatePVector(void)=0;
    4949                virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
    5050                virtual void   DeleteMaterials(void)=0;
     51                virtual void   ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    5152                virtual void   EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0;
    5253                virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0;
    5354                virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0;

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

     
    7878                void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
    7979                void        DeleteMaterials(void);
    8080                void        Delta18oParameterization(void);
     81                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    8182                void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
    8283                void        FindParam(int* pvalue,int paramenum);
    8384                void        FindParam(IssmDouble* pvalue,int paramenum);

  • ../trunk-jpl/src/c/classes/Elements/Penta.cpp

     
    12091209        return this->element_type;
    12101210}
    12111211/*}}}*/
    1212 /*FUNCTION Penta::GetElementSizes{{{*/
    1213 void Penta::GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){
     1212/*FUNCTION Penta::ElementSizes{{{*/
     1213void Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){
    12141214
    12151215        IssmDouble xyz_list[NUMVERTICES][3];
    12161216        IssmDouble xmin,ymin,zmin;
    … …  
    39803980                /*Artificial diffusivity*/
    39813981                if(stabilization==1){
    39823982                        /*Build K: */
    3983                         GetElementSizes(&hx,&hy,&hz);
     3983                        ElementSizes(&hx,&hy,&hz);
    39843984                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
    39853985                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
    39863986
    … …  
    42094209                /*Artifficial diffusivity*/
    42104210                if(stabilization==1){
    42114211                        /*Build K: */
    4212                         GetElementSizes(&hx,&hy,&hz);
     4212                        ElementSizes(&hx,&hy,&hz);
    42134213                        vel=sqrt(vx*vx + vy*vy + vz*vz)+1.e-14;
    42144214                        h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) + pow(hz*vz/vel,2));
    42154215

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

     
    7474                void   ComputeStressTensor();
    7575                void   Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
    7676                void   DeleteMaterials(void){_error_("not implemented yet");};
     77                void   ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    7778                void   FindParam(int* pvalue,int paramenum);
    7879                void   FindParam(IssmDouble* pvalue,int paramenum);
    7980                int    FiniteElement(void);
    … …  
    220221                void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
    221222                IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
    222223                int            GetElementType(void);
    223                 void           GetElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    224224                Input*         GetInput(int inputenum);
    225225                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    226226                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);

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

     
    8383                ElementVector* CreatePVectorL2Projection(void);
    8484                void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("not implemented yet");};
    8585                void        Delta18oParameterization(void){_error_("not implemented yet");};
     86                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    8687                void        EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};
    8788                IssmDouble  EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
    8889                IssmDouble  EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
Note: See TracBrowser for help on using the repository browser.