Changeset 15485

Timestamp:

07/11/13 15:44:54 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: mass matrix is now dynamically allocated for P2 elements

File:

• issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (modified) (4 diffs)

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -241 +241 @@
 ElementMatrix* Tria::CreateMassMatrix(void){
 
-        /*constants: */
-        const int    numdof=NDOF1*NUMVERTICES;
-
         /* Intermediaries */
         IssmDouble  D,Jdet;
         IssmDouble  xyz_list[NUMVERTICES][3];
-        IssmDouble  basis[numdof];
-        GaussTria  *gauss = NULL;
-
-        /*Initialize Element matrix*/
-        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
-
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = this->NumberofNodes();
+        int numdof   = numnodes*NDOF2;
+
+        /*Initialize Element matrix and vectors*/
+        ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
+        IssmDouble*    basis = xNew<IssmDouble>(numdof);
+
+        /*Retrieve all inputs and parameters*/
         GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
         /* Start looping on the number of gaussian points: */
-        gauss=new GaussTria(2);
+        GaussTria* gauss=new GaussTria(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){
 
                 gauss->GaussPoint(ig);
 
-                GetNodalFunctions(&basis[0],gauss);
+                GetNodalFunctions(basis,gauss);
                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
                 D=gauss->weight*Jdet;
 
-                TripleMultiply(&basis[0],1,3,1,
+                TripleMultiply(basis,1,numnodes,1,
                                         &D,1,1,0,
-                                        &basis[0],1,3,0,
+                                        basis,1,numnodes,0,
                                         &Ke->values[0],1);
         }
@@ -273 +274 @@
         /*Clean up and return*/
         delete gauss;
+        xDelete<IssmDouble>(basis);
         return Ke;
 }
@@ -2781 +2783 @@
         int numdof   = numnodes*NDOF2;
 
-        /*Initialize Element matrix, vectors and Gaussian points*/
+        /*Initialize Element matrix and vectors*/
         ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,MacAyealApproximationEnum);
         IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
@@ -2797 +2799 @@
 
         /* Start  looping on the number of gaussian points: */
-        GaussTria*     gauss  = new GaussTria(2);
+        GaussTria* gauss  = new GaussTria(2);
         for(int ig=gauss->begin();ig<gauss->end();ig++){