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

Timestamp:

07/25/13 14:54:51 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: spcvx vy and vy need to be scaled before they are marshalled, removed unit conversion in CreateConstraintsDiagnosticHoriz and preparing HO with
quadratic elements

File:

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

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

@@ -45 +45 @@
 /*FUNCTION PentaRef::SetElementType{{{*/
 void PentaRef::SetElementType(int type,int type_counter){
-
-        _assert_(type==P1Enum || type==P1DGEnum);
 
         /*initialize element type*/
@@ -143 +141 @@
          * where h is the interpolation function for node i.
          *
-         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
-         */
-
-        IssmDouble dbasis[3][NUMNODESP1];
-
-        /*Get dbasis in actual coordinate system: */
-        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+         * We assume B has been allocated already, of size: 5x(NDOF2*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++){
-                B[NDOF2*NUMNODESP1*0+NDOF2*i+0] = dbasis[0][i];
-                B[NDOF2*NUMNODESP1*0+NDOF2*i+1] = 0.;
-
-                B[NDOF2*NUMNODESP1*1+NDOF2*i+0] = 0.;
-                B[NDOF2*NUMNODESP1*1+NDOF2*i+1] = dbasis[1][i];
-
-                B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
-                B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
-
-                B[NDOF2*NUMNODESP1*3+NDOF2*i+0] = .5*dbasis[2][i];
-                B[NDOF2*NUMNODESP1*3+NDOF2*i+1] = 0.;
-
-                B[NDOF2*NUMNODESP1*4+NDOF2*i+0] = 0.;
-                B[NDOF2*NUMNODESP1*4+NDOF2*i+1] = .5*dbasis[2][i];
-        }
-
+        for(int i=0;i<numnodes;i++){
+                B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
+                B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
+                B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
+                B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
+                B[NDOF2*numnodes*2+NDOF2*i+0] = .5*dbasis[1*numnodes+i];
+                B[NDOF2*numnodes*2+NDOF2*i+1] = .5*dbasis[0*numnodes+i];
+                B[NDOF2*numnodes*3+NDOF2*i+0] = .5*dbasis[2*numnodes+i];
+                B[NDOF2*numnodes*3+NDOF2*i+1] = 0.;
+                B[NDOF2*numnodes*4+NDOF2*i+0] = 0.;
+                B[NDOF2*numnodes*4+NDOF2*i+1] = .5*dbasis[2*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
 /*FUNCTION PentaRef::GetBprimeHO {{{*/
-void PentaRef::GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){
+void PentaRef::GetBprimeHO(IssmDouble* B,IssmDouble* xyz_list,GaussPenta* gauss){
         /*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
          * For node i, Bi can be expressed in the actual coordinate system
          * by: 
          *       Bi=[ 2*dh/dx     dh/dy   ]
-         *                [   dh/dx    2*dh/dy  ]
-         *                [ dh/dy      dh/dx    ]
-         *                [ dh/dz         0     ]
-         *                [  0         dh/dz    ]
+         *          [   dh/dx    2*dh/dy  ]
+         *          [ dh/dy      dh/dx    ]
+         *          [ dh/dz         0     ]
+         *          [  0         dh/dz    ]
          * where h is the interpolation function for node i.
          *
-         * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
-         */
-        IssmDouble dbasis[3][NUMNODESP1];
-
-        /*Get dbasis in actual coordinate system: */
-        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss_coord);
+         * We assume B has been allocated already, of size: 5x(NDOF2*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 BPrime: */
-        for(int i=0;i<NUMNODESP1;i++){
-                B[NDOF2*NUMNODESP1*0+NDOF2*i+0]=2.*dbasis[0][i]; 
-                B[NDOF2*NUMNODESP1*0+NDOF2*i+1]=dbasis[1][i];
-
-                B[NDOF2*NUMNODESP1*1+NDOF2*i+0]=dbasis[0][i];
-                B[NDOF2*NUMNODESP1*1+NDOF2*i+1]=2.*dbasis[1][i];
-
-                B[NDOF2*NUMNODESP1*2+NDOF2*i+0]=dbasis[1][i]; 
-                B[NDOF2*NUMNODESP1*2+NDOF2*i+1]=dbasis[0][i]; 
-
-                B[NDOF2*NUMNODESP1*3+NDOF2*i+0]=dbasis[2][i]; 
-                B[NDOF2*NUMNODESP1*3+NDOF2*i+1]=0.;
-
-                B[NDOF2*NUMNODESP1*4+NDOF2*i+0]=0.;
-                B[NDOF2*NUMNODESP1*4+NDOF2*i+1]=dbasis[2][i]; 
-        }
+        for(int i=0;i<numnodes;i++){
+                B[NDOF2*numnodes*0+NDOF2*i+0]=2.*dbasis[0*numnodes+i]; 
+                B[NDOF2*numnodes*0+NDOF2*i+1]=dbasis[1*numnodes+i];
+                B[NDOF2*numnodes*1+NDOF2*i+0]=dbasis[0*numnodes+i];
+                B[NDOF2*numnodes*1+NDOF2*i+1]=2.*dbasis[1*numnodes+i];
+                B[NDOF2*numnodes*2+NDOF2*i+0]=dbasis[1*numnodes+i]; 
+                B[NDOF2*numnodes*2+NDOF2*i+1]=dbasis[0*numnodes+i]; 
+                B[NDOF2*numnodes*3+NDOF2*i+0]=dbasis[2*numnodes+i]; 
+                B[NDOF2*numnodes*3+NDOF2*i+1]=0.;
+                B[NDOF2*numnodes*4+NDOF2*i+0]=0.;
+                B[NDOF2*numnodes*4+NDOF2*i+1]=dbasis[2*numnodes+i]; 
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -1183 +1184 @@
                 dbasis[numnodes*2+i]=Jinv[2][0]*dbasis_ref[0*numnodes+i]+Jinv[2][1]*dbasis_ref[1*numnodes+i]+Jinv[2][2]*dbasis_ref[2*numnodes+i];
         }
+
+        /*Clean up*/
+        xDelete<IssmDouble>(dbasis_ref);
 
 }
@@ -1622 +1626 @@
         switch(this->element_type){
                 case P1Enum:   return NUMNODESP1;
+                case P2Enum:   return NUMNODESP2;
                 case MINIEnum: return NUMNODESP1b;
                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");