Changeset 5854

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:

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

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

@@ -687 +687 @@
 void  Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){
 
-        /*retrive parameters: */
+        /*retrieve parameters: */
         int analysis_type;
+        ElementMatrix* Ke=NULL;
         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 
@@ -722 +723 @@
                         break;
                 case MeltingAnalysisEnum:
-                        CreateKMatrixMelting( Kgg);
+                        Ke=CreateKMatrixMelting();
+                        if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
+                        delete Ke;
                         break;
                 default:
@@ -2321 +2324 @@
           bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
           the stiffness matrix. */
-        if (!IsOnBed()) return NULL;
+        if (!IsOnBed() || IsOnWater()) return NULL;
 
         /*Depth Averaging B*/
@@ -2821 +2824 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrixMelting {{{1*/
-void  Penta::CreateKMatrixMelting(Mat Kgg){
-
-        Tria* tria=NULL;
-
-        /*If on water, skip: */
-        if(IsOnWater())return;
-
-        if (!IsOnBed()){
-                return;
-        }
-        else{
-
-                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-                tria->CreateKMatrixMelting(Kgg);
-                delete tria->matice; delete tria;
-                return;
-        }
+ElementMatrix* Penta::CreateKMatrixMelting(void){
+
+        if (!IsOnBed() || IsOnWater()) return NULL;
+
+        Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+        ElementMatrix* Ke=tria->CreateKMatrixMelting();
+
+        delete tria->matice; delete tria;
+        return Ke;
 }
 /*}}}*/

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

@@ -140 +140 @@
                 ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
                 void      CreateKMatrixDiagnosticVert( Mat Kgg);
-                void      CreateKMatrixMelting(Mat Kggg);
+                ElementMatrix* CreateKMatrixMelting(void);
                 void      CreateKMatrixPrognostic(Mat Kggg);
                 void      CreateKMatrixSlope(Mat Kggg);

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

@@ -3154 +3154 @@
 /*}}}*/
 /*FUNCTION Tria::CreateKMatrixMelting {{{1*/
-void  Tria::CreateKMatrixMelting(Mat Kgg){
-
-        /*indexing: */
-        int i,j;
+ElementMatrix* Tria::CreateKMatrixMelting(void){
 
         const int  numdof=NUMVERTICES*NDOF1;
-        int*       doflist=NULL;
-
-        /*Grid data: */
+        int i,j,ig;
         double     xyz_list[NUMVERTICES][3];
-
-        /*Material constants */
         double     heatcapacity,latentheat;
-
-        /* gaussian points: */
-        int     ig;
         GaussTria *gauss=NULL;
-
-        /*matrices: */
         double     Jdet;
         double     D_scalar;
@@ -3180 +3168 @@
         double     Ke_gaussian[numdof][numdof]={0.0};
 
-        /*Recover constants of ice */
+        /*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);
         latentheat=matpar->GetLatentHeat();
         heatcapacity=matpar->GetHeatCapacity();
-
-        /* Get node coordinates and dof list: */
-        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
 
         /* Start looping on the number of gauss  (nodes on the bedrock) */
@@ -3195 +3183 @@
                 gauss->GaussPoint(ig);
 
+                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
-
-                /*Get L matrix : */
-                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
-                /*Calculate DL on gauss point */
                 D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
 
-                /*  Do the triple product tL*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[0][0],0);
 
-                for(i=0;i<NUMVERTICES;i++){
-                        for(j=0;j<NUMVERTICES;j++){
-                                K_terms[i][j]+=Ke_gaussian[i][j];
-                        }
-                }
-        }
-
-        /*Add Ke_gg to global matrix Kgg: */
-        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
+                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j];
+        }
 
         /*Clean up and return*/
         delete gauss;
-        xfree((void**)&doflist);
+        return Ke;
 }
 /*}}}*/

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

@@ -133 +133 @@
                 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
                 void      CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);
-                void      CreateKMatrixMelting(Mat Kgg);
+                ElementMatrix* CreateKMatrixMelting(void);
                 ElementMatrix* CreateKMatrixPrognostic(void);
                 ElementMatrix* CreateKMatrixPrognostic_CG(void);