Changeset 15642

Timestamp:

07/26/13 11:31:45 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: extending higher order elements to other solutions

File:

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

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

@@ -144 +144 @@
          * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
          */
-
 
         /*Fetch number of nodes for this finite element*/
@@ -531 +530 @@
          */
 
-        /*Same thing in the actual coordinate system: */
-        IssmDouble basis[6];
-
-        /*Get dh1dh2dh3 in actual coordinates system : */
-        GetNodalFunctionsP1(basis, gauss);
-
-        /*Build B': */
-        for(int i=0;i<NUMNODESP1;i++){
-                B_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = basis[i];
-                B_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = basis[i];
-                B_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = basis[i];
-        }
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* basis=xNew<IssmDouble>(numnodes);
+        GetNodalFunctions(basis,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B_advec[numnodes*0+i] = basis[i];
+                B_advec[numnodes*1+i] = basis[i];
+                B_advec[numnodes*2+i] = basis[i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(basis);
 }
 /*}}}*/
@@ -555 +558 @@
          * where h is the interpolation function for node i.
          *
-         * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
-         */
-
-        /*Same thing in the actual coordinate system: */
-        IssmDouble dbasis[3][NUMNODESP1];
-
-        /*Get dh1dh2dh3 in actual coordinates system : */
-        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
-
-        /*Build B': */
-        for(int i=0;i<NUMNODESP1;i++){
-                B_conduct[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
-                B_conduct[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
-                B_conduct[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
-        }
+         * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+        /*Build B: */
+        for(int i=0;i<numnodes;i++){
+                B_conduct[numnodes*0+i] = dbasis[0*numnodes+i];
+                B_conduct[numnodes*1+i] = dbasis[1*numnodes+i];
+                B_conduct[numnodes*2+i] = dbasis[2*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -577 +584 @@
                 where hi is the interpolation function for node i.*/
 
-        /*Get dbasis in actual coordinate system: */
-        IssmDouble dbasis[3][NUMNODESP1];
-        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
         /*Build B: */
-        for(int i=0;i<NUMNODESP1;i++){
-                B[i] = dbasis[2][i];  
-        }
-
+        for(int i=0;i<numnodes;i++){
+                B[i] = dbasis[2*numnodes+i];  
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -598 +610 @@
          * where h is the interpolation function for node i.
          *
-         * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
-         */
-
-        /*Get nodal function derivatives in actual coordinates system : */
-        IssmDouble dbasis[3][NUMNODESP1];
-        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
+         * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Get nodal functions derivatives*/
+        IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
         /*Build B': */
-        for(int i=0;i<NUMNODESP1;i++){
-                Bprime_advec[NDOF1*NUMNODESP1*0+NDOF1*i] = dbasis[0][i];
-                Bprime_advec[NDOF1*NUMNODESP1*1+NDOF1*i] = dbasis[1][i];
-                Bprime_advec[NDOF1*NUMNODESP1*2+NDOF1*i] = dbasis[2][i];
-        }
+        for(int i=0;i<numnodes;i++){
+                Bprime_advec[numnodes*0+i] = dbasis[0*numnodes+i];
+                Bprime_advec[numnodes*1+i] = dbasis[1*numnodes+i];
+                Bprime_advec[numnodes*2+i] = dbasis[2*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -616 +634 @@
 void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
 
-        GetNodalFunctionsP1(B, gauss);
+        GetNodalFunctions(B,gauss);
 
 }
@@ -631 +649 @@
          ** We assume B has been allocated already, of size: 2 (2 x numnodes)
          **/
-
 
         /*Fetch number of nodes for this finite element*/