Changeset 5872


Ignore:
Timestamp:
09/17/10 16:27:59 (15 years ago)
Author:
Mathieu Morlighem
Message:

added some element matrices

Location:
issm/trunk/src/c/objects/Elements
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r5863 r5872  
    20042004/*}}}*/
    20052005/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattyn{{{1*/
    2006 void Penta::CreateKMatrixCouplingMacAyealPattyn( Mat Kgg){
    2007 
    2008         this->CreateKMatrixCouplingMacAyealPattynViscous(Kgg);
    2009         this->CreateKMatrixCouplingMacAyealPattynFriction(Kgg);
     2006ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattyn(void){
     2007
     2008        /*compute all stiffness matrices for this element*/
     2009        ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealPattynViscous();
     2010        ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction();
     2011        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
     2012
     2013        /*clean-up and return*/
     2014        delete Ke1;
     2015        delete Ke2;
     2016        return Ke;
    20102017}
    20112018/*}}}*/
    20122019/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynViscous{{{1*/
    2013 void Penta::CreateKMatrixCouplingMacAyealPattynViscous( Mat Kgg){
    2014 
    2015         /* local declarations */
    2016         int             i,j;
    2017 
    2018         /* node data: */
     2020ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynViscous(void){
     2021
    20192022        const int    NUMVERTICESm=3;  //MacAyealNUMVERTICES
    20202023        const int    numdofm=2*NUMVERTICESm;
    20212024        const int    NUMVERTICESp=6; //Pattyn NUMVERTICES
    20222025        const int    numdofp=2*NUMVERTICESp;
     2026        const int    numdoftotal=2*NDOF2*NUMVERTICES;
     2027        int             i,j,ig;
    20232028        double       xyz_list[NUMVERTICESp][3];
    2024         int*         doflistm=NULL;
    2025         int*         doflistp=NULL;
    2026 
    2027         /* 3d gaussian points: */
    2028         int     ig;
    20292029        GaussPenta *gauss=NULL;
    20302030        GaussTria  *gauss_tria=NULL;
    2031 
    2032         /* material data: */
    20332031        double viscosity; //viscosity
    20342032        double oldviscosity; //viscosity
    20352033        double newviscosity; //viscosity
    2036 
    2037         /* strain rate: */
    20382034        double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    20392035        double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
    2040 
    2041         /* matrices: */
    20422036        double B[3][numdofp];
    20432037        double Bprime[3][numdofm];
    … …  
    20472041        double DL[2][2]={0.0}; //for basal drag
    20482042        double DL_scalar;
    2049 
    2050         /* local element matrices: */
    20512043        double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix
    2052         double Ke_gg_transp[numdofm][numdofp]={0.0}; //local element stiffness matrix
    20532044        double Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point.
    20542045        double Jdet;
    2055 
    2056         /*friction: */
    20572046        double  alpha2_list[3];
    20582047        double  alpha2;
    2059 
    2060         /*parameters: */
    20612048        double viscosity_overshoot;
    20622049
    2063         /*Collapsed formulation: */
    2064         Tria*  tria     =NULL;
    2065         Penta* pentabase=NULL;
    2066 
    2067         /*retrieve some parameters: */
    2068         this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
    2069 
    2070         /*If on water, skip stiffness: */
    2071         if(IsOnWater())return;
    2072 
    20732050        /*Find penta on bed as pattyn must be coupled to the dofs on the bed: */
    2074         pentabase=GetBasalElement();
    2075         tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2051        Penta* pentabase=GetBasalElement();
     2052        Tria* tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2053
     2054        /*Initialize Element matrix and return if necessary*/
     2055        if(IsOnWater()) return NULL;
     2056        ElementMatrix* Ke1=pentabase->NewElementMatrix(MacAyealApproximationEnum);
     2057        ElementMatrix* Ke2=this->NewElementMatrix(PattynApproximationEnum);
     2058        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
     2059        delete Ke1; delete Ke2;
    20762060
    20772061        /* Get node coordinates and dof list: */
    20782062        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESp);
    2079         tria->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);  //Pattyn dof list
    2080         GetDofList(&doflistp,PattynApproximationEnum,GsetEnum); //MacAyeal dof list
    2081 
    2082         /*Retrieve all inputs we will be needing: */
     2063        this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
    20832064        Input* vx_input=inputs->GetInput(VxEnum);       ISSMASSERT(vx_input);
    20842065        Input* vy_input=inputs->GetInput(VyEnum);       ISSMASSERT(vy_input);
    … …  
    20942075                gauss->SynchronizeGaussTria(gauss_tria);
    20952076
    2096                 /*Get strain rate from velocity: */
    20972077                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    20982078                this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
    2099 
    2100                 /*Get viscosity: */
    21012079                matice->GetViscosity3d(&viscosity, &epsilon[0]);
    21022080                matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
    21032081
    2104                 /*Get B and Bprime matrices: */
    21052082                GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss);
    21062083                tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
    21072084
    2108                 /* Get Jacobian determinant: */
    21092085                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    2110 
    2111                 /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant
    2112                   onto this scalar matrix, so that we win some computational time: */
    2113 
    21142086                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    21152087                D_scalar=2*newviscosity*gauss->weight*Jdet;
    21162088                for (i=0;i<3;i++) D[i][i]=D_scalar;
    21172089
    2118                 /*  Do the triple product tB*D*Bprime: */
    21192090                TripleMultiply( &B[0][0],3,numdofp,1,
    21202091                                        &D[0][0],3,3,0,
    … …  
    21222093                                        &Ke_gg_gaussian[0][0],0);
    21232094
    2124                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
    21252095                for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    21262096        }
    2127 
    2128         /*Add Ke_gg and its transpose to global matrix Kgg: */
    2129         MatrixTranspose(&Ke_gg_transp[0][0],&Ke_gg[0][0],12,6);
    2130         MatSetValues(Kgg,numdofp,doflistp,numdofm,doflistm,(const double*)Ke_gg,ADD_VALUES);
    2131         MatSetValues(Kgg,numdofm,doflistm,numdofp,doflistp,(const double*)Ke_gg_transp,ADD_VALUES);
    2132 
    2133         xfree((void**)&doflistm);
    2134         xfree((void**)&doflistp);
     2097        for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
     2098        for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j][i];
     2099
     2100        /*Clean-up and return*/
    21352101        delete tria->matice; delete tria;
    21362102        delete gauss;
    21372103        delete gauss_tria;
     2104        return Ke;
    21382105}
    21392106/*}}}*/
    21402107/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynFriction{{{1*/
    2141 void Penta::CreateKMatrixCouplingMacAyealPattynFriction( Mat Kgg){
     2108ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynFriction(void){
    21422109
    21432110        /*Initialize Element matrix and return if necessary*/
    2144         if(IsOnWater() || IsOnShelf() || !IsOnBed()) return;
     2111        if(IsOnWater() || IsOnShelf() || !IsOnBed()) return NULL;
    21452112
    21462113        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    21472114        ElementMatrix* Ke=tria->CreateKMatrixCouplingMacAyealPattynFriction();
    2148         if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
    2149         delete Ke;
    21502115        delete tria->matice; delete tria;
     2116
     2117        return Ke;
    21512118}
    21522119/*}}}*/
    21532120/*FUNCTION Penta::CreateKMatrixCouplingPattynStokes{{{1*/
    2154 void Penta::CreateKMatrixCouplingPattynStokes( Mat Kgg){
     2121ElementMatrix* Penta::CreateKMatrixCouplingPattynStokes(void){
    21552122
    21562123        int        i,j;
    … …  
    22332200                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
    22342201                        delete Ke;
    2235                         CreateKMatrixCouplingMacAyealPattyn(Kgg);
     2202                        Ke=CreateKMatrixCouplingMacAyealPattyn();
     2203                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     2204                        delete Ke;
    22362205                        break;
    22372206                case PattynStokesApproximationEnum:
    … …  
    22422211                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
    22432212                        delete Ke;
    2244                         CreateKMatrixCouplingPattynStokes( Kgg);
     2213                        Ke=CreateKMatrixCouplingPattynStokes();
     2214                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     2215                        delete Ke;
    22452216                        break;
    22462217                default:

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

    r5863 r5872  
    123123                ElementMatrix* CreateKMatrixBalancedthickness(void);
    124124                ElementMatrix* CreateKMatrixBalancedvelocities(void);
    125                 void      CreateKMatrixCouplingMacAyealPattyn( Mat Kgg);
    126                 void      CreateKMatrixCouplingMacAyealPattynViscous( Mat Kgg);
    127                 void      CreateKMatrixCouplingMacAyealPattynFriction( Mat Kgg);
    128                 void      CreateKMatrixCouplingPattynStokes( Mat Kgg);
     125                ElementMatrix* CreateKMatrixCouplingMacAyealPattyn(void);
     126                ElementMatrix* CreateKMatrixCouplingMacAyealPattynViscous(void);
     127                ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
     128                ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
    129129                void      CreateKMatrixDiagnosticHoriz( Mat Kgg);
    130130                ElementMatrix* CreateKMatrixDiagnosticHutter(void);

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

    r5865 r5872  
    29432943        }
    29442944
    2945 
    29462945        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j];
    29472946        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j];
Note: See TracChangeset for help on using the changeset viewer.