Changeset 15485


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r15483 r15485  
    241241ElementMatrix* Tria::CreateMassMatrix(void){
    242242
    243         /*constants: */
    244         const int    numdof=NDOF1*NUMVERTICES;
    245 
    246243        /* Intermediaries */
    247244        IssmDouble  D,Jdet;
    248245        IssmDouble  xyz_list[NUMVERTICES][3];
    249         IssmDouble  basis[numdof];
    250         GaussTria  *gauss = NULL;
    251 
    252         /*Initialize Element matrix*/
    253         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    254 
     246
     247        /*Fetch number of nodes and dof for this finite element*/
     248        int numnodes = this->NumberofNodes();
     249        int numdof   = numnodes*NDOF2;
     250
     251        /*Initialize Element matrix and vectors*/
     252        ElementMatrix* Ke    = new ElementMatrix(nodes,numnodes,this->parameters,NoneApproximationEnum);
     253        IssmDouble*    basis = xNew<IssmDouble>(numdof);
     254
     255        /*Retrieve all inputs and parameters*/
    255256        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    256257
    257258        /* Start looping on the number of gaussian points: */
    258         gauss=new GaussTria(2);
     259        GaussTria* gauss=new GaussTria(2);
    259260        for(int ig=gauss->begin();ig<gauss->end();ig++){
    260261
    261262                gauss->GaussPoint(ig);
    262263
    263                 GetNodalFunctions(&basis[0],gauss);
     264                GetNodalFunctions(basis,gauss);
    264265                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    265266                D=gauss->weight*Jdet;
    266267
    267                 TripleMultiply(&basis[0],1,3,1,
     268                TripleMultiply(basis,1,numnodes,1,
    268269                                        &D,1,1,0,
    269                                         &basis[0],1,3,0,
     270                                        basis,1,numnodes,0,
    270271                                        &Ke->values[0],1);
    271272        }
     
    273274        /*Clean up and return*/
    274275        delete gauss;
     276        xDelete<IssmDouble>(basis);
    275277        return Ke;
    276278}
     
    27812783        int numdof   = numnodes*NDOF2;
    27822784
    2783         /*Initialize Element matrix, vectors and Gaussian points*/
     2785        /*Initialize Element matrix and vectors*/
    27842786        ElementMatrix* Ke     = new ElementMatrix(nodes,numnodes,this->parameters,MacAyealApproximationEnum);
    27852787        IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
     
    27972799
    27982800        /* Start  looping on the number of gaussian points: */
    2799         GaussTria*     gauss  = new GaussTria(2);
     2801        GaussTria* gauss  = new GaussTria(2);
    28002802        for(int ig=gauss->begin();ig<gauss->end();ig++){
    28012803
Note: See TracChangeset for help on using the changeset viewer.