Changeset 11332


Ignore:
Timestamp:
02/06/12 10:55:48 (13 years ago)
Author:
Mathieu Morlighem
Message:

Added Newton's method for MacAyeal

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

    r11322 r11332  
    62576257}
    62586258/*}}}*/
    6259 /*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
    6260 ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
    6261 
    6262         /*Constants*/
    6263         const int    numdof=NDOF2*NUMVERTICES;
    6264 
    6265         /*Intermediaries */
    6266         int        i,j,ig;
    6267         int        approximation;
    6268         double     xyz_list[NUMVERTICES][3];
    6269         double     Jdet;
    6270         double     eps1dotdphii,eps1dotdphij;
    6271         double     eps2dotdphii,eps2dotdphij;
    6272         double     mu_prime;
    6273         double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    6274         double     eps1[3],eps2[3];
    6275         double     phi[NUMVERTICES];
    6276         double     dphi[3][NUMVERTICES];
    6277         GaussPenta *gauss=NULL;
    6278 
    6279         /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
    6280         ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
    6281 
    6282         /*Retrieve all inputs and parameters*/
    6283         inputs->GetInputValue(&approximation,ApproximationEnum);
    6284         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    6285         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    6286         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    6287 
    6288         /* Start  looping on the number of gaussian points: */
    6289         gauss=new GaussPenta(5,5);
    6290         for (ig=gauss->begin();ig<gauss->end();ig++){
    6291 
    6292                 gauss->GaussPoint(ig);
    6293 
    6294                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6295                 GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
    6296 
    6297                 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    6298                 matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    6299                 eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    6300                 eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    6301                 eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
    6302 
    6303                 for(i=0;i<6;i++){
    6304                         for(j=0;j<6;j++){
    6305                                 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
    6306                                 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
    6307                                 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
    6308                                 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
    6309 
    6310                                 Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
    6311                                 Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
    6312                                 Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
    6313                                 Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
    6314                         }
    6315                 }
    6316         }
    6317 
    6318         /*Transform Coordinate System*/
    6319         TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
    6320 
    6321         /*Clean up and return*/
    6322         delete gauss;
    6323         return Ke;
    6324 }
    6325 /*}}}*/
    63266259/*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{1*/
    63276260ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){
    … …  
    74917424}
    74927425/*}}}*/
     7426/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
     7427ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
     7428
     7429        /*Constants*/
     7430        const int    numdof=NDOF2*NUMVERTICES;
     7431
     7432        /*Intermediaries */
     7433        int        i,j,ig;
     7434        int        approximation;
     7435        double     xyz_list[NUMVERTICES][3];
     7436        double     Jdet;
     7437        double     eps1dotdphii,eps1dotdphij;
     7438        double     eps2dotdphii,eps2dotdphij;
     7439        double     mu_prime;
     7440        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     7441        double     eps1[3],eps2[3];
     7442        double     phi[NUMVERTICES];
     7443        double     dphi[3][NUMVERTICES];
     7444        GaussPenta *gauss=NULL;
     7445
     7446        /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
     7447        ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
     7448
     7449        /*Retrieve all inputs and parameters*/
     7450        inputs->GetInputValue(&approximation,ApproximationEnum);
     7451        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     7452        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     7453        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     7454
     7455        /* Start  looping on the number of gaussian points: */
     7456        gauss=new GaussPenta(5,5);
     7457        for (ig=gauss->begin();ig<gauss->end();ig++){
     7458
     7459                gauss->GaussPoint(ig);
     7460
     7461                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7462                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
     7463
     7464                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     7465                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     7466                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
     7467                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
     7468                eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
     7469
     7470                for(i=0;i<6;i++){
     7471                        for(j=0;j<6;j++){
     7472                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
     7473                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
     7474                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
     7475                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
     7476
     7477                                Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     7478                                Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     7479                                Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     7480                                Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     7481                        }
     7482                }
     7483        }
     7484
     7485        /*Transform Coordinate System*/
     7486        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
     7487
     7488        /*Clean up and return*/
     7489        delete gauss;
     7490        return Ke;
     7491}
     7492/*}}}*/
    74937493/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    74947494void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){

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

    r11322 r11332  
    886886        delete gauss;
    887887        return pe;
     888}
     889/*}}}*/
     890/*FUNCTION Tria::CreateJacobianMatrix{{{1*/
     891void  Tria::CreateJacobianMatrix(Mat Jff){
     892
     893        /*retrieve parameters: */
     894        ElementMatrix* Ke=NULL;
     895        int analysis_type;
     896        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     897
     898        /*Checks in debugging {{{2*/
     899        _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
     900        /*}}}*/
     901
     902        /*Skip if water element*/
     903        if(IsOnWater()) return;
     904
     905        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     906        switch(analysis_type){
     907#ifdef _HAVE_DIAGNOSTIC_
     908                case DiagnosticHorizAnalysisEnum:
     909                        Ke=CreateJacobianDiagnosticMacayeal();
     910                        break;
     911#endif
     912                default:
     913                        _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
     914        }
     915
     916        /*Add to global matrix*/
     917        if(Ke){
     918                Ke->AddToGlobal(Jff);
     919                delete Ke;
     920        }
    888921}
    889922/*}}}*/
    … …  
    30283061        delete gauss;
    30293062        return pe;
     3063}
     3064/*}}}*/
     3065/*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/
     3066ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
     3067
     3068        /*Constants*/
     3069        const int    numdof=NDOF2*NUMVERTICES;
     3070
     3071        /*Intermediaries */
     3072        int        i,j,ig;
     3073        double     xyz_list[NUMVERTICES][3];
     3074        double     Jdet,thickness;
     3075        double     eps1dotdphii,eps1dotdphij;
     3076        double     eps2dotdphii,eps2dotdphij;
     3077        double     mu_prime;
     3078        double     epsilon[3];/* epsilon=[exx,eyy,exy];*/
     3079        double     eps1[2],eps2[2];
     3080        double     phi[NUMVERTICES];
     3081        double     dphi[2][NUMVERTICES];
     3082        GaussTria *gauss=NULL;
     3083
     3084        /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
     3085        ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
     3086
     3087        /*Retrieve all inputs and parameters*/
     3088        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3089        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     3090        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     3091        Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
     3092
     3093        /* Start  looping on the number of gaussian points: */
     3094        gauss=new GaussTria(2);
     3095        for (ig=gauss->begin();ig<gauss->end();ig++){
     3096
     3097                gauss->GaussPoint(ig);
     3098
     3099                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
     3100                GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
     3101
     3102                thickness_input->GetInputValue(&thickness, gauss);
     3103                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     3104                matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     3105                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
     3106                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
     3107
     3108                for(i=0;i<3;i++){
     3109                        for(j=0;j<3;j++){
     3110                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
     3111                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
     3112                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
     3113                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
     3114
     3115                                Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
     3116                                Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
     3117                                Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
     3118                                Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
     3119                        }
     3120                }
     3121        }
     3122
     3123        /*Transform Coordinate System*/
     3124        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
     3125
     3126        /*Clean up and return*/
     3127        delete gauss;
     3128        return Ke;
    30303129}
    30313130/*}}}*/

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

    r11322 r11332  
    8484                void   CreateKMatrix(Mat Kff, Mat Kfs,Vec df);
    8585                void   CreatePVector(Vec pf);
    86                 void   CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
     86                void   CreateJacobianMatrix(Mat Jff);
    8787                int    GetNodeIndex(Node* node);
    8888                int    Sid();
    … …  
    208208                ElementVector* CreatePVectorDiagnosticMacAyeal(void);
    209209                ElementVector* CreatePVectorDiagnosticHutter(void);
     210                ElementMatrix* CreateJacobianDiagnosticMacayeal(void);
    210211                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    211212                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);

  • issm/trunk-jpl/src/c/objects/Loads/Icefront.cpp

    r10576 r11332  
    361361}
    362362/*}}}*/
     363/*FUNCTION Icefront::CreateJacobianMatrix{{{1*/
     364void  Icefront::CreateJacobianMatrix(Mat Jff){
     365        this->CreateKMatrix(Jff,NULL);
     366}
     367/*}}}1*/
    363368/*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/
    364369void  Icefront::PenaltyCreateKMatrix(Mat Kff, Mat Kfs, double kmax){
    … …  
    373378}
    374379/*}}}*/
     380/*FUNCTION Icefront::PenaltyCreateJacobianMatrix{{{1*/
     381void  Icefront::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
     382        this->PenaltyCreateKMatrix(Jff,NULL,kmax);
     383}
     384/*}}}1*/
    375385/*FUNCTION Icefront::InAnalysis{{{1*/
    376386bool Icefront::InAnalysis(int in_analysis_type){

  • issm/trunk-jpl/src/c/objects/Loads/Icefront.h

    r11327 r11332  
    7676                void  CreateKMatrix(Mat Kff, Mat Kfs);
    7777                void  CreatePVector(Vec pf);
    78                 void  CreateJacobianMatrix(Mat Jff){_error_("Not implemented yet");};
    79                 void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax){_error_("Not implemented yet");};
     78                void  CreateJacobianMatrix(Mat Jff);
    8079                void  PenaltyCreateKMatrix(Mat Kff, Mat kfs, double kmax);
    8180                void  PenaltyCreatePVector(Vec pf, double kmax);
     81                void  PenaltyCreateJacobianMatrix(Mat Jff,double kmax);
    8282                bool  InAnalysis(int analysis_type);
    8383                /*}}}*/

  • issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp

    r11327 r11332  
    216216/*FUNCTION Penpair::CreateJacobianMatrix{{{1*/
    217217void  Penpair::CreateJacobianMatrix(Mat Jff){
    218         /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/
    219         /*No loads applied, do nothing: */
    220         return;
    221 
     218        this->CreateKMatrix(Jff,NULL);
    222219}
    223220/*}}}1*/

  • issm/trunk-jpl/src/c/objects/Materials/Matice.cpp

    r11322 r11332  
    549549                eyz=epsilon[4];
    550550                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy + exz*exz + eyz*eyz;
     551
     552                mu_prime=(1-n)/(2*n) * mu/eff2;
     553        }
     554
     555        /*Assign output pointers:*/
     556        *pmu_prime=mu_prime;
     557}
     558/*}}}*/
     559/*FUNCTION Matice::GetViscosity2dDerivativeEpsSquare{{{1*/
     560void  Matice::GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* epsilon){
     561
     562        /*output: */
     563        double mu_prime;
     564        double mu,n,eff2;
     565
     566        /*input strain rate: */
     567        double exx,eyy,exy,exz;
     568
     569        /*Get visocisty and n*/
     570        GetViscosity2d(&mu,epsilon);
     571        n=GetN();
     572
     573        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
     574                mu_prime=0.5*pow((double)10,(double)14);
     575        }
     576        else{
     577                /*Retrive strain rate components: */
     578                exx=epsilon[0];
     579                eyy=epsilon[1];
     580                exy=epsilon[2];
     581                eff2 = exx*exx + eyy*eyy + exx*eyy + exy*exy ;
    551582
    552583                mu_prime=(1-n)/(2*n) * mu/eff2;

  • issm/trunk-jpl/src/c/objects/Materials/Matice.h

    r11322 r11332  
    6969                void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
    7070                void   GetViscosityDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
     71                void   GetViscosity2dDerivativeEpsSquare(double* pmu_prime, double* pepsilon);
    7172                double GetB();
    7273                double GetBbar();

  • issm/trunk-jpl/src/c/solvers/solver_newton.cpp

    r11327 r11332  
    6161                convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters);
    6262                MatFree(&Kff);VecFree(&pf);
    63                 count++;
    64                 if(converged==true)break;
     63                if(converged==true) break;
    6564                if(count>=max_nonlinear_iterations){
    6665                        _printf_(true,"   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
    … …  
    7170                SystemMatricesx(&Kff,NULL,&pf,NULL,&kmax,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
    7271                CreateJacobianMatrixx(&Jff,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,kmax);
    73                 CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);
     72                CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type);VecScale(ys,0.);
    7473                VecDuplicate(pf,&pJf);
    7574                MatMultPatch(Kff,uf,pJf); MatFree(&Kff);
    … …  
    8180                Mergesolutionfromftogx(&ug,uf,ys,femmodel->nodes,femmodel->parameters);VecFree(&ys);
    8281                InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ug);
     82
     83                count++;
    8384        }
    8485
Note: See TracChangeset for help on using the changeset viewer.