Changeset 5849 for issm/trunk


Ignore:
Timestamp:
09/16/10 14:26:44 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added ElementMatrix for Pattyn

Location:
issm/trunk/src/c
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r5848 r5849  
    22002200                        break;
    22012201                case PattynApproximationEnum:
    2202                         CreateKMatrixDiagnosticPattyn( Kgg);
     2202                        Ke=CreateKMatrixDiagnosticPattyn();
     2203                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     2204                        delete Ke;
    22032205                        break;
    22042206                case StokesApproximationEnum:
    … …  
    22112213                case MacAyealPattynApproximationEnum:
    22122214                        CreateKMatrixDiagnosticMacAyeal3d( Kgg);
    2213                         CreateKMatrixDiagnosticPattyn( Kgg);
    2214                         CreateKMatrixCouplingMacAyealPattyn( Kgg);
     2215                        Ke=CreateKMatrixDiagnosticPattyn();
     2216                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     2217                        delete Ke;
     2218                        CreateKMatrixCouplingMacAyealPattyn(Kgg);
    22152219                        break;
    22162220                case PattynStokesApproximationEnum:
    2217                         CreateKMatrixDiagnosticPattyn( Kgg);
     2221                        Ke=CreateKMatrixDiagnosticPattyn();
     2222                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     2223                        delete Ke;
    22182224                        CreateKMatrixDiagnosticStokes( Kgg);
    22192225                        CreateKMatrixCouplingPattynStokes( Kgg);
    … …  
    24152421/*}}}*/
    24162422/*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{1*/
    2417 void Penta::CreateKMatrixDiagnosticPattyn( Mat Kgg){
     2423ElementMatrix* Penta::CreateKMatrixDiagnosticPattyn(void){
     2424
     2425        /*compute all stiffness matrices for this element*/
     2426        ElementMatrix* Ke1=CreateKMatrixDiagnosticPattynViscous();
     2427        ElementMatrix* Ke2=CreateKMatrixDiagnosticPattynFriction();
     2428        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     2429
     2430        /*clean-up and return*/
     2431        delete Ke1;
     2432        delete Ke2;
     2433        return Ke;
     2434}
     2435/*}}}*/
     2436/*FUNCTION Penta::CreateKMatrixDiagnosticPattynViscous{{{1*/
     2437ElementMatrix* Penta::CreateKMatrixDiagnosticPattynViscous(void){
    24182438
    24192439        /* local declarations */
    2420         int             i,j;
    2421 
    2422         /* node data: */
    24232440        const int    numdof=2*NUMVERTICES;
     2441        int             i,j,ig;
    24242442        double       xyz_list[NUMVERTICES][3];
    2425         int*         doflist=NULL;
    2426 
    2427         /* 3d gaussian points: */
    2428         int     ig;
    24292443        GaussPenta *gauss=NULL;
    2430 
    2431         /* material data: */
    24322444        double viscosity; //viscosity
    24332445        double oldviscosity; //viscosity
    24342446        double newviscosity; //viscosity
    2435 
    2436         /* strain rate: */
     2447        double viscosity_overshoot;
    24372448        double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    24382449        double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    2439 
    2440         /* matrices: */
    24412450        double B[5][numdof];
    24422451        double Bprime[5][numdof];
    … …  
    24462455        double DL[2][2]={0.0}; //for basal drag
    24472456        double DL_scalar;
    2448 
    2449         /* local element matrices: */
    24502457        double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix
    24512458        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    24522459        double Jdet;
    2453 
    2454         /*slope: */
    24552460        double  slope[2]={0.0};
    24562461        double  slope_magnitude;
    2457 
    2458         /*friction: */
    24592462        double  alpha2_list[3];
    24602463        double  alpha2;
    2461 
    24622464        double MAXSLOPE=.06; // 6 %
    24632465        double MOUNTAINKEXPONENT=10;
    2464 
    2465         /*parameters: */
    2466         double viscosity_overshoot;
    2467 
    2468         /*Collapsed formulation: */
    24692466        Tria*  tria=NULL;
    24702467
    2471         /*retrieve some parameters: */
     2468        /*Initialize Element matrix and return if necessary*/
     2469        if(IsOnWater()) return NULL;
     2470        ElementMatrix* Ke=this->NewElementMatrix(PattynApproximationEnum);
     2471
     2472        /*Retrieve all inputs and parameters*/
     2473        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    24722474        this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
    2473 
    2474         /*If on water, skip stiffness: */
    2475         if(IsOnWater())return;
    2476 
    2477         /*Implement standard penta element: */
    2478 
    2479         /* Get node coordinates and dof list: */
    2480         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2481         GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    2482 
    2483         /*Retrieve all inputs we will be needing: */
    24842475        Input* vx_input=inputs->GetInput(VxEnum);       ISSMASSERT(vx_input);
    24852476        Input* vy_input=inputs->GetInput(VyEnum);       ISSMASSERT(vy_input);
    … …  
    24932484                gauss->GaussPoint(ig);
    24942485
    2495                 /*Get strain rate from velocity: */
    24962486                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    24972487                this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    2498 
    2499                 /*Get viscosity: */
    25002488                matice->GetViscosity3d(&viscosity, &epsilon[0]);
    25012489                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    25022490
    2503                 /*Get B and Bprime matrices: */
    25042491                GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
    25052492                GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
    25062493
    2507                 /* Get Jacobian determinant: */
    25082494                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2509 
    2510                 /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
    2511                   onto this scalar matrix, so that we win some computational time: */
    25122495
    25132496                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    … …  
    25152498                for (i=0;i<5;i++) D[i][i]=D_scalar;
    25162499
    2517                 /*  Do the triple product tB*D*Bprime: */
    25182500                TripleMultiply( &B[0][0],5,numdof,1,
    25192501                                        &D[0][0],5,5,0,
    … …  
    25212503                                        &Ke_gg_gaussian[0][0],0);
    25222504
    2523                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
    2524                 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    2525         }
    2526 
    2527 
    2528         /*Add Ke_gg to global matrix Kgg: */
    2529         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    2530 
    2531         //Deal with 2d friction at the bedrock interface
    2532         if(IsOnBed() && !IsOnShelf()){
    2533                 /*Build a tria element using the 3 grids of the base of the penta. Then use
    2534                  * the tria functionality to build a friction stiffness matrix on these 3
    2535                  * grids: */
    2536 
    2537                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2538                 tria->CreateKMatrixDiagnosticPattynFriction(Kgg);
    2539                 delete tria->matice; delete tria;
     2505                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
    25402506        }
    25412507
    25422508        /*Clean up and return*/
    25432509        delete gauss;
    2544         xfree((void**)&doflist);
     2510        return Ke;
     2511}
     2512/*}}}*/
     2513/*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{1*/
     2514ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){
     2515
     2516        /*Initialize Element matrix and return if necessary*/
     2517        if(IsOnWater() || IsOnShelf() || !IsOnBed()) return NULL;
     2518
     2519        /*Build a tria element using the 3 grids of the base of the penta. Then use
     2520         * the tria functionality to build a friction stiffness matrix on these 3
     2521         * grids: */
     2522        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2523        ElementMatrix* Ke=tria->CreateKMatrixDiagnosticPattynFriction();
     2524        delete tria->matice; delete tria;
     2525
     2526        /*clean-up and return*/
     2527        return Ke;
    25452528}
    25462529/*}}}*/
    … …  
    40844067        /*First, figure out size of doflist: */
    40854068        numdof=0;
    4086         for(i=0;i<3;i++){
     4069        for(i=0;i<NUMVERTICES;i++){
    40874070                ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    40884071                numdof+=ndof_list[i];
    … …  
    40954078                /*Populate: */
    40964079                count=0;
    4097                 for(i=0;i<3;i++){
     4080                for(i=0;i<NUMVERTICES;i++){
    40984081                        nodes[i]->GetDofList(&doflist[count],approximation_enum,setenum);
    40994082                        count+=ndof_list[i];

  • issm/trunk/src/c/objects/Elements/Penta.h

    r5848 r5849  
    129129                ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
    130130                void      CreateKMatrixDiagnosticMacAyeal3d(Mat Kgg);
    131                 void      CreateKMatrixDiagnosticPattyn( Mat Kgg);
     131                ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
     132                ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);
     133                ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
    132134                void      CreateKMatrixDiagnosticStokes( Mat Kgg);
    133135                void      CreateKMatrixDiagnosticVert( Mat Kgg);

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

    r5844 r5849  
    29782978/*}}}*/
    29792979/*FUNCTION Tria::CreateKMatrixDiagnosticPattynFriction {{{1*/
    2980 void  Tria::CreateKMatrixDiagnosticPattynFriction(Mat Kgg){
    2981 
    2982         /* local declarations */
     2980ElementMatrix* Tria::CreateKMatrixDiagnosticPattynFriction(void){
     2981
     2982        const int numdof   = NDOF2*NUMVERTICES;
    29832983        int       i,j;
     2984        int     ig;
    29842985        int       analysis_type;
    2985 
    2986         /* node data: */
    2987         const int numdof   = 2 *NUMVERTICES;
    29882986        double    xyz_list[NUMVERTICES][3];
    2989         int*      doflist=NULL;
    29902987        int       numberofdofspernode=2;
    2991 
    2992         /* gaussian points: */
    2993         int     ig;
    29942988        GaussTria *gauss=NULL;
    2995 
    2996         /* matrices: */
    29972989        double L[2][numdof];
    29982990        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
    29992991        double DL_scalar;
    3000 
    3001         /* local element matrices: */
    30022992        double Ke_gg[numdof][numdof]={0.0};
    30032993        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    3004        
    30052994        double Jdet;
    3006        
    3007         /*slope: */
    30082995        double  slope[2]={0.0,0.0};
    30092996        double  slope_magnitude;
    3010 
    3011         /*friction: */
    30122997        Friction *friction = NULL;
    30132998        double    alpha2;
    3014 
    30152999        double MAXSLOPE=.06; // 6 %
    30163000        double MOUNTAINKEXPONENT=10;
    3017 
    3018         /*inputs: */
    30193001        int  drag_type;
    30203002
    3021         /*retrive parameters: */
     3003        /*Initialize Element matrix and return if necessary*/
     3004        if(IsOnWater() || IsOnShelf()) return NULL;
     3005        ElementMatrix* Ke=this->NewElementMatrix(PattynApproximationEnum);
     3006
     3007        /*Retrieve all inputs and parameters*/
     3008        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    30223009        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3023 
    3024         /*retrieve inputs :*/
    30253010        inputs->GetParameterValue(&drag_type,DragTypeEnum);
    30263011        Input* surface_input=inputs->GetInput(SurfaceEnum); ISSMASSERT(surface_input);
    … …  
    30283013        Input* vy_input=inputs->GetInput(VyEnum);           ISSMASSERT(vy_input);
    30293014        Input* vz_input=inputs->GetInput(VzEnum);           ISSMASSERT(vz_input);
    3030        
    3031         /* Get node coordinates and dof list: */
    3032         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3033         GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
    3034 
    3035         if (IsOnShelf()){
    3036                 /*no friction, do nothing*/
    3037                 return;
    3038         }
    30393015
    30403016        /*build friction object, used later on: */
    … …  
    30483024                gauss->GaussPoint(ig);
    30493025
    3050                 /*Friction: */
    30513026                friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
    30523027
    … …  
    30603035                }
    30613036
    3062                 /* Get Jacobian determinant: */
     3037                GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
    30633038                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    3064 
    3065                 /*Get L matrix: */
    3066                 GetL(&L[0][0], &xyz_list[0][0], gauss,numberofdofspernode);
    3067 
    30683039               
    30693040                DL_scalar=alpha2*gauss->weight*Jdet;
    3070                 for (i=0;i<2;i++){
    3071                         DL[i][i]=DL_scalar;
    3072                 }
     3041                for (i=0;i<2;i++) DL[i][i]=DL_scalar;
    30733042               
    3074                 /*  Do the triple producte tL*D*L: */
    30753043                TripleMultiply( &L[0][0],2,numdof,1,
    30763044                                        &DL[0][0],2,2,0,
    … …  
    30783046                                        &Ke_gg_gaussian[0][0],0);
    30793047
    3080                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    3081 
    3082         }
    3083 
    3084         /*Add Ke_gg to global matrix Kgg: */
    3085         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     3048                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
     3049        }
    30863050
    30873051        /*Clean up and return*/
    30883052        delete gauss;
    30893053        delete friction;
    3090         xfree((void**)&doflist);
     3054        return Ke;
    30913055}       
    30923056/*}}}*/
    … …  
    51185082        /*First, figure out size of doflist: */
    51195083        numdof=0;
    5120         for(i=0;i<3;i++){
     5084        for(i=0;i<NUMVERTICES;i++){
    51215085                ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    51225086                numdof+=ndof_list[i];
    … …  
    51295093                /*Populate: */
    51305094                count=0;
    5131                 for(i=0;i<3;i++){
     5095                for(i=0;i<NUMVERTICES;i++){
    51325096                        nodes[i]->GetDofList(&doflist[count],approximation_enum,setenum);
    51335097                        count+=ndof_list[i];

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

    r5843 r5849  
    130130                void    CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs);
    131131                void      CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg);
    132                 void      CreateKMatrixDiagnosticPattynFriction(Mat Kgg);
     132                ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
    133133                ElementMatrix* CreateKMatrixDiagnosticHutter(void);
    134134                void      CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);

  • issm/trunk/src/c/shared/Elements/Paterson.cpp

    r3332 r5849  
    3232                B=pow((double)10,(double)8)*(-0.000292866376675*pow(T+50,3)+ 0.011672640664130*pow(T+50,2)  -0.325004442485481*(T+50)+  6.524779401948101);
    3333        }
    34         if((T>=-45.0) && (T<=-40.0)){
     34        else if((T>=-45.0) && (T<=-40.0)){
    3535                B=pow((double)10,(double)8)*(-0.000292866376675*pow(T+45,3)+ 0.007279645014004*pow(T+45,2)  -0.230243014094813*(T+45)+  5.154964909039554);
    3636        }
    37         if((T>=-40.0) && (T<=-35.0)){
     37        else if((T>=-40.0) && (T<=-35.0)){
    3838                B=pow((double)10,(double)8)*(0.000072737147457*pow(T+40,3)+  0.002886649363879*pow(T+40,2)  -0.179411542205399*(T+40)+  4.149132666831214);
    3939        }
    40         if((T>=-35.0) && (T<=-30.0)){
     40        else if((T>=-35.0) && (T<=-30.0)){
    4141                B=pow((double)10,(double)8)*(-0.000086144770023*pow(T+35,3)+ 0.003977706575736*pow(T+35,2)  -0.145089762507325*(T+35)+  3.333333333333331);
    4242        }
    43         if((T>=-30.0) && (T<=-25.0)){
     43        else if((T>=-30.0) && (T<=-25.0)){
    4444                B=pow((double)10,(double)8)*(-0.000043984685769*pow(T+30,3)+ 0.002685535025386*pow(T+30,2)  -0.111773554501713*(T+30)+  2.696559088937191);
    4545        }
    46         if((T>=-25.0) && (T<=-20.0)){
     46        else if((T>=-25.0) && (T<=-20.0)){
    4747                B=pow((double)10,(double)8)*(-0.000029799523463*pow(T+25,3)+ 0.002025764738854*pow(T+25,2)  -0.088217055680511*(T+25)+  2.199331606342181);
    4848        }
    49         if((T>=-20.0) && (T<=-15.0)){
     49        else if((T>=-20.0) && (T<=-15.0)){
    5050                B=pow((double)10,(double)8)*(0.000136920904777*pow(T+20,3)+  0.001578771886910*pow(T+20,2)  -0.070194372551690*(T+20)+  1.805165505978111);
    5151        }
    52         if((T>=-15.0) && (T<=-10.0)){
     52        else if((T>=-15.0) && (T<=-10.0)){
    5353                B=pow((double)10,(double)8)*(-0.000899763781026*pow(T+15,3)+ 0.003632585458564*pow(T+15,2)  -0.044137585824322*(T+15)+  1.510778053489523);
    5454        }
    55         if((T>=-10.0) && (T<=-5.0)){
     55        else if((T>=-10.0) && (T<=-5.0)){
    5656                B=pow((double)10,(double)8)*(0.001676964325070*pow(T+10,3)-  0.009863871256831*pow(T+10,2)  -0.075294014815659*(T+10)+  1.268434288203714);
    5757        }
    58         if((T>=-5.0) && (T<=-2.0)){
     58        else if((T>=-5.0) && (T<=-2.0)){
    5959                B=pow((double)10,(double)8)*(-0.003748937622487*pow(T+5,3)+0.015290593619213*pow(T+5,2)  -0.048160403003748*(T+5)+  0.854987973338348);
    6060        }
    61         if(T>=-2.0){
     61        else if(T>=-2.0){
    6262                B=pow((double)10,(double)8)*(-0.003748937622488*pow(T+2,3)-0.018449844983174*pow(T+2,2)  -0.057638157095631*(T+2)+  0.746900791092860);
    6363        }
Note: See TracChangeset for help on using the changeset viewer.