Changeset 5879

Timestamp:

09/17/10 17:42:12 (15 years ago)

Author:

Mathieu Morlighem

Message:

final clean

Location:

issm/trunk/src/c/objects/Elements

Files:

• Penta.cpp (modified) (13 diffs)
• Penta.h (modified) (1 diff)
• Tria.cpp (modified) (4 diffs)
• Tria.h (modified) (1 diff)

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

@@ -700 +700 @@
                 case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
                         Ke=CreateKMatrixDiagnosticHoriz();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 case DiagnosticHutterAnalysisEnum:
                         Ke=CreateKMatrixDiagnosticHutter();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 case DiagnosticVertAnalysisEnum:
                         Ke=CreateKMatrixDiagnosticVert();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
                         Ke=CreateKMatrixSlope();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 case PrognosticAnalysisEnum:
                         Ke=CreateKMatrixPrognostic();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 case BalancedthicknessAnalysisEnum:
                         Ke=CreateKMatrixBalancedthickness();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 case BalancedvelocitiesAnalysisEnum:
                         Ke=CreateKMatrixBalancedvelocities();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 case ThermalAnalysisEnum:
-                        CreateKMatrixThermal( Kgg);
+                        Ke=CreateKMatrixThermal();
                         break;
                 case MeltingAnalysisEnum:
                         Ke=CreateKMatrixMelting();
-                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
-                        delete Ke;
                         break;
                 default:
                         ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
+        }
+
+        /*Add to global matrix*/
+        if(Ke){
+                Ke->AddToGlobal(Kgg,Kff,Kfs);
+                delete Ke;
         }
 }
@@ -2840 +2830 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrixThermal {{{1*/
-void  Penta::CreateKMatrixThermal(Mat Kgg){
-
-        /* local declarations */
+ElementMatrix* Penta::CreateKMatrixThermal(void){
+
+        /*compute all stiffness matrices for this element*/
+        ElementMatrix* Ke1=CreateKMatrixThermalVolume();
+        ElementMatrix* Ke2=CreateKMatrixThermalShelf();
+        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+        printf("-------------- file: Penta.cpp line: %i\n",__LINE__); 
+        Ke->Echo();
+
+        /*clean-up and return*/
+        delete Ke1;
+        delete Ke2;
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixThermalVolume {{{1*/
+ElementMatrix* Penta::CreateKMatrixThermalVolume(void){
+
         int i,j;
+        int     ig;
         int found=0;
-
-        /* node data: */
         const int    numdof=NDOF1*NUMVERTICES;
         double       xyz_list[NUMVERTICES][3];
         int*         doflist=NULL;
-
-        /* gaussian points: */
-        int     ig;
         GaussPenta *gauss=NULL;
-
         double  K[2][2]={0.0};
-
         double  u,v,w;
-
-        /*matrices: */
-        double     K_terms[numdof][numdof]={0.0};
         double     Ke_gaussian_conduct[numdof][numdof];
         double     Ke_gaussian_advec[numdof][numdof];
@@ -2881 +2877 @@
         double     tBD_artdiff[3][numdof];
         double     tLD[numdof];
-
         double     Jdet;
-
-        /*Material properties: */
         double     gravity,rho_ice,rho_water;
         double     heatcapacity,thermalconductivity;
         double     mixed_layer_capacity,thermal_exchange_velocity;
-
-        /*parameters: */
         double dt,epsvel;
         bool   artdiff;
-
-        /*Collapsed formulation: */
         Tria*  tria=NULL;
 
-        /*If on water, skip: */
-        if(IsOnWater())return;
-
-        /* Get node coordinates and dof list: */
+        /*Initialize Element matrix and return if necessary*/
+        if(IsOnWater()) return NULL;
+        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
+
+        /*Retrieve all inputs and parameters*/
         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-        // /*recovre material parameters: */
         rho_water=matpar->GetRhoWater();
         rho_ice=matpar->GetRhoIce();
@@ -2909 +2896 @@
         heatcapacity=matpar->GetHeatCapacity();
         thermalconductivity=matpar->GetThermalConductivity();
-
-        /*retrieve some parameters: */
         this->parameters->FindParam(&dt,DtEnum);
         this->parameters->FindParam(&artdiff,ArtDiffEnum);
         this->parameters->FindParam(&epsvel,EpsVelEnum);
-
         Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
         Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
@@ -2929 +2913 @@
                 /*Conduction: */
 
-                /*Get B_conduct matrix: */
                 GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); 
 
-                /*Build D: */
                 D_scalar=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
-
                 if(dt) D_scalar=D_scalar*dt;
 
@@ -2941 +2922 @@
                 D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
 
-                /*  Do the triple product B'*D*B: */
                 MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
                 MatrixMultiply(&tBD_conduct[0][0],numdof,3,0,&B_conduct[0][0],3,numdof,0,&Ke_gaussian_conduct[0][0],0);
@@ -2947 +2927 @@
                 /*Advection: */
 
-                /*Get B_advec and Bprime_advec matrices: */
                 GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); 
                 GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); 
 
-                //Build the D matrix
                 vx_input->GetParameterValue(&u, gauss);
                 vy_input->GetParameterValue(&v, gauss);
@@ -2957 +2935 @@
 
                 D_scalar=gauss->weight*Jdet;
-
                 if(dt) D_scalar=D_scalar*dt;
 
@@ -2964 +2941 @@
                 D[2][0]=0;         D[2][1]=0;         D[2][2]=D_scalar*w;
 
-                /*  Do the triple product B'*D*Bprime: */
                 MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
                 MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
 
                 /*Transient: */
+
                 if(dt){
                         GetNodalFunctionsP1(&L[0], gauss);
@@ -2974 +2951 @@
                         D_scalar=D_scalar;
 
-                        /*  Do the triple product L'*D*L: */
                         MatrixMultiply(&L[0],numdof,1,0,&D_scalar,1,1,0,&tLD[0],0);
                         MatrixMultiply(&tLD[0],numdof,1,0,&L[0],1,numdof,0,&Ke_gaussian_transient[0][0],0);
@@ -2983 +2959 @@
 
                 /*Artifficial diffusivity*/
+
                 if(artdiff){
                         /*Build K: */
@@ -2990 +2967 @@
                         K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
 
-                        /*Get B_artdiff: */
                         GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); 
 
-                        /*  Do the triple product B'*K*B: */
                         MatrixMultiply(&B_artdiff[0][0],2,numdof,1,&K[0][0],2,2,0,&tBD_artdiff[0][0],0);
                         MatrixMultiply(&tBD_artdiff[0][0],numdof,2,0,&B_artdiff[0][0],2,numdof,0,&Ke_gaussian_artdiff[0][0],0);
@@ -3001 +2976 @@
                 }
 
-                /*Add Ke_gaussian to pKe: */
-                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) K_terms[i][j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
-        }
-
-        /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
-
-        //Ice/ocean heat exchange flux on ice shelf base 
-        if(IsOnBed() && IsOnShelf()){
-
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrixThermal(Kgg);
-                delete tria->matice; delete tria;
-        }
-        
-        /*Free ressources:*/
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];
+        }
+
+        /*Clean up and return*/
         delete gauss;
-        xfree((void**)&doflist);
-
+        return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixThermalShelf {{{1*/
+ElementMatrix* Penta::CreateKMatrixThermalShelf(void){
+
+        if (!IsOnBed() || !IsOnShelf() || IsOnWater()) return NULL;
+
+        /*Call Tria function*/
+        Tria* tria=(Tria*)SpawnTria(0,1,2);
+        ElementMatrix* Ke=tria->CreateKMatrixThermal();
+        delete tria->matice; delete tria;
+
+        return Ke;
 }
 /*}}}*/

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

@@ -147 +147 @@
                 ElementMatrix* CreateKMatrixPrognostic(void);
                 ElementMatrix* CreateKMatrixSlope(void);
-                void      CreateKMatrixThermal(Mat Kggg);
+                ElementMatrix* CreateKMatrixThermal(void);
+                ElementMatrix* CreateKMatrixThermalVolume(void);
+                ElementMatrix* CreateKMatrixThermalShelf(void);
                 void      CreatePVectorBalancedthickness( Vec pg);
                 void      CreatePVectorBalancedvelocities( Vec pg);

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

@@ -3408 +3408 @@
 /*}}}*/
 /*FUNCTION Tria::CreateKMatrixThermal {{{1*/
-void  Tria::CreateKMatrixThermal(Mat Kgg){
-
+ElementMatrix* Tria::CreateKMatrixThermal(void){
+
+        const int    numdof=NDOF1*NUMVERTICES;
         int i,j;
-        
-        /* node data: */
-        const int    numdof=NDOF1*NUMVERTICES;
+        int     ig;
         double       xyz_list[NUMVERTICES][3];
-        int*         doflist=NULL;
-
         double mixed_layer_capacity;
         double thermal_exchange_velocity;
@@ -3422 +3419 @@
         double rho_ice;
         double heatcapacity;
-
-        int     ig;
         GaussTria *gauss=NULL;
-
-        /*matrices: */
         double  Jdet;
         double  K_terms[numdof][numdof]={0.0};
@@ -3433 +3426 @@
         double     tl1l2l3D[3];
         double  D_scalar;
-
-        /*parameters: */
         double dt;
 
-        /*retrieve some parameters: */
+        /*Initialize Element matrix and return if necessary*/
+        if(IsOnWater() || !IsOnShelf()) return NULL;
+        ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
+
+        /*Retrieve all inputs and parameters*/
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
         this->parameters->FindParam(&dt,DtEnum);
-
-        /* Get node coordinates and dof list: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-        //recover material parameters
         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
@@ -3457 +3447 @@
                 gauss->GaussPoint(ig);
                 
-                //Get the Jacobian determinant
-                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
-                
-                /*Get nodal functions values: */
                 GetNodalFunctions(&l1l2l3[0], gauss);
                                 
-                /*Calculate DL on gauss point */
+                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0], gauss);
                 D_scalar=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
-                if(dt){
-                        D_scalar=dt*D_scalar;
-                }
-
-                /*  Do the triple product tL*D*L: */
+                if(dt) D_scalar=dt*D_scalar;
+
                 MatrixMultiply(&l1l2l3[0],numdof,1,0,&D_scalar,1,1,0,&tl1l2l3D[0],0);
                 MatrixMultiply(&tl1l2l3D[0],numdof,1,0,&l1l2l3[0],1,numdof,0,&Ke_gaussian[0][0],0);
 
-                for(i=0;i<3;i++){
-                        for(j=0;j<3;j++){
-                                K_terms[i][j]+=Ke_gaussian[i][j];
-                        }
-                }
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
         }
         
-        /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
-
-
         /*Clean up and return*/
         delete gauss;
-        xfree((void**)&doflist);
+        return Ke;
 }
 /*}}}*/

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

@@ -138 +138 @@
                 ElementMatrix* CreateKMatrixPrognostic_DG(void);
                 ElementMatrix* CreateKMatrixSlope(void);
-                void      CreateKMatrixThermal(Mat Kgg);
+                ElementMatrix* CreateKMatrixThermal(void);
                 void      CreatePVectorBalancedthickness_DG(Vec pg);
                 void      CreatePVectorBalancedthickness_CG(Vec pg);