source: issm/oecreview/Archive/16554-17801/ISSM-17390-17391.diff@ 17802

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

Added archives
File size: 15.2 KB

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

     
    2525                ElementMatrix* CreateJacobianMatrix(Element* element);
    2626                ElementMatrix* CreateKMatrix(Element* element);
    2727                ElementVector* CreatePVector(Element* element);
    28                 void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
    29                 void GetBprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
     28                void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     29                void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    3030                void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3131                void InputUpdateFromSolution(IssmDouble* solution,Element* element);
    3232                void UpdateConstraints(FemModel* femmodel);

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

     
    116116}/*}}}*/
    117117ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/
    118118
     119        /* Check if ice in element */
     120        if(!element->IsIceInElement()) return NULL;
     121
    119122        /*Intermediaries*/
    120         int      meshtype;
    121         Element* basalelement;
     123        Element*    basalelement;
     124        int         meshtype,dim;
     125        int         stabilization;
     126        IssmDouble  Jdet,dt,D_scalar,h;
     127        IssmDouble  vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
     128        IssmDouble *xyz_list  = NULL;
    122129
    123         /*Get basal element*/
     130        /*Get problem dimension and basal element*/
    124131        element->FindParam(&meshtype,MeshTypeEnum);
    125132        switch(meshtype){
    126133                case Mesh2DhorizontalEnum:
    127134                        basalelement = element;
     135                        dim = 2;
    128136                        break;
    129137                case Mesh3DEnum:
    130138                        if(!element->IsOnBed()) return NULL;
    131139                        basalelement = element->SpawnBasalElement();
     140                        dim = 2;
    132141                        break;
    133142                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    134143        }
    135144
    136         /*Intermediaries */
    137         int         stabilization;
    138         IssmDouble  Jdet,dt,D_scalar;
    139         IssmDouble  vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
    140         IssmDouble *xyz_list  = NULL;
    141 
    142145        /*Fetch number of nodes and dof for this finite element*/
    143146        int numnodes = basalelement->GetNumberOfNodes();
    144147
    145148        /*Initialize Element vector*/
    146         ElementMatrix* Ke     = basalelement->NewElementMatrix(NoneApproximationEnum);
     149        ElementMatrix* Ke     = basalelement->NewElementMatrix();
    147150        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    148         IssmDouble*    B      = xNew<IssmDouble>(2*numnodes);
    149         IssmDouble*    Bprime = xNew<IssmDouble>(2*numnodes);
    150         IssmDouble     D[2][2]={0.};
     151        IssmDouble*    B      = xNew<IssmDouble>(dim*numnodes);
     152        IssmDouble*    Bprime = xNew<IssmDouble>(dim*numnodes);
     153        IssmDouble*    D      = xNewZeroInit<IssmDouble>(dim*dim);
    151154
    152155        /*Retrieve all inputs and parameters*/
    153156        basalelement->GetVerticesCoordinates(&xyz_list);
    … …  
    160163                vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input);
    161164        }
    162165        else{
    163                 vxaverage_input=basalelement->GetInput(VxAverageEnum); _assert_(vxaverage_input);
    164                 vyaverage_input=basalelement->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     166                if(dim==1){
     167                        vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
     168                }
     169                if(dim==2){
     170                        vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
     171                        vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
     172                }
    165173        }
    166         IssmDouble h=basalelement->CharacteristicLength();
     174        h=basalelement->CharacteristicLength();
    167175
    168176        /* Start  looping on the number of gaussian points: */
    169177        Gauss* gauss=basalelement->NewGauss(2);
    … …  
    172180
    173181                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    174182                basalelement->NodalFunctions(basis,gauss);
    175                 D_scalar=gauss->weight*Jdet;
    176 
     183               
    177184                vxaverage_input->GetInputValue(&vx,gauss);
    178                 vyaverage_input->GetInputValue(&vy,gauss);
    179185                vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
    180                 vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     186                if(dim==2){
     187                        vyaverage_input->GetInputValue(&vy,gauss);
     188                        vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     189                }
    181190
     191                D_scalar=gauss->weight*Jdet;
    182192                TripleMultiply(basis,1,numnodes,1,
    183193                                        &D_scalar,1,1,0,
    184194                                        basis,1,numnodes,0,
    185195                                        &Ke->values[0],1);
    186196
    187                 GetB(B,element,xyz_list,gauss);
    188                 GetBprime(Bprime,element,xyz_list,gauss);
     197                GetB(B,basalelement,dim,xyz_list,gauss);
     198                GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
    189199
    190200                dvxdx=dvx[0];
    191                 dvydy=dvy[1];
     201                if(dim==2) dvydy=dvy[1];
    192202                D_scalar=dt*gauss->weight*Jdet;
    193203
    194                 D[0][0]=D_scalar*dvxdx;
    195                 D[1][1]=D_scalar*dvydy;
    196                 TripleMultiply(B,2,numnodes,1,
    197                                         &D[0][0],2,2,0,
    198                                         B,2,numnodes,0,
     204                D[0*dim+0]=D_scalar*dvxdx;
     205                if(dim==2) D[1*dim+1]=D_scalar*dvydy;
     206
     207                TripleMultiply(B,dim,numnodes,1,
     208                                        D,dim,dim,0,
     209                                        B,dim,numnodes,0,
    199210                                        &Ke->values[0],1);
    200211
    201                 D[0][0]=D_scalar*vx;
    202                 D[1][1]=D_scalar*vy;
    203                 TripleMultiply(B,2,numnodes,1,
    204                                         &D[0][0],2,2,0,
    205                                         Bprime,2,numnodes,0,
     212                D[0*dim+0]=D_scalar*vx;
     213                if(dim==2) D[1*dim+1]=D_scalar*vy;
     214
     215                TripleMultiply(B,dim,numnodes,1,
     216                                        D,dim,dim,0,
     217                                        Bprime,dim,numnodes,0,
    206218                                        &Ke->values[0],1);
    207219
    208220                if(stabilization==2){
    209                         /*Streamline upwinding*/
    210                         vel=sqrt(vx*vx+vy*vy)+1.e-8;
    211                         D[0][0]=h/(2*vel)*vx*vx;
    212                         D[1][0]=h/(2*vel)*vy*vx;
    213                         D[0][1]=h/(2*vel)*vx*vy;
    214                         D[1][1]=h/(2*vel)*vy*vy;
     221                        if(dim==1){
     222                                vel=fabs(vx)+1.e-8;
     223                                D[0]=h/(2*vel)*vx*vx;
     224                        }
     225                        else{
     226                                /*Streamline upwinding*/
     227                                vel=sqrt(vx*vx+vy*vy)+1.e-8;
     228                                D[0*dim+0]=h/(2*vel)*vx*vx;
     229                                D[1*dim+0]=h/(2*vel)*vy*vx;
     230                                D[0*dim+1]=h/(2*vel)*vx*vy;
     231                                D[1*dim+1]=h/(2*vel)*vy*vy;
     232                        }
    215233                }
    216234                else if(stabilization==1){
    217                         /*SSA*/
    218235                        vxaverage_input->GetInputAverage(&vx);
    219                         vyaverage_input->GetInputAverage(&vy);
    220                         D[0][0]=h/2.0*fabs(vx);
    221                         D[0][1]=0.;
    222                         D[1][0]=0.;
    223                         D[1][1]=h/2.0*fabs(vy);
     236                        if(dim==2) vyaverage_input->GetInputAverage(&vy);
     237                        D[0*dim+0]=h/2.0*fabs(vx);
     238                        if(dim==2) D[1*dim+1]=h/2.0*fabs(vy);
    224239                }
    225240                if(stabilization==1 || stabilization==2){
    226                         D[0][0]=D_scalar*D[0][0];
    227                         D[1][0]=D_scalar*D[1][0];
    228                         D[0][1]=D_scalar*D[0][1];
    229                         D[1][1]=D_scalar*D[1][1];
    230                         TripleMultiply(Bprime,2,numnodes,1,
    231                                                 &D[0][0],2,2,0,
    232                                                 Bprime,2,numnodes,0,
     241                        if(dim==1) D[0]=D_scalar*D[0];
     242                        else{
     243                                D[0*dim+0]=D_scalar*D[0*dim+0];
     244                                D[1*dim+0]=D_scalar*D[1*dim+0];
     245                                D[0*dim+1]=D_scalar*D[0*dim+1];
     246                                D[1*dim+1]=D_scalar*D[1*dim+1];
     247                        }
     248
     249                        TripleMultiply(Bprime,dim,numnodes,1,
     250                                                D,dim,dim,0,
     251                                                Bprime,dim,numnodes,0,
    233252                                                &Ke->values[0],1);
    234253                }
    235254
    … …  
    240259        xDelete<IssmDouble>(basis);
    241260        xDelete<IssmDouble>(B);
    242261        xDelete<IssmDouble>(Bprime);
     262        xDelete<IssmDouble>(D);
    243263        delete gauss;
    244264        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    245265        return Ke;
    246266}/*}}}*/
    247267ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
    248268
     269        /* Check if ice in element */
     270        if(!element->IsIceInElement()) return NULL;
     271
    249272        /*Intermediaries*/
    250273        int      meshtype;
    251274        Element* basalelement;
     275        IssmDouble  Jdet,dt;
     276        IssmDouble  f,damage;
     277        IssmDouble* xyz_list = NULL;
    252278
    253279        /*Get basal element*/
    254280        element->FindParam(&meshtype,MeshTypeEnum);
    … …  
    263289                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    264290        }
    265291
    266         /*Intermediaries */
    267         IssmDouble  Jdet,dt;
    268         IssmDouble  f,damage;
    269         IssmDouble* xyz_list = NULL;
    270 
    271292        /*Fetch number of nodes and dof for this finite element*/
    272         int numnodes = element->GetNumberOfNodes();
     293        int numnodes = basalelement->GetNumberOfNodes();
    273294
    274295        /*Initialize Element vector and other vectors*/
    275         ElementVector* pe    = element->NewElementVector();
     296        ElementVector* pe    = basalelement->NewElementVector();
    276297        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    277298
    278299        /*Retrieve all inputs and parameters*/
    279         element->GetVerticesCoordinates(&xyz_list);
    280         element->FindParam(&dt,TimesteppingTimeStepEnum);
    281         this->CreateDamageFInput(element);
    282         Input* damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input);
    283         Input* damagef_input = element->GetInput(DamageFEnum);    _assert_(damagef_input);
     300        basalelement->GetVerticesCoordinates(&xyz_list);
     301        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
     302        this->CreateDamageFInput(basalelement);
     303        Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input);
     304        Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input);
    284305
     306
    285307        /* Start  looping on the number of gaussian points: */
    286         Gauss* gauss=element->NewGauss(2);
     308        Gauss* gauss=basalelement->NewGauss(2);
    287309        for(int ig=gauss->begin();ig<gauss->end();ig++){
    288310                gauss->GaussPoint(ig);
    289311
    290                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    291                 element->NodalFunctions(basis,gauss);
     312                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     313                basalelement->NodalFunctions(basis,gauss);
    292314
    293315                damaged_input->GetInputValue(&damage,gauss);
    294316                damagef_input->GetInputValue(&f,gauss);
    295317
    296                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
     318                for(int i=0;i<numnodes;i++){
     319                        pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
     320                }
    297321        }
    298322
    299323        /*Clean up and return*/
    … …  
    303327        delete gauss;
    304328        return pe;
    305329}/*}}}*/
    306 void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     330void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    307331        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    308332         * For node i, Bi can be expressed in the actual coordinate system
    309333         * by:
    … …  
    323347
    324348        /*Build B: */
    325349        for(int i=0;i<numnodes;i++){
    326                 B[numnodes*0+i] = basis[i];
    327                 B[numnodes*1+i] = basis[i];
     350                for(int j=0;j<dim;j++){
     351                        B[numnodes*j+i] = basis[i];
     352                }
    328353        }
    329354
    330355        /*Clean-up*/
    331356        xDelete<IssmDouble>(basis);
    332357}/*}}}*/
    333 void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     358void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    334359        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    335360         * For node i, Bi' can be expressed in the actual coordinate system
    336361         * by:
    … …  
    345370        int numnodes = element->GetNumberOfNodes();
    346371
    347372        /*Get nodal functions derivatives*/
    348         IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
     373        IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    349374        element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    350375
    351376        /*Build B': */
    352377        for(int i=0;i<numnodes;i++){
    353                 Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
    354                 Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
     378                for(int j=0;j<dim;j++){
     379                        Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
     380                }
    355381        }
    356382
    357383        /*Clean-up*/
    … …  
    363389}/*}}}*/
    364390void DamageEvolutionAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    365391
     392        int meshtype;
    366393        IssmDouble  max_damage;
    367394        int                     *doflist = NULL;
     395        Element*   basalelement=NULL;
    368396
     397        element->FindParam(&meshtype,MeshTypeEnum);
     398        if(meshtype!=Mesh2DhorizontalEnum){
     399                if(!element->IsOnBed()) return;
     400                basalelement=element->SpawnBasalElement();
     401        }
     402        else{
     403                basalelement = element;
     404        }
    369405        /*Fetch number of nodes and dof for this finite element*/
    370         int numnodes = element->GetNumberOfNodes();
     406        int numnodes = basalelement->GetNumberOfNodes();
    371407
    372408        /*Fetch dof list and allocate solution vector*/
    373         element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    374         IssmDouble* values    = xNew<IssmDouble>(numnodes);
     409        basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     410        IssmDouble* newdamage = xNew<IssmDouble>(numnodes);
    375411
    376412        /*Get user-supplied max_damage: */
    377         element->FindParam(&max_damage,DamageMaxDamageEnum);
     413        basalelement->FindParam(&max_damage,DamageMaxDamageEnum);
    378414
    379415        /*Use the dof list to index into the solution vector: */
    380416        for(int i=0;i<numnodes;i++){
    381                 values[i]=solution[doflist[i]];
     417                newdamage[i]=solution[doflist[i]];
    382418                /*Check solution*/
    383                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
     419                if(xIsNan<IssmDouble>(newdamage[i])) _error_("NaN found in solution vector");
    384420                /*Enforce D < max_damage and D > 0 */
    385                 if(values[i]>max_damage) values[i]=max_damage;
    386                 else if(values[i]<0.)    values[i]=0.;
     421                if(newdamage[i]>max_damage) newdamage[i]=max_damage;
     422                else if(newdamage[i]<0.)    newdamage[i]=0.;
    387423        }
    388424
    389425        /*Get all inputs and parameters*/
    390         element->AddInput(DamageDbarEnum,values,P1Enum);
     426        basalelement->AddBasalInput(DamageDEnum,newdamage,P1Enum);
    391427
    392428        /*Free ressources:*/
    393         xDelete<IssmDouble>(values);
     429        xDelete<IssmDouble>(newdamage);
    394430        xDelete<int>(doflist);
     431        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    395432}/*}}}*/
    396433void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    397434        /*Default, do nothing*/
    … …  
    426463        Input* sigma_xx_input  = element->GetInput(StressTensorxxEnum);     _assert_(sigma_xx_input);
    427464        Input* sigma_xy_input  = element->GetInput(StressTensorxyEnum);     _assert_(sigma_xy_input);
    428465        Input* sigma_yy_input  = element->GetInput(StressTensoryyEnum);     _assert_(sigma_yy_input);
    429         Input* damage_input    = element->GetInput(DamageDbarEnum); _assert_(damage_input);
     466        Input* damage_input    = element->GetInput(DamageDEnum); _assert_(damage_input);
    430467
    431468        /*retrieve the desired type of equivalent stress*/
    432469        element->FindParam(&equivstress,DamageEquivStressEnum);
    … …  
    445482                s_xy=sigma_xy/(1.-damage);
    446483                s_yy=sigma_yy/(1.-damage);
    447484
    448                 if(equivstress==1){ /* von Mises */
    449                         J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy);
    450                         Chi=sqrt(3.0)*J2s;
     485                if(equivstress==0){ /* von Mises */
     486                        Chi=sqrt(s_xx*s_xx - s_xx*s_yy + s_yy*s_yy + 3*s_xy*s_xy);
    451487                }
    452488                Psi=Chi-stress_threshold;
    453489                PosPsi=max(Psi,0.);
    454                 NegPsi=max(-Psi,0.);
     490                NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
    455491
    456492                f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3);
    457493        }

  • ../trunk-jpl/src/c/cores/damage_core.cpp

     
    1818        int    numoutputs          = 0;
    1919        char   **requested_outputs = NULL;
    2020
    21         if(VerboseSolution()) _printf0_("   computing damage\n");
    2221       
    2322        //first recover parameters common to all solutions
    2423        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    … …  
    3231                ResetConstraintsx(femmodel);
    3332        }
    3433
     34        if(VerboseSolution()) _printf0_("   computing damage\n");
    3535        femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
    3636        solutionsequence_linear(femmodel);
    3737

  • ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp

     
    712712        }
    713713
    714714        //First recover damage  using the element: */
    715         element->GetInputValue(&damage,node,DamageDbarEnum);
     715        element->GetInputValue(&damage,node,DamageDEnum);
    716716
    717717        //Recover our data:
    718718        parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);
Note: See TracBrowser for help on using the repository browser.