Changeset 26132


Ignore:
Timestamp:
03/23/21 11:58:46 (4 years ago)
Author:
tsantos
Message:

CHG: working on mono layer HO

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26090 r26132  
    27372737
    27382738        /*Intermediaries*/
    2739         IssmDouble  viscosity,Jdet;
     2739        IssmDouble  viscosity,Jdet,thickness,n;
    27402740        IssmDouble *xyz_list = NULL;
    27412741
    … …  
    27482748        /*Initialize Element matrix and vectors*/
    27492749        ElementMatrix* Ke     = basalelement->NewElementMatrix(MLHOApproximationEnum);
    2750         IssmDouble*    dbasis = xNew<IssmDouble>(3*numnodes);
     2750        IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes); // like SSA
     2751        IssmDouble*    basis  = xNew<IssmDouble>(numnodes); // like SSA
    27512752
    27522753        /*Retrieve all inputs and parameters*/
    27532754        element->GetVerticesCoordinates(&xyz_list);
    2754         Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
    2755         Input* vx_input      = element->GetInput(VxEnum);      _assert_(vx_input);
    2756         Input* vy_input      = element->GetInput(VyEnum);      _assert_(vy_input);
     2755        Input* thickness_input= element->GetInput(ThicknessEnum); _assert_(thickness_input);
     2756        Input* surface_input  = element->GetInput(SurfaceEnum);   _assert_(surface_input);
     2757        Input* vx_input       = element->GetInput(VxEnum);        _assert_(vx_input);
     2758        Input* vy_input       = element->GetInput(VyEnum);        _assert_(vy_input);
     2759        Input* n_input                   = element->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
    27572760
    27582761        /* Start  looping on the number of gaussian points: */
    … …  
    27642767                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    27652768                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss_base);
    2766 
     2769                element->NodalFunctions(basis, gauss);
     2770
     2771                thickness_input->GetInputValue(&thickness, gauss);
     2772                n_input->GetInputValue(&n,gauss);
    27672773                element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
    2768 
    2769                 for(int i=0;i<numnodes;i++){
     2774               
     2775                for(int i=0;i<numnodes;i++){//shape functions on tria element
    27702776                        for(int j=0;j<numnodes;j++){
    2771                                 Ke->values[2*i*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
     2777                                Ke->values[(i+0)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
    27722778                                                        4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
    2773                                                         );
    2774                                 Ke->values[2*i*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
     2779                                                        );//K11
     2780                                Ke->values[(i+0)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
     2781                                                        4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
     2782                                                        )*(n+1)/(n+2);//K12
     2783                                Ke->values[(i+0)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
    27752784                                                        2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
    2776                                                         );
    2777                                 Ke->values[(2*i+1)*2*numnodes+2*j] += gauss->weight*Jdet*viscosity*(
     2785                                                        );//K13
     2786                                Ke->values[(i+0)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
     2787                     2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
     2788                     )*(n+1)/(n+2);//K14
     2789                               
     2790                                Ke->values[(i+3)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
     2791                                                        4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
     2792                                                        )*(n+1)/(n+2);//K21
     2793                                Ke->values[(i+3)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
     2794                                                        4.*dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
     2795                                                        )*2*pow(n+1,2)/((2*n+3)*(n+2))
     2796                                                        +
     2797                                                        gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
     2798                                                        ;//K22
     2799                                Ke->values[(i+3)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
     2800                                                        2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
     2801                                                        )*(n+1)/(n+2);//K23
     2802                                Ke->values[(i+3)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
     2803                                                        2.*dbasis[1*numnodes+j]*dbasis[0*numnodes+i] + dbasis[0*numnodes+j]*dbasis[1*numnodes+i]
     2804                                                        )*2*pow(n+1,2)/((2*n+3)*(n+2));//K24
     2805                               
     2806                                Ke->values[(i+6)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
    27782807                                                        2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
    2779                                                         );
    2780                                 Ke->values[(2*i+1)*2*numnodes+2*j+1] += gauss->weight*Jdet*viscosity*(
    2781                                                         dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + 4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i]
    2782                                                         );
     2808                                                        );//K31
     2809                                Ke->values[(i+6)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
     2810                                                        2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
     2811                                                        )*(n+1)/(n+2);//K32
     2812                                Ke->values[(i+6)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
     2813                                                        4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
     2814                                                        );//K33
     2815                                Ke->values[(i+6)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
     2816                                                        4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
     2817                                                        )*(n+1)/(n+2);//K34
     2818
     2819                                Ke->values[(i+9)*2*2*numnodes+j+0] += gauss->weight*Jdet*viscosity*thickness*(
     2820                     2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
     2821                     )*(n+1)/(n+2);//K41
     2822                                Ke->values[(i+9)*2*2*numnodes+j+3] += gauss->weight*Jdet*viscosity*thickness*(
     2823                     2.*dbasis[0*numnodes+j]*dbasis[1*numnodes+i] + dbasis[1*numnodes+j]*dbasis[0*numnodes+i]
     2824                     )*2*pow(n+1,2)/((2*n+3)*(n+2));//K42
     2825                                Ke->values[(i+9)*2*2*numnodes+j+6] += gauss->weight*Jdet*viscosity*thickness*(
     2826                                                        4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
     2827                                                        )*(n+1)/(n+2);//K43
     2828                                Ke->values[(i+9)*2*2*numnodes+j+9] += gauss->weight*Jdet*viscosity*thickness*(
     2829                                                        4.*dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[0*numnodes+j]*dbasis[0*numnodes+i]
     2830                                                        )*2*pow(n+1,2)/((2*n+3)*(n+2))
     2831                                                        +
     2832                                                        gauss->weight*Jdet*viscosity*basis[j]*basis[i]*pow(n+1,2)/(thickness*(2*n+1))
     2833                                                        ;//K44
    27832834                        }
    27842835                }
    … …  
    27942845        xDelete<IssmDouble>(xyz_list);
    27952846        xDelete<IssmDouble>(dbasis);
     2847        xDelete<IssmDouble>(basis);
    27962848        return Ke;
    27972849}/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.