Changeset 17225

Timestamp:

02/06/14 21:28:00 (11 years ago)

Author:

seroussi

Message:

BUG: fixed problem in the case of no initial field

Location:

issm/trunk-jpl/src/c/classes

Files:

• Elements/PentaRef.cpp (modified) (5 diffs)
• Elements/PentaRef.h (modified) (3 diffs)
• IoModel.cpp (modified) (1 diff)

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

@@ -178 +178 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetNodalFunctions{{{*/
-void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in){
+/*FUNCTION PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss){{{*/
+void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss){
+        /*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+        _assert_(basis);
+        GetNodalFunctions(basis,gauss,this->element_type);
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement){{{*/
+void PentaRef::GetNodalFunctions(IssmDouble* basis,Gauss* gauss_in,int finiteelement){
         /*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -191 +200 @@
         IssmDouble zeta=gauss->coord4;
 
-        switch(this->element_type){
+        switch(finiteelement){
+                case P0Enum: 
+                        basis[0]=1.;
+                        return;
                 case P1Enum: case P1DGEnum:
                         basis[0]=gauss->coord1*(1.-zeta)/2.;
@@ -408 +420 @@
 
         switch(this->element_type){
+                case P0Enum: 
+                        /*Zero derivative*/
+                        dbasis[NUMNODESP0*0+0]   = 0.;
+                        dbasis[NUMNODESP0*1+0]   = 0.;
+                        dbasis[NUMNODESP0*2+0]   = 0.;
+                        return;
                 case P1Enum: case P1DGEnum:
                         /*Nodal function 1*/
@@ -996 +1014 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetInputValue{{{*/
+/*FUNCTION PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss){{{*/
 void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss){
-        /*P1 interpolation on Gauss point*/
-
-        /*intermediary*/
-        IssmDouble basis[6];
-
-        /*nodal functions: */
-        GetNodalFunctionsP1(&basis[0],gauss);
-
-        /*Assign output pointers:*/
-        *pvalue=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2]+basis[3]*plist[3]+basis[4]*plist[4]+basis[5]*plist[5];
+
+        GetInputValue(pvalue,plist,gauss,this->element_type);
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){{{*/
+void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,Gauss* gauss,int finiteelement){
+
+        /*Output*/
+        IssmDouble value =0.;
+
+        /*Fetch number of nodes for this finite element*/
+        int numnodes = this->NumberofNodes(finiteelement);
+
+        /*Get nodal functions*/
+        IssmDouble* basis=xNew<IssmDouble>(numnodes);
+        GetNodalFunctions(basis, gauss,finiteelement);
+
+        /*Calculate parameter for this Gauss point*/
+        for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];
+
+        /*Assign output pointer*/
+        xDelete<IssmDouble>(basis);
+        *pvalue = value;
 
 }
@@ -1022 +1054 @@
          *   p is a vector of size 3x1 already allocated.
          */
-        IssmDouble dbasis[3][NUMNODESP1];
-
-        /*Get nodal funnctions derivatives in actual coordinate system: */
-        GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
-
-        /*Assign output*/
-        p[0]=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2]+plist[3]*dbasis[0][3]+plist[4]*dbasis[0][4]+plist[5]*dbasis[0][5];
-        p[1]=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2]+plist[3]*dbasis[1][3]+plist[4]*dbasis[1][4]+plist[5]*dbasis[1][5];
-        p[2]=plist[0]*dbasis[2][0]+plist[1]*dbasis[2][1]+plist[2]*dbasis[2][2]+plist[3]*dbasis[2][3]+plist[4]*dbasis[2][4]+plist[5]*dbasis[2][5];
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::NumberofNodes{{{*/
+
+        /*Output*/
+        IssmDouble dpx=0.;
+        IssmDouble dpy=0.;
+        IssmDouble dpz=0.;
+
+        /*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);
+
+        /*Calculate parameter for this Gauss point*/
+        for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
+        for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];
+        for(int i=0;i<numnodes;i++) dpz += dbasis[2*numnodes+i]*plist[i];
+
+        /*Assign values*/
+        xDelete<IssmDouble>(dbasis);
+        p[0]=dpx;
+        p[1]=dpy;
+        p[2]=dpz;
+
+}
+/*}}}*/
+/*FUNCTION PentaRef::NumberofNodes(){{{*/
 int PentaRef::NumberofNodes(void){
 
-        switch(this->element_type){
+        return this->NumberofNodes(this->element_type);
+}
+/*}}}*/
+/*FUNCTION PentaRef::NumberofNodes(int finiteelement){{{*/
+int PentaRef::NumberofNodes(int finiteelement){
+
+        switch(finiteelement){
                 case P0Enum:                return NUMNODESP0;
                 case P1Enum:                return NUMNODESP1;

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

@@ -23 +23 @@
                 /*Numerics*/
                 void GetNodalFunctions(IssmDouble* basis, Gauss* gauss);
+                void GetNodalFunctions(IssmDouble* basis, Gauss* gauss,int finiteelement);
                 void GetNodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss);
                 void GetNodalFunctionsPressure(IssmDouble* basis, Gauss* gauss);
@@ -43 +44 @@
                 void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss);
                 void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss);
+                void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss,int finiteelement);
                 void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, Gauss* gauss);
 
@@ -48 +50 @@
                 void SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement);
                 int  NumberofNodes(void);
+                int  NumberofNodes(int finiteelement);
                 int  NumberofNodesVelocity(void);
                 int  NumberofNodesPressure(void);

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -1482 +1482 @@
         for(int i=0;i<elements->Size();i++){
                 Element* element=dynamic_cast<Element*>(elements->GetObjectByOffset(i));
-                if(!doublearray) element->InputUpdateFromConstant(default_value,vector_enum); 
+                if(!doublearray) element->AddInput(vector_enum,&default_value,P0Enum); 
                 else             element->InputCreate(doublearray,this,M,N,vector_layout,vector_enum,code);//we need i to index into elements.
         }