Changeset 15312

Timestamp:

06/22/13 07:37:35 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: changed BMacAyeeal and BprimeMacAyeal for quadratic elements

File:

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

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

@@ -81 +81 @@
          * For node i, Bi can be expressed in the actual coordinate system
          * by: 
-         *       Bi=[ dh/dx           0    ]
-         *          [   0           dh/dy  ]
-         *          [ 1/2*dh/dy  1/2*dh/dx ]
-         * where h is the interpolation function for node i.
-         *
-         * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
-         */
-
-        int i;
-        IssmDouble dbasis[NDOF2][NUMNODESP1];
-
-        /*Get dh1dh2dh3 in actual coordinate system: */
-        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
+         *       Bi=[ dN/dx           0    ]
+         *          [   0           dN/dy  ]
+         *          [ 1/2*dN/dy  1/2*dN/dx ]
+         * where N is the interpolation function for node i.
+         *
+         * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Get nodal functions*/
+        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
         /*Build B: */
-        for (i=0;i<NUMNODESP1;i++){
-                *(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
-                *(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0;
-                *(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0;
-                *(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
-                *(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i]; 
-                *(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][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]; 
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -205 +209 @@
          * For node i, Bi' can be expressed in the actual coordinate system
          * by: 
-         *       Bi_prime=[ 2*dh/dx    dh/dy ]
-         *                [   dh/dx  2*dh/dy ]
-         *                [   dh/dy    dh/dx ]
-         * where h is the interpolation function for node i.
-         *
-         * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODESP1)
-         */
-
-        /*Same thing in the actual coordinate system: */
-        IssmDouble dbasis[NDOF2][NUMNODESP1];
-
-        /*Get dh1dh2dh3 in actual coordinates system : */
-        GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
+         *       Bi_prime=[ 2*dN/dx    dN/dy ]
+         *                [   dN/dx  2*dN/dy ]
+         *                [   dN/dy    dN/dx ]
+         * where hNis the interpolation function for node i.
+         *
+         * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+         */
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Get nodal functions*/
+        IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+        GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
         /*Build B': */
         for (int i=0;i<NUMNODESP1;i++){
-                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i)   = 2.*dbasis[0][i];
-                *(Bprime+NDOF2*NUMNODESP1*0+NDOF2*i+1) = dbasis[1][i];
-                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i)   = dbasis[0][i];
-                *(Bprime+NDOF2*NUMNODESP1*1+NDOF2*i+1) = 2.*dbasis[1][i];
-                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i)   = dbasis[1][i];
-                *(Bprime+NDOF2*NUMNODESP1*2+NDOF2*i+1) = dbasis[0][i];
-        }
+                Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
+                Bprime[NDOF2*NUMNODESP1*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
+                Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+0] =    dbasis[0*numnodes+i];
+                Bprime[NDOF2*NUMNODESP1*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];
+                Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+0] =    dbasis[1*numnodes+i];
+                Bprime[NDOF2*NUMNODESP1*2+NDOF2*i+1] =    dbasis[0*numnodes+i];
+        }
+
+        /*Clean-up*/
+        xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
@@ -425 +433 @@
 
         /*Get nodal functions*/
-        IssmDouble* BasisFunctions=xNew<IssmDouble>(numnodes);
-        GetNodalFunctions(BasisFunctions,gauss);
+        IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
+        GetNodalFunctions(triabasis,gauss);
 
         switch(this->element_type){
                 case P1Enum:
                 case P1DGEnum:
-                        basis[0]=BasisFunctions[index1];
-                        basis[1]=BasisFunctions[index2];
-                        xDelete<IssmDouble>(BasisFunctions);
+                        basis[0]=triabasis[index1];
+                        basis[1]=triabasis[index2];
+                        xDelete<IssmDouble>(triabasis);
                         return;
                 case P2Enum:
                         _assert_(index2<index1);
-                        basis[0]=BasisFunctions[index1];
-                        basis[1]=BasisFunctions[index2];
-                        basis[2]=BasisFunctions[3+index2-1];
-                        xDelete<IssmDouble>(BasisFunctions);
+                        basis[0]=triabasis[index1];
+                        basis[1]=triabasis[index2];
+                        basis[2]=triabasis[3+index2-1];
+                        xDelete<IssmDouble>(triabasis);
                         return;
                 default:
@@ -446 +454 @@
         }
 
-        xDelete<IssmDouble>(BasisFunctions);
+        /*Clean up*/
+        xDelete<IssmDouble>(triabasis);
 }
 /*}}}*/