Changeset 5854


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

ElementMatrix for melting

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

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5853 r5854  
    687687void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
    688688
    689         /*retrive parameters: */
     689        /*retrieve parameters: */
    690690        int analysis_type;
     691        ElementMatrix* Ke=NULL;
    691692        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    692693
     
    722723                        break;
    723724                case MeltingAnalysisEnum:
    724                         CreateKMatrixMelting( Kgg);
     725                        Ke=CreateKMatrixMelting();
     726                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
     727                        delete Ke;
    725728                        break;
    726729                default:
     
    23212324          bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
    23222325          the stiffness matrix. */
    2323         if (!IsOnBed()) return NULL;
     2326        if (!IsOnBed() || IsOnWater()) return NULL;
    23242327
    23252328        /*Depth Averaging B*/
     
    28212824/*}}}*/
    28222825/*FUNCTION Penta::CreateKMatrixMelting {{{1*/
    2823 void  Penta::CreateKMatrixMelting(Mat Kgg){
    2824 
    2825         Tria* tria=NULL;
    2826 
    2827         /*If on water, skip: */
    2828         if(IsOnWater())return;
    2829 
    2830         if (!IsOnBed()){
    2831                 return;
    2832         }
    2833         else{
    2834 
    2835                 tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    2836                 tria->CreateKMatrixMelting(Kgg);
    2837                 delete tria->matice; delete tria;
    2838                 return;
    2839         }
     2826ElementMatrix* Penta::CreateKMatrixMelting(void){
     2827
     2828        if (!IsOnBed() || IsOnWater()) return NULL;
     2829
     2830        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
     2831        ElementMatrix* Ke=tria->CreateKMatrixMelting();
     2832
     2833        delete tria->matice; delete tria;
     2834        return Ke;
    28402835}
    28412836/*}}}*/
  • issm/trunk/src/c/objects/Elements/Penta.h

    r5853 r5854  
    140140                ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
    141141                void      CreateKMatrixDiagnosticVert( Mat Kgg);
    142                 void      CreateKMatrixMelting(Mat Kggg);
     142                ElementMatrix* CreateKMatrixMelting(void);
    143143                void      CreateKMatrixPrognostic(Mat Kggg);
    144144                void      CreateKMatrixSlope(Mat Kggg);
  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r5852 r5854  
    31543154/*}}}*/
    31553155/*FUNCTION Tria::CreateKMatrixMelting {{{1*/
    3156 void  Tria::CreateKMatrixMelting(Mat Kgg){
    3157 
    3158         /*indexing: */
    3159         int i,j;
     3156ElementMatrix* Tria::CreateKMatrixMelting(void){
    31603157
    31613158        const int  numdof=NUMVERTICES*NDOF1;
    3162         int*       doflist=NULL;
    3163 
    3164         /*Grid data: */
     3159        int i,j,ig;
    31653160        double     xyz_list[NUMVERTICES][3];
    3166 
    3167         /*Material constants */
    31683161        double     heatcapacity,latentheat;
    3169 
    3170         /* gaussian points: */
    3171         int     ig;
    31723162        GaussTria *gauss=NULL;
    3173 
    3174         /*matrices: */
    31753163        double     Jdet;
    31763164        double     D_scalar;
     
    31803168        double     Ke_gaussian[numdof][numdof]={0.0};
    31813169
    3182         /*Recover constants of ice */
     3170        /*Initialize Element matrix and return if necessary*/
     3171        if(IsOnWater()) return NULL;
     3172        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
     3173
     3174        /*Retrieve all inputs and parameters*/
     3175        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    31833176        latentheat=matpar->GetLatentHeat();
    31843177        heatcapacity=matpar->GetHeatCapacity();
    3185 
    3186         /* Get node coordinates and dof list: */
    3187         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3188         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    3189 
    31903178
    31913179        /* Start looping on the number of gauss  (nodes on the bedrock) */
     
    31953183                gauss->GaussPoint(ig);
    31963184
     3185                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    31973186                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
    3198 
    3199                 /*Get L matrix : */
    3200                 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    3201 
    3202                 /*Calculate DL on gauss point */
    32033187                D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
    32043188
    3205                 /*  Do the triple product tL*D*L: */
    32063189                MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
    32073190                MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian[0][0],0);
    32083191
    3209                 for(i=0;i<NUMVERTICES;i++){
    3210                         for(j=0;j<NUMVERTICES;j++){
    3211                                 K_terms[i][j]+=Ke_gaussian[i][j];
    3212                         }
    3213                 }
    3214         }
    3215 
    3216         /*Add Ke_gg to global matrix Kgg: */
    3217         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
     3192                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
     3193        }
    32183194
    32193195        /*Clean up and return*/
    32203196        delete gauss;
    3221         xfree((void**)&doflist);
     3197        return Ke;
    32223198}
    32233199/*}}}*/
  • issm/trunk/src/c/objects/Elements/Tria.h

    r5852 r5854  
    133133                ElementMatrix* CreateKMatrixDiagnosticHutter(void);
    134134                void      CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);
    135                 void      CreateKMatrixMelting(Mat Kgg);
     135                ElementMatrix* CreateKMatrixMelting(void);
    136136                ElementMatrix* CreateKMatrixPrognostic(void);
    137137                ElementMatrix* CreateKMatrixPrognostic_CG(void);
Note: See TracChangeset for help on using the changeset viewer.