Changeset 15326

Timestamp:

06/22/13 21:18:51 (12 years ago)

Author:

Mathieu Morlighem

Message:

CHG: cosmetics'

Location:

issm/trunk-jpl/src/c/classes/Elements

Files:

• PentaRef.cpp (modified) (3 diffs)
• PentaRef.h (modified) (1 diff)
• TriaRef.cpp (modified) (1 diff)
• TriaRef.h (modified) (2 diffs)

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

@@ -986 +986 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctions{{{*/
+void PentaRef::GetNodalFunctions(IssmDouble* basis,GaussPenta* gauss){
+        /*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+        _assert_(basis);
+
+        switch(this->element_type){
+                case P1Enum:
+                case P1DGEnum:
+                        basis[0]=gauss->coord1*(1.-gauss->coord4)/2.;
+                        basis[1]=gauss->coord2*(1.-gauss->coord4)/2.;
+                        basis[2]=gauss->coord3*(1.-gauss->coord4)/2.;
+                        basis[3]=gauss->coord1*(1.+gauss->coord4)/2.;
+                        basis[4]=gauss->coord2*(1.+gauss->coord4)/2.;
+                        basis[5]=gauss->coord3*(1.+gauss->coord4)/2.;
+                        return;
+                case MINIEnum:
+                        basis[0]=gauss->coord1*(1.-gauss->coord4)/2.;
+                        basis[1]=gauss->coord2*(1.-gauss->coord4)/2.;
+                        basis[2]=gauss->coord3*(1.-gauss->coord4)/2.;
+                        basis[3]=gauss->coord1*(1.+gauss->coord4)/2.;
+                        basis[4]=gauss->coord2*(1.+gauss->coord4)/2.;
+                        basis[5]=gauss->coord3*(1.+gauss->coord4)/2.;
+                        basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+gauss->coord4)*(1.-gauss->coord4);
+                        return;
+                default:
+                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+        }
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivatives{{{*/
+void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
+
+        /*This routine returns the values of the nodal functions derivatives  (with respect to the 
+         * actual coordinate system): */
+        IssmDouble    Jinv[3][3];
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes();
+
+        /*Get nodal functions derivatives in reference triangle*/
+        IssmDouble* dbasis_ref=xNew<IssmDouble>(3*numnodes);
+        GetNodalFunctionsDerivativesReference(dbasis_ref,gauss); 
+
+        /*Get Jacobian invert: */
+        GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
+
+        /*Build dbasis: 
+         *
+         * [dhi/dx]= Jinv'*[dhi/dr]
+         * [dhi/dy]        [dhi/ds]
+         * [dhi/dz]        [dhi/dzeta]
+         */
+
+        for(int i=0;i<numnodes;i++){
+                dbasis[numnodes*0+i]=Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i]+Jinv[0][2]*dbasis_ref[2*numnodes+i];
+                dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i]+Jinv[1][2]*dbasis_ref[2*numnodes+i];
+                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];
+        }
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctionsDerivativesReference{{{*/
+void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussPenta* gauss){
+
+        /*This routine returns the values of the nodal functions derivatives  (with respect to the 
+         * natural coordinate system) at the gaussian point. */
+
+        _assert_(dbasis && gauss);
+
+        /*Get current coordinates in reference element*/
+        IssmDouble r=gauss->coord2-gauss->coord1;
+        IssmDouble s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0);
+        IssmDouble zeta=gauss->coord4;
+
+        switch(this->element_type){
+                case P1Enum: case P1DGEnum:
+                        /*Nodal function 1*/
+                        dbasis[NUMNODESP1*0+0]=-0.5*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESP1*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESP1*2+0]=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 2*/
+                        dbasis[NUMNODESP1*0+1]=0.5*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESP1*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESP1*2+1]=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 3*/
+                        dbasis[NUMNODESP1*0+2]=0;
+                        dbasis[NUMNODESP1*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESP1*2+2]=-0.5*(SQRT3/3.0*s+ONETHIRD);
+                        /*Nodal function 4*/
+                        dbasis[NUMNODESP1*0+3]=-0.5*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESP1*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESP1*2+3]=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 5*/
+                        dbasis[NUMNODESP1*0+4]=0.5*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESP1*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESP1*2+4]=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 6*/
+                        dbasis[NUMNODESP1*0+5]=0;
+                        dbasis[NUMNODESP1*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESP1*2+5]=0.5*(SQRT3/3.0*s+ONETHIRD);
+                case MINIEnum:
+                        /*Nodal function 1*/
+                        dbasis[NUMNODESMINI*0+0]=-0.5*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESMINI*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESMINI*2+0]=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 2*/
+                        dbasis[NUMNODESMINI*0+1]=0.5*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESMINI*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESMINI*2+1]=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 3*/
+                        dbasis[NUMNODESMINI*0+2]=0;
+                        dbasis[NUMNODESMINI*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
+                        dbasis[NUMNODESMINI*2+2]=-0.5*(SQRT3/3.0*s+ONETHIRD);
+                        /*Nodal function 4*/
+                        dbasis[NUMNODESMINI*0+3]=-0.5*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESMINI*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESMINI*2+3]=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 5*/
+                        dbasis[NUMNODESMINI*0+4]=0.5*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESMINI*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESMINI*2+4]=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
+                        /*Nodal function 6*/
+                        dbasis[NUMNODESMINI*0+5]=0;
+                        dbasis[NUMNODESMINI*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
+                        dbasis[NUMNODESMINI*2+5]=0.5*(SQRT3/3.0*s+ONETHIRD);
+                        /*Nodal function 7*/
+                        dbasis[NUMNODESMINI*0+6]=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
+                        dbasis[NUMNODESMINI*1+6]=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
+                        dbasis[NUMNODESMINI*2+6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
+                        return;
+                default:
+                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+        }
+
+
+}
+/*}}}*/
 /*FUNCTION PentaRef::GetNodalFunctionsMINI{{{*/
 void PentaRef::GetNodalFunctionsMINI(IssmDouble* l1l7, GaussPenta* gauss){
@@ -1136 +1274 @@
         *(dl1dl6+NUMNODESP1*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
         *(dl1dl6+NUMNODESP1*2+0)=-0.5*A1;
+        //dbasis[NUMNODESP1*0+0]=-0.5*(1.0-zeta)/2.0;
+        //dbasis[NUMNODESP1*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
+        //dbasis[NUMNODESP1*2+0]=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
 
         /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
@@ -1250 +1391 @@
         switch(this->element_type){
                 case P1Enum:   return NUMNODESP1;
+                case MINIEnum: return NUMNODESMINI;
                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
         }

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

@@ -22 +22 @@
 
                 /*Numerics*/
+                void GetNodalFunctions(IssmDouble* basis, GaussPenta* gauss);
+                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,GaussPenta* gauss);
+                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussPenta* gauss);
                 void GetNodalFunctionsP1(IssmDouble* l1l6, GaussPenta* gauss);
                 void GetNodalFunctionsMINI(IssmDouble* l1l7, GaussPenta* gauss);

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

@@ -491 +491 @@
         /*This routine returns the values of the nodal functions derivatives  (with respect to the 
          * actual coordinate system): */
-        IssmDouble    Jinv[NDOF2][NDOF2];
+        IssmDouble    Jinv[2][2];
 
         /*Fetch number of nodes for this finite element*/

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

@@ -27 +27 @@
                 void GetBprimeMacAyeal(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
                 void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
-                void GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss);
-                void GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss);
+                void GetBprimePrognostic(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
+                void GetBPrognostic(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
                 void GetBHydro(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
                 void GetBMacAyealFriction(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss);
@@ -39 +39 @@
                 void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
                 void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
-                void GetNodalFunctionsDerivatives(IssmDouble* basis,IssmDouble* xyz_list, GaussTria* gauss);
-                void GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss);
+                void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss);
+                void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
                 void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
                 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);