Changeset 5743


Ignore:
Timestamp:
09/10/10 10:20:08 (15 years ago)
Author:
Mathieu Morlighem
Message:

removed all double gauss

Location:
issm/trunk/src/c
Files:
25 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Container/Inputs.cpp

    r5682 r5743  
    4343
    4444/*Object management*/
    45 /*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type){{{1*/
    46 void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type){
     45/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/
     46void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){
    4747
    4848        vector<Object*>::iterator object;
    … …  
    7272/*}}}*/
    7373/*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/
    74 void Inputs::GetParameterValue(double* pvalue,GaussTria* gauss, int enum_type){
     74void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){
    7575
    7676        vector<Object*>::iterator object;
    … …  
    9999}
    100100/*}}}*/
    101 /*FUNCTION Inputs::GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type){{{1*/
    102 void Inputs::GetParameterValue(double* pvalue,GaussPenta* gauss, int enum_type){
    103 
    104         vector<Object*>::iterator object;
    105         Input* input=NULL;
    106         bool   found=false;
    107 
    108         /*Go through inputs and check whether any input with the same name is already in: */
    109         for ( object=objects.begin() ; object < objects.end(); object++ ){
    110 
    111                 input=(Input*)(*object);
    112                 if (input->EnumType()==enum_type){
    113                         found=true;
    114                         break;
    115                 }
    116         }
    117 
    118         if (!found){
    119                 /*we could not find an input with the correct enum type. No defaults values were provided,
    120                  * error out: */
    121                 ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));
    122         }
    123 
    124         /*Ok, we have an input if we made it here, request the input to return the values: */
    125         input->GetParameterValue(pvalue,gauss);
    126 
    127 }
    128 /*}}}*/
    129 /*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/
    130 void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){
    131 
    132         vector<Object*>::iterator object;
    133         Input* input=NULL;
    134         bool   found=false;
    135 
    136         /*Go through inputs and check whether any input with the same name is already in: */
    137         for ( object=objects.begin() ; object < objects.end(); object++ ){
    138 
    139                 input=(Input*)(*object);
    140                 if (input->EnumType()==enum_type){
    141                         found=true;
    142                         break;
    143                 }
    144         }
    145 
    146         if (!found){
    147                 /*we could not find an input with the correct enum type. Return the default value: */
    148                 *pvalue=defaultvalue;
    149         }
    150         else{
    151                 input->GetParameterValue(pvalue,gauss);
    152         }
    153 }
    154 /*}}}*/
    155101/*FUNCTION Inputs::GetParameterValue(bool* pvalue,int enum-type){{{1*/
    156102void Inputs::GetParameterValue(bool* pvalue,int enum_type){
    … …  
    263209        input->GetParameterAverage(pvalue);
    264210
    265 }
    266 /*}}}*/
    267 /*FUNCTION Inputs::GetParameterDerivativeValue{{{1*/
    268 void Inputs::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type){
    269 
    270         vector<Object*>::iterator object;
    271         Input* input=NULL;
    272         bool   found=false;
    273 
    274         /*Go through inputs and check whether any input with the same name is already in: */
    275         for ( object=objects.begin() ; object < objects.end(); object++ ){
    276 
    277                 input=(Input*)(*object);
    278                 if (input->EnumType()==enum_type){
    279                         found=true;
    280                         break;
    281                 }
    282         }
    283 
    284         if (!found){
    285                 /*we could not find an input with the correct enum type. No defaults values were provided,
    286                  * error out: */
    287                 ISSMERROR("could not find input with enum type %i (%s)",enum_type,EnumToString(enum_type));
    288         }
    289 
    290         /*Ok, we have an input if we made it here, request the input to return the value: */
    291         input->GetParameterDerivativeValue(derivativevalues,xyz_list,gauss);
    292211}
    293212/*}}}*/

  • issm/trunk/src/c/Container/Inputs.h

    r5682 r5743  
    4848                void GetParameterValue(int* pvalue,int enum_type);
    4949                void GetParameterValue(double* pvalue,int enum_type);
    50                 void GetParameterValue(double* pvalue,double* gauss,int enum_type);
    5150                void GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type);
    5251                void GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type);
    53                 void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue);
    54                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type);
    5552                /*}}}*/
    5653

  • issm/trunk/src/c/objects/Elements/Element.h

    r5727 r5743  
    3131                virtual void   CreatePVector(Vec pg)=0;
    3232                virtual void   GetSolutionFromInputs(Vec solution)=0;
    33                 virtual void   GetNodes(void** nodes)=0;
    3433                virtual int    GetNodeIndex(Node* node)=0;
    3534                virtual bool   GetShelf()=0;

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r5736 r5743  
    884884        this->results=new Results();
    885885
    886 }
    887 /*}}}*/
    888 /*FUNCTION Penta::GetNodes {{{1*/
    889 void  Penta::GetNodes(void** vpnodes){
    890 
    891         int i;
    892         Node** pnodes=NULL;
    893        
    894         /*virtual object: */
    895         pnodes=(Node**)vpnodes;
    896 
    897         for(i=0;i<6;i++){
    898                 pnodes[i]=nodes[i];
    899         }
    900886}
    901887/*}}}*/

  • issm/trunk/src/c/objects/Elements/Penta.h

    r5736 r5743  
    7878                void   CreatePVector(Vec pg);
    7979                void   DeleteResults(void);
    80                 void   GetNodes(void** nodes);
    8180                int    GetNodeIndex(Node* node);
    8281                bool   GetOnBed();
    … …  
    149148                void      GetDofList(int** pdoflist,int approximation_enum=0);
    150149                void      GetDofList1(int* doflist);
    151                 int       GetElementType(void);
    152                 void      GetNodalFunctions(double* l1l6, double* gauss_coord);
    153                 void      GetNodalFunctionsDerivatives(double* dh1dh6,double* xyz_list, double* gauss_coord);
    154                 void      GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord);
    155                 void      GetNodalFunctionsDerivativesReferenceStokes(double* dl1dl7,double* gauss_coord);
    156                 void      GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord);
    157                 void      GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
     150                int     GetElementType(void);
    158151                void    GetParameterListOnVertices(double* pvalue,int enumtype);
    159152                void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);

  • issm/trunk/src/c/objects/Elements/PentaRef.cpp

    r5726 r5743  
    5353/*Reference Element numerics*/
    5454/*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/
    55 void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss){
     55void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){
    5656        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    5757         * For grid i, Bi can be expressed in the actual coordinate system
    … …  
    9191/*}}}*/
    9292/*FUNCTION PentaRef::GetBPattyn {{{1*/
    93 void PentaRef::GetBPattyn(double* B, double* xyz_list, double* gauss){
     93void PentaRef::GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss){
    9494        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    9595         * For grid i, Bi can be expressed in the actual coordinate system
    … …  
    136136/*}}}*/
    137137/*FUNCTION PentaRef::GetBprimePattyn {{{1*/
    138 void PentaRef::GetBprimePattyn(double* B, double* xyz_list, double* gauss_coord){
     138void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){
    139139        /*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    140140         * For grid i, Bi can be expressed in the actual coordinate system
    … …  
    180180/*}}}*/
    181181/*FUNCTION PentaRef::GetBStokes {{{1*/
    182 void PentaRef::GetBStokes(double* B, double* xyz_list, double* gauss){
     182void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){
    183183
    184184        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
    … …  
    251251/*}}}*/
    252252/*FUNCTION PentaRef::GetBprimeStokes {{{1*/
    253 void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss){
     253void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){
    254254        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    255255         *      For grid i, Bi' can be expressed in the actual coordinate system
    … …  
    323323/*}}}*/
    324324/*FUNCTION PentaRef::GetBArtdiff {{{1*/
    325 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, double* gauss){
     325void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){
    326326        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    327327         * For grid i, Bi' can be expressed in the actual coordinate system
    … …  
    353353/*}}}*/
    354354/*FUNCTION PentaRef::GetBAdvec{{{1*/
    355 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, double* gauss){
     355void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){
    356356        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    357357         * For grid i, Bi' can be expressed in the actual coordinate system
    … …  
    385385/*}}}*/
    386386/*FUNCTION PentaRef::GetBConduct{{{1*/
    387 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, double* gauss){
     387void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){
    388388        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    389389         * For grid i, Bi' can be expressed in the actual coordinate system
    … …  
    417417/*}}}*/
    418418/*FUNCTION PentaRef::GetBVert{{{1*/
    419 void PentaRef::GetBVert(double* B, double* xyz_list, double* gauss){
     419void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){
    420420        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
    421421                where hi is the interpolation function for grid i.*/
    … …  
    437437/*}}}*/
    438438/*FUNCTION PentaRef::GetBprimeAdvec{{{1*/
    439 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, double* gauss){
     439void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){
    440440        /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    441441         * For grid i, Bi' can be expressed in the actual coordinate system
    … …  
    469469/*}}}*/
    470470/*FUNCTION PentaRef::GetBprimeVert{{{1*/
    471 void PentaRef::GetBprimeVert(double* B, double* xyz_list, double* gauss){
     471void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){
    472472        /* Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i*/
    473473
    … …  
    477477/*}}}*/
    478478/*FUNCTION PentaRef::GetLStokes {{{1*/
    479 void PentaRef::GetLStokes(double* LStokes, double* gauss){
     479void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){
    480480        /*
    481481         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    … …  
    507507
    508508        /*Get l1l2l3 in actual coordinate system: */
    509         l1l2l3[0]=gauss[0];
    510         l1l2l3[1]=gauss[1];
    511         l1l2l3[2]=gauss[2];
     509        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     510        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     511        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    512512
    513513        /*Build LStokes: */
    … …  
    574574/*}}}*/
    575575/*FUNCTION PentaRef::GetLprimeStokes {{{1*/
    576 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss){
     576void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
    577577
    578578        /*
    … …  
    604604
    605605        /*Get l1l2l3 in actual coordinate system: */
    606         l1l2l3[0]=gauss[0];
    607         l1l2l3[1]=gauss[1];
    608         l1l2l3[2]=gauss[2];
     606        l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
     607        l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
     608        l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    609609
    610610        GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
    … …  
    673673/*}}}*/
    674674/*FUNCTION PentaRef::GetJacobian {{{1*/
    675 void PentaRef::GetJacobian(double* J, double* xyz_list,double* gauss){
    676 
    677         const int NDOF3=3;
    678         int i,j;
    679 
    680         /*The Jacobian is constant over the element, discard the gaussian points.
    681          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    682 
    683         double A1,A2,A3; //area coordinates
    684         double xi,eta,zi; //parametric coordinates
    685 
    686         double x1,x2,x3,x4,x5,x6;
    687         double y1,y2,y3,y4,y5,y6;
    688         double z1,z2,z3,z4,z5,z6;
    689 
    690         /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
    691         A1=gauss[0];
    692         A2=gauss[1];
    693         A3=gauss[2];
    694 
    695         xi=A2-A1;
    696         eta=SQRT3*A3;
    697         zi=gauss[3];
    698 
    699         x1=*(xyz_list+3*0+0);
    700         x2=*(xyz_list+3*1+0);
    701         x3=*(xyz_list+3*2+0);
    702         x4=*(xyz_list+3*3+0);
    703         x5=*(xyz_list+3*4+0);
    704         x6=*(xyz_list+3*5+0);
    705 
    706         y1=*(xyz_list+3*0+1);
    707         y2=*(xyz_list+3*1+1);
    708         y3=*(xyz_list+3*2+1);
    709         y4=*(xyz_list+3*3+1);
    710         y5=*(xyz_list+3*4+1);
    711         y6=*(xyz_list+3*5+1);
    712 
    713         z1=*(xyz_list+3*0+2);
    714         z2=*(xyz_list+3*1+2);
    715         z3=*(xyz_list+3*2+2);
    716         z4=*(xyz_list+3*3+2);
    717         z5=*(xyz_list+3*4+2);
    718         z6=*(xyz_list+3*5+2);
    719 
    720         *(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
    721         *(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
    722         *(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);
    723 
    724         *(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
    725         *(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
    726         *(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);
    727 
    728         *(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
    729         *(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
    730         *(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
    731 
    732 }
    733 /*}}}*/
    734 /*FUNCTION PentaRef::GetJacobianDeterminant {{{1*/
    735 void PentaRef::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss){
    736         /*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take
    737          * the determinant of it: */
    738         double J[3][3];
    739 
    740         /*Get Jacobian*/
    741         GetJacobian(&J[0][0],xyz_list,gauss);
    742 
    743         /*Get Determinant*/
    744         Matrix3x3Determinant(Jdet,&J[0][0]);
    745         if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
    746 
    747 }
    748 /*}}}*/
    749 /*FUNCTION PentaRef::GetJacobianInvert {{{1*/
    750 void PentaRef::GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss){
    751 
    752         /*Jacobian*/
    753         double J[3][3];
    754 
    755         /*Call Jacobian routine to get the jacobian:*/
    756         GetJacobian(&J[0][0], xyz_list, gauss);
    757 
    758         /*Invert Jacobian matrix: */
    759         Matrix3x3Invert(Jinv,&J[0][0]);
    760 }
    761 /*}}}*/
    762 /*FUNCTION PentaRef::GetNodalFunctionsMINI{{{1*/
    763 void PentaRef::GetNodalFunctionsMINI(double* l1l7, double* gauss){
    764         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    765 
    766         l1l7[0]=gauss[0]*(1.0-gauss[3])/2.0;
    767         l1l7[1]=gauss[1]*(1.0-gauss[3])/2.0;
    768         l1l7[2]=gauss[2]*(1.0-gauss[3])/2.0;
    769         l1l7[3]=gauss[0]*(1.0+gauss[3])/2.0;
    770         l1l7[4]=gauss[1]*(1.0+gauss[3])/2.0;
    771         l1l7[5]=gauss[2]*(1.0+gauss[3])/2.0;
    772         l1l7[6]=27*gauss[0]*gauss[1]*gauss[2]*(1.0+gauss[3])*(1.0-gauss[3]);
    773 
    774 }
    775 /*}}}*/
    776 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{1*/
    777 void PentaRef::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss){
    778 
    779         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    780          * actual coordinate system): */
    781 
    782         int       i;
    783         const int numgrids = 7;
    784         double    dh1dh7_ref[3][numgrids];
    785         double    Jinv[3][3];
    786 
    787         /*Get derivative values with respect to parametric coordinate system: */
    788         GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss);
    789 
    790         /*Get Jacobian invert: */
    791         GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    792 
    793         /*Build dh1dh6:
    794          *
    795          * [dhi/dx]= Jinv'*[dhi/dr]
    796          * [dhi/dy]        [dhi/ds]
    797          * [dhi/dz]        [dhi/dzeta]
    798          */
    799 
    800         for (i=0;i<numgrids;i++){
    801                 *(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
    802                 *(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
    803                 *(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
    804         }
    805 
    806 }
    807 /*}}}*/
    808 /*FUNCTION PentaRef::GetNodalFunctionsMINIDerivativesReference{{{1*/
    809 void PentaRef::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss){
    810 
    811         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    812          * natural coordinate system) at the gaussian point. */
    813 
    814         int    numgrids=7; //six plus bubble grids
    815 
    816         double r=gauss[1]-gauss[0];
    817         double s=-3.0/SQRT3*(gauss[0]+gauss[1]-2.0/3.0);
    818         double zeta=gauss[3];
    819 
    820         /*First nodal function: */
    821         *(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
    822         *(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
    823         *(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
    824 
    825         /*Second nodal function: */
    826         *(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
    827         *(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
    828         *(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
    829 
    830         /*Third nodal function: */
    831         *(dl1dl7+numgrids*0+2)=0;
    832         *(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
    833         *(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);
    834 
    835         /*Fourth nodal function: */
    836         *(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
    837         *(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
    838         *(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);
    839 
    840         /*Fith nodal function: */
    841         *(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
    842         *(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
    843         *(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);
    844 
    845         /*Sixth nodal function: */
    846         *(dl1dl7+numgrids*0+5)=0;
    847         *(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
    848         *(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);
    849 
    850         /*Seventh nodal function: */
    851         *(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
    852         *(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
    853         *(dl1dl7+numgrids*2+6)=27*gauss[0]*gauss[1]*gauss[2]*(-2.0*zeta);
    854 
    855 }
    856 /*}}}*/
    857 /*FUNCTION PentaRef::GetNodalFunctionsP1 {{{1*/
    858 void PentaRef::GetNodalFunctionsP1(double* l1l6, double* gauss){
    859         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    860 
    861         l1l6[0]=gauss[0]*(1-gauss[3])/2.0;
    862         l1l6[1]=gauss[1]*(1-gauss[3])/2.0;
    863         l1l6[2]=gauss[2]*(1-gauss[3])/2.0;
    864         l1l6[3]=gauss[0]*(1+gauss[3])/2.0;
    865         l1l6[4]=gauss[1]*(1+gauss[3])/2.0;
    866         l1l6[5]=gauss[2]*(1+gauss[3])/2.0;
    867 
    868 }
    869 /*}}}*/
    870 /*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{1*/
    871 void PentaRef::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss){
    872 
    873         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    874          * actual coordinate system): */
    875         int       i;
    876         const int NDOF3    = 3;
    877         const int numgrids = 6;
    878         double    dh1dh6_ref[NDOF3][numgrids];
    879         double    Jinv[NDOF3][NDOF3];
    880 
    881         /*Get derivative values with respect to parametric coordinate system: */
    882         GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss);
    883 
    884         /*Get Jacobian invert: */
    885         GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    886 
    887         /*Build dh1dh3:
    888          *
    889          * [dhi/dx]= Jinv*[dhi/dr]
    890          * [dhi/dy]       [dhi/ds]
    891          * [dhi/dz]       [dhi/dn]
    892          */
    893 
    894         for (i=0;i<numgrids;i++){
    895                 *(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
    896                 *(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
    897                 *(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
    898         }
    899 
    900 }
    901 /*}}}*/
    902 /*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{1*/
    903 void PentaRef::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss){
    904 
    905         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    906          * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
    907 
    908         const int numgrids=6;
    909         double A1,A2,A3,z;
    910 
    911         A1=gauss[0]; ISSMASSERT(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
    912         A2=gauss[1]; ISSMASSERT(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
    913         A3=gauss[2]; ISSMASSERT(A3>=0 && A3<=1);//third area coordinate value  In term of xi and eta: A3=y/SQRT3;
    914         z=gauss[3];  ISSMASSERT(z>=-1 &&  z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2
    915 
    916         /*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
    917         *(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
    918         *(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
    919         *(dl1dl6+numgrids*2+0)=-0.5*A1;
    920 
    921         /*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
    922         *(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
    923         *(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
    924         *(dl1dl6+numgrids*2+1)=-0.5*A2;
    925 
    926         /*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
    927         *(dl1dl6+numgrids*0+2)=0.0;
    928         *(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
    929         *(dl1dl6+numgrids*2+2)=-0.5*A3;
    930 
    931         /*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
    932         *(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
    933         *(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
    934         *(dl1dl6+numgrids*2+3)=0.5*A1;
    935 
    936         /*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
    937         *(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
    938         *(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
    939         *(dl1dl6+numgrids*2+4)=0.5*A2;
    940 
    941         /*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
    942         *(dl1dl6+numgrids*0+5)=0.0;
    943         *(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
    944         *(dl1dl6+numgrids*2+5)=0.5*A3;
    945 }
    946 /*}}}*/
    947 /*FUNCTION PentaRef::GetParameterValue{{{1*/
    948 void PentaRef::GetParameterValue(double* pvalue,double* plist,double* gauss){
    949         /*P1 interpolation on Gauss point*/
    950 
    951         /*intermediary*/
    952         double l1l6[6];
    953 
    954         /*nodal functions: */
    955         GetNodalFunctionsP1(&l1l6[0],gauss);
    956 
    957         /*Assign output pointers:*/
    958         *pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[5]*plist[5];
    959 
    960 }
    961 /*}}}*/
    962 /*FUNCTION PentaRef::GetParameterDerivativeValue{{{1*/
    963 void PentaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss){
    964         /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
    965          *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
    966          *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
    967          *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
    968          *
    969          *   p is a vector of size 3x1 already allocated.
    970          */
    971         double dh1dh6[3][6];
    972 
    973         /*Get nodal funnctions derivatives in actual coordinate system: */
    974         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
    975 
    976         /*Assign output*/
    977         p[0]=plist[0]*dh1dh6[0][0]+plist[1]*dh1dh6[0][1]+plist[2]*dh1dh6[0][2]+plist[3]*dh1dh6[0][3]+plist[4]*dh1dh6[0][4]+plist[5]*dh1dh6[0][5];
    978         p[1]=plist[0]*dh1dh6[1][0]+plist[1]*dh1dh6[1][1]+plist[2]*dh1dh6[1][2]+plist[3]*dh1dh6[1][3]+plist[4]*dh1dh6[1][4]+plist[5]*dh1dh6[1][5];
    979         p[2]=plist[0]*dh1dh6[2][0]+plist[1]*dh1dh6[2][1]+plist[2]*dh1dh6[2][2]+plist[3]*dh1dh6[2][3]+plist[4]*dh1dh6[2][4]+plist[5]*dh1dh6[2][5];
    980 
    981 }
    982 /*}}}*/
    983 
    984 /*FUNCTION PentaRef::GetBMacAyealPattyn {{{1*/
    985 void PentaRef::GetBMacAyealPattyn(double* B, double* xyz_list, GaussPenta* gauss){
    986         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    987          * For grid i, Bi can be expressed in the actual coordinate system
    988          * by:
    989          *       Bi=[ dh/dx          0      ]
    990          *          [   0           dh/dy   ]
    991          *          [ 1/2*dh/dy  1/2*dh/dx  ]
    992          * where h is the interpolation function for grid i.
    993          *
    994          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    995          */
    996 
    997         int i;
    998         const int numgrids=6;
    999         const int NDOF3=3;
    1000         const int NDOF2=2;
    1001 
    1002         double dh1dh6[NDOF3][numgrids];
    1003 
    1004         /*Get dh1dh6 in actual coordinate system: */
    1005         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
    1006 
    1007         /*Build B: */
    1008         for (i=0;i<numgrids;i++){
    1009                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
    1010                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
    1011 
    1012                 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
    1013                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
    1014 
    1015                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
    1016                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
    1017 
    1018         }
    1019 
    1020 }
    1021 /*}}}*/
    1022 /*FUNCTION PentaRef::GetBPattyn {{{1*/
    1023 void PentaRef::GetBPattyn(double* B, double* xyz_list, GaussPenta* gauss){
    1024         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    1025          * For grid i, Bi can be expressed in the actual coordinate system
    1026          * by:
    1027          *       Bi=[ dh/dx          0      ]
    1028          *          [   0           dh/dy   ]
    1029          *          [ 1/2*dh/dy  1/2*dh/dx  ]
    1030          *          [ 1/2*dh/dz      0      ]
    1031          *          [  0         1/2*dh/dz  ]
    1032          * where h is the interpolation function for grid i.
    1033          *
    1034          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    1035          */
    1036 
    1037         int i;
    1038         const int numgrids=6;
    1039         const int NDOF3=3;
    1040         const int NDOF2=2;
    1041 
    1042         double dh1dh6[NDOF3][numgrids];
    1043 
    1044         /*Get dh1dh6 in actual coordinate system: */
    1045         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
    1046 
    1047         /*Build B: */
    1048         for (i=0;i<numgrids;i++){
    1049                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i];
    1050                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;
    1051 
    1052                 *(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
    1053                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];
    1054 
    1055                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i];
    1056                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i];
    1057 
    1058                 *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i];
    1059                 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    1060 
    1061                 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    1062                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i];
    1063         }
    1064 
    1065 }
    1066 /*}}}*/
    1067 /*FUNCTION PentaRef::GetBprimePattyn {{{1*/
    1068 void PentaRef::GetBprimePattyn(double* B, double* xyz_list, GaussPenta* gauss_coord){
    1069         /*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
    1070          * For grid i, Bi can be expressed in the actual coordinate system
    1071          * by:
    1072          *       Bi=[ 2*dh/dx     dh/dy   ]
    1073          *                [   dh/dx    2*dh/dy  ]
    1074          *                [ dh/dy      dh/dx    ]
    1075          *                [ dh/dz         0     ]
    1076          *                [  0         dh/dz    ]
    1077          * where h is the interpolation function for grid i.
    1078          *
    1079          * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    1080          */
    1081 
    1082         int i;
    1083         const int NDOF3=3;
    1084         const int NDOF2=2;
    1085         const int numgrids=6;
    1086 
    1087         double dh1dh6[NDOF3][numgrids];
    1088 
    1089         /*Get dh1dh6 in actual coordinate system: */
    1090         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
    1091 
    1092         /*Build BPrime: */
    1093         for (i=0;i<numgrids;i++){
    1094                 *(B+NDOF2*numgrids*0+NDOF2*i)=2.0*dh1dh6[0][i];
    1095                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh6[1][i];
    1096 
    1097                 *(B+NDOF2*numgrids*1+NDOF2*i)=dh1dh6[0][i];
    1098                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
    1099 
    1100                 *(B+NDOF2*numgrids*2+NDOF2*i)=dh1dh6[1][i];
    1101                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh6[0][i];
    1102 
    1103                 *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh6[2][i];
    1104                 *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    1105 
    1106                 *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    1107                 *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh6[2][i];
    1108         }
    1109 }
    1110 /*}}}*/
    1111 /*FUNCTION PentaRef::GetBStokes {{{1*/
    1112 void PentaRef::GetBStokes(double* B, double* xyz_list, GaussPenta* gauss){
    1113 
    1114         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID.
    1115          * For grid i, Bi can be expressed in the actual coordinate system
    1116          * by:          Bi=[ dh/dx          0              0       0  ]
    1117          *                                      [   0           dh/dy           0       0  ]
    1118          *                                      [   0             0           dh/dy     0  ]
    1119          *                                      [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
    1120          *                                      [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
    1121          *                                      [   0          1/2*dh/dz    1/2*dh/dy   0  ]
    1122          *                                      [   0             0             0       h  ]
    1123          *                                      [ dh/dx         dh/dy         dh/dz     0  ]
    1124          *      where h is the interpolation function for grid i.
    1125          *      Same thing for Bb except the last column that does not exist.
    1126          */
    1127 
    1128         int i;
    1129         const int calculationdof=3;
    1130         const int numgrids=6;
    1131         int DOFPERGRID=4;
    1132 
    1133         double dh1dh7[calculationdof][numgrids+1];
    1134         double l1l6[numgrids];
    1135 
    1136 
    1137         /*Get dh1dh7 in actual coordinate system: */
    1138         GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
    1139         GetNodalFunctionsP1(l1l6, gauss);
    1140 
    1141         /*Build B: */
    1142         for (i=0;i<numgrids+1;i++){
    1143                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
    1144                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
    1145                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
    1146                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
    1147                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
    1148                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
    1149                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
    1150                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
    1151                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
    1152                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i];
    1153                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i];
    1154                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
    1155                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
    1156                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
    1157                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
    1158                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
    1159                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
    1160                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
    1161                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
    1162                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
    1163                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
    1164                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
    1165                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
    1166                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
    1167         }
    1168 
    1169         for (i=0;i<numgrids;i++){ //last column not for the bubble function
    1170                 *(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
    1171                 *(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
    1172                 *(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
    1173                 *(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
    1174                 *(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
    1175                 *(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
    1176                 *(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
    1177                 *(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
    1178         }
    1179 
    1180 }
    1181 /*}}}*/
    1182 /*FUNCTION PentaRef::GetBprimeStokes {{{1*/
    1183 void PentaRef::GetBprimeStokes(double* B_prime, double* xyz_list, GaussPenta* gauss){
    1184         /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    1185          *      For grid i, Bi' can be expressed in the actual coordinate system
    1186          *      by:
    1187          *                              Bi'=[  dh/dx   0          0       0]
    1188          *                                       [   0      dh/dy      0       0]
    1189          *                                       [   0      0         dh/dz    0]
    1190          *                                       [  dh/dy   dh/dx      0       0]
    1191          *                                       [  dh/dz   0        dh/dx     0]
    1192          *                                       [   0      dh/dz    dh/dy     0]
    1193          *                                       [  dh/dx   dh/dy    dh/dz     0]
    1194          *                                       [   0      0          0       h]
    1195          *      where h is the interpolation function for grid i.
    1196          *
    1197          *      Same thing for the bubble fonction except that there is no fourth column
    1198          */
    1199 
    1200         int i;
    1201         const int calculationdof=3;
    1202         const int numgrids=6;
    1203         int DOFPERGRID=4;
    1204 
    1205         double dh1dh7[calculationdof][numgrids+1];
    1206         double l1l6[numgrids];
    1207 
    1208         /*Get dh1dh7 in actual coordinate system: */
    1209         GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
    1210 
    1211         GetNodalFunctionsP1(l1l6, gauss);
    1212 
    1213         /*B_primeuild B_prime: */
    1214         for (i=0;i<numgrids+1;i++){
    1215                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B_prime[0][DOFPERGRID*i]=dh1dh6[0][i];
    1216                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
    1217                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
    1218                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
    1219                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
    1220                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
    1221                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
    1222                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
    1223                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
    1224                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=dh1dh7[1][i];
    1225                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=dh1dh7[0][i];
    1226                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
    1227                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=dh1dh7[2][i];
    1228                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
    1229                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=dh1dh7[0][i];
    1230                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
    1231                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=dh1dh7[2][i];
    1232                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=dh1dh7[1][i];
    1233                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=dh1dh7[0][i];
    1234                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=dh1dh7[1][i];
    1235                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=dh1dh7[2][i];
    1236                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=0;
    1237                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=0;
    1238                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=0;
    1239         }
    1240 
    1241         for (i=0;i<numgrids;i++){ //last column not for the bubble function
    1242                 *(B_prime+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
    1243                 *(B_prime+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
    1244                 *(B_prime+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
    1245                 *(B_prime+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
    1246                 *(B_prime+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
    1247                 *(B_prime+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
    1248                 *(B_prime+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=0;
    1249                 *(B_prime+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=l1l6[i];
    1250         }
    1251 
    1252 }
    1253 /*}}}*/
    1254 /*FUNCTION PentaRef::GetBArtdiff {{{1*/
    1255 void PentaRef::GetBArtdiff(double* B_artdiff, double* xyz_list, GaussPenta* gauss){
    1256         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    1257          * For grid i, Bi' can be expressed in the actual coordinate system
    1258          * by:
    1259          *       Bi_artdiff=[ dh/dx ]
    1260          *                 [ dh/dy ]
    1261          * where h is the interpolation function for grid i.
    1262          *
    1263          * We assume B has been allocated already, of size: 2x(DOFPERGRID*numgrids)
    1264          */
    1265 
    1266         int i;
    1267         const int calculationdof=3;
    1268         const int numgrids=6;
    1269         int DOFPERGRID=1;
    1270 
    1271         /*Same thing in the actual coordinate system: */
    1272         double dh1dh6[calculationdof][numgrids];
    1273 
    1274         /*Get dh1dh6 in actual coordinates system : */
    1275         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
    1276 
    1277         /*Build B': */
    1278         for (i=0;i<numgrids;i++){
    1279                 *(B_artdiff+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
    1280                 *(B_artdiff+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
    1281         }
    1282 }
    1283 /*}}}*/
    1284 /*FUNCTION PentaRef::GetBAdvec{{{1*/
    1285 void PentaRef::GetBAdvec(double* B_advec, double* xyz_list, GaussPenta* gauss){
    1286         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    1287          * For grid i, Bi' can be expressed in the actual coordinate system
    1288          * by:
    1289          *       Bi_advec =[ h ]
    1290          *                 [ h ]
    1291          *                 [ h ]
    1292          * where h is the interpolation function for grid i.
    1293          *
    1294          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    1295          */
    1296 
    1297         int i;
    1298         int calculationdof=3;
    1299         int numgrids=6;
    1300         int DOFPERGRID=1;
    1301 
    1302         /*Same thing in the actual coordinate system: */
    1303         double l1l6[6];
    1304 
    1305         /*Get dh1dh2dh3 in actual coordinates system : */
    1306         GetNodalFunctionsP1(l1l6, gauss);
    1307 
    1308         /*Build B': */
    1309         for (i=0;i<numgrids;i++){
    1310                 *(B_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=l1l6[i];
    1311                 *(B_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=l1l6[i];
    1312                 *(B_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=l1l6[i];
    1313         }
    1314 }
    1315 /*}}}*/
    1316 /*FUNCTION PentaRef::GetBConduct{{{1*/
    1317 void PentaRef::GetBConduct(double* B_conduct, double* xyz_list, GaussPenta* gauss){
    1318         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    1319          * For grid i, Bi' can be expressed in the actual coordinate system
    1320          * by:
    1321          *       Bi_conduct=[ dh/dx ]
    1322          *                  [ dh/dy ]
    1323          *                  [ dh/dz ]
    1324          * where h is the interpolation function for grid i.
    1325          *
    1326          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    1327          */
    1328 
    1329         int i;
    1330         const int calculationdof=3;
    1331         const int numgrids=6;
    1332         int DOFPERGRID=1;
    1333 
    1334         /*Same thing in the actual coordinate system: */
    1335         double dh1dh6[calculationdof][numgrids];
    1336 
    1337         /*Get dh1dh2dh3 in actual coordinates system : */
    1338         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
    1339 
    1340         /*Build B': */
    1341         for (i=0;i<numgrids;i++){
    1342                 *(B_conduct+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
    1343                 *(B_conduct+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
    1344                 *(B_conduct+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
    1345         }
    1346 }
    1347 /*}}}*/
    1348 /*FUNCTION PentaRef::GetBVert{{{1*/
    1349 void PentaRef::GetBVert(double* B, double* xyz_list, GaussPenta* gauss){
    1350         /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
    1351                 where hi is the interpolation function for grid i.*/
    1352 
    1353         int i;
    1354         const int NDOF3=3;
    1355         const int numgrids=6;
    1356         double dh1dh6[NDOF3][numgrids];
    1357 
    1358         /*Get dh1dh6 in actual coordinate system: */
    1359         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
    1360 
    1361         /*Build B: */
    1362         for (i=0;i<numgrids;i++){
    1363                 B[i]=dh1dh6[2][i]; 
    1364         }
    1365 
    1366 }
    1367 /*}}}*/
    1368 /*FUNCTION PentaRef::GetBprimeAdvec{{{1*/
    1369 void PentaRef::GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss){
    1370         /*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*DOFPERGRID.
    1371          * For grid i, Bi' can be expressed in the actual coordinate system
    1372          * by:
    1373          *       Biprime_advec=[ dh/dx ]
    1374          *                     [ dh/dy ]
    1375          *                     [ dh/dz ]
    1376          * where h is the interpolation function for grid i.
    1377          *
    1378          * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    1379          */
    1380 
    1381         int i;
    1382         const int calculationdof=3;
    1383         const int numgrids=6;
    1384         int DOFPERGRID=1;
    1385 
    1386         /*Same thing in the actual coordinate system: */
    1387         double dh1dh6[calculationdof][numgrids];
    1388 
    1389         /*Get dh1dh2dh3 in actual coordinates system : */
    1390         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
    1391 
    1392         /*Build B': */
    1393         for (i=0;i<numgrids;i++){
    1394                 *(Bprime_advec+DOFPERGRID*numgrids*0+DOFPERGRID*i)=dh1dh6[0][i];
    1395                 *(Bprime_advec+DOFPERGRID*numgrids*1+DOFPERGRID*i)=dh1dh6[1][i];
    1396                 *(Bprime_advec+DOFPERGRID*numgrids*2+DOFPERGRID*i)=dh1dh6[2][i];
    1397         }
    1398 }
    1399 /*}}}*/
    1400 /*FUNCTION PentaRef::GetBprimeVert{{{1*/
    1401 void PentaRef::GetBprimeVert(double* B, double* xyz_list, GaussPenta* gauss){
    1402         /* Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for grid i*/
    1403 
    1404         GetNodalFunctionsP1(B, gauss);
    1405 
    1406 }
    1407 /*}}}*/
    1408 /*FUNCTION PentaRef::GetLStokes {{{1*/
    1409 void PentaRef::GetLStokes(double* LStokes, GaussPenta* gauss){
    1410         /*
    1411          * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    1412          * For grid i, Li can be expressed in the actual coordinate system
    1413          * by:
    1414          *       Li=[ h    0    0   0]
    1415          *                    [ 0    h    0   0]
    1416          *                    [ 0    0    h   0]
    1417          *                    [ 0    0    h   0]
    1418          *                    [ h    0    0   0]
    1419          *                    [ 0    h    0   0]
    1420          *                    [ h    0    0   0]
    1421          *                    [ 0    h    0   0]
    1422          *                    [ 0    0    h   0]
    1423          *                    [ 0    0    h   0]
    1424          *                    [ 0    0    h   0]
    1425          *                    [ h    0    0   0]
    1426          *                    [ 0    h    0   0]
    1427          *                    [ 0    0    h   0]
    1428          * where h is the interpolation function for grid i.
    1429          */
    1430 
    1431         int i;
    1432         const int numgrids2d=3;
    1433         int num_dof=4;
    1434 
    1435         double l1l2l3[numgrids2d];
    1436 
    1437 
    1438         /*Get l1l2l3 in actual coordinate system: */
    1439         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    1440         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    1441         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    1442 
    1443         /*Build LStokes: */
    1444         for (i=0;i<3;i++){
    1445                 *(LStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
    1446                 *(LStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
    1447                 *(LStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
    1448                 *(LStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
    1449                 *(LStokes+num_dof*numgrids2d*1+num_dof*i)=0;
    1450                 *(LStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
    1451                 *(LStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
    1452                 *(LStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
    1453                 *(LStokes+num_dof*numgrids2d*2+num_dof*i)=0;
    1454                 *(LStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
    1455                 *(LStokes+num_dof*numgrids2d*2+num_dof*i+2)=l1l2l3[i];
    1456                 *(LStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
    1457                 *(LStokes+num_dof*numgrids2d*3+num_dof*i)=0;
    1458                 *(LStokes+num_dof*numgrids2d*3+num_dof*i+1)=0;
    1459                 *(LStokes+num_dof*numgrids2d*3+num_dof*i+2)=l1l2l3[i];
    1460                 *(LStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
    1461                 *(LStokes+num_dof*numgrids2d*4+num_dof*i)=l1l2l3[i];
    1462                 *(LStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
    1463                 *(LStokes+num_dof*numgrids2d*4+num_dof*i+2)=0;
    1464                 *(LStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
    1465                 *(LStokes+num_dof*numgrids2d*5+num_dof*i)=0;
    1466                 *(LStokes+num_dof*numgrids2d*5+num_dof*i+1)=l1l2l3[i];
    1467                 *(LStokes+num_dof*numgrids2d*5+num_dof*i+2)=0;
    1468                 *(LStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
    1469                 *(LStokes+num_dof*numgrids2d*6+num_dof*i)=l1l2l3[i];
    1470                 *(LStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
    1471                 *(LStokes+num_dof*numgrids2d*6+num_dof*i+2)=0;
    1472                 *(LStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
    1473                 *(LStokes+num_dof*numgrids2d*7+num_dof*i)=0;
    1474                 *(LStokes+num_dof*numgrids2d*7+num_dof*i+1)=l1l2l3[i];
    1475                 *(LStokes+num_dof*numgrids2d*7+num_dof*i+2)=0;
    1476                 *(LStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
    1477                 *(LStokes+num_dof*numgrids2d*8+num_dof*i)=0;
    1478                 *(LStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
    1479                 *(LStokes+num_dof*numgrids2d*8+num_dof*i+2)=l1l2l3[i];
    1480                 *(LStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
    1481                 *(LStokes+num_dof*numgrids2d*9+num_dof*i)=0;
    1482                 *(LStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
    1483                 *(LStokes+num_dof*numgrids2d*9+num_dof*i+2)=l1l2l3[i];
    1484                 *(LStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
    1485                 *(LStokes+num_dof*numgrids2d*10+num_dof*i)=0;
    1486                 *(LStokes+num_dof*numgrids2d*10+num_dof*i+1)=0;
    1487                 *(LStokes+num_dof*numgrids2d*10+num_dof*i+2)=l1l2l3[i];
    1488                 *(LStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
    1489                 *(LStokes+num_dof*numgrids2d*11+num_dof*i)=l1l2l3[i];
    1490                 *(LStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
    1491                 *(LStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
    1492                 *(LStokes+num_dof*numgrids2d*11+num_dof*i+3)=0;
    1493                 *(LStokes+num_dof*numgrids2d*12+num_dof*i)=0;
    1494                 *(LStokes+num_dof*numgrids2d*12+num_dof*i+1)=l1l2l3[i];
    1495                 *(LStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
    1496                 *(LStokes+num_dof*numgrids2d*12+num_dof*i+3)=0;
    1497                 *(LStokes+num_dof*numgrids2d*13+num_dof*i)=0;
    1498                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
    1499                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
    1500                 *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
    1501 
    1502         }
    1503 }
    1504 /*}}}*/
    1505 /*FUNCTION PentaRef::GetLprimeStokes {{{1*/
    1506 void PentaRef::GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss){
    1507 
    1508         /*
    1509          * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    1510          * For grid i, Lpi can be expressed in the actual coordinate system
    1511          * by:
    1512          *       Lpi=[ h    0    0   0]
    1513          *                     [ 0    h    0   0]
    1514          *                     [ h    0    0   0]
    1515          *                     [ 0    h    0   0]
    1516          *                     [ 0    0    h   0]
    1517          *                     [ 0    0    h   0]
    1518          *                     [ 0    0  dh/dz 0]
    1519          *                     [ 0    0  dh/dz 0]
    1520          *                     [ 0    0  dh/dz 0]
    1521          *                     [dh/dz 0  dh/dx 0]
    1522          *                     [ 0 dh/dz dh/dy 0]
    1523          *           [ 0    0    0   h]
    1524          *           [ 0    0    0   h]
    1525          *           [ 0    0    0   h]
    1526          * where h is the interpolation function for grid i.
    1527          */
    1528         int i;
    1529         const int numgrids2d=3;
    1530         int num_dof=4;
    1531 
    1532         double l1l2l3[numgrids2d];
    1533         double dh1dh6[3][6];
    1534 
    1535         /*Get l1l2l3 in actual coordinate system: */
    1536         l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
    1537         l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
    1538         l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
    1539 
    1540         GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
    1541 
    1542         /*Build LprimeStokes: */
    1543         for (i=0;i<3;i++){
    1544                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
    1545                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+1)=0;
    1546                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+2)=0;
    1547                 *(LprimeStokes+num_dof*numgrids2d*0+num_dof*i+3)=0;
    1548                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i)=0;
    1549                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+1)=l1l2l3[i];
    1550                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+2)=0;
    1551                 *(LprimeStokes+num_dof*numgrids2d*1+num_dof*i+3)=0;
    1552                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i)=l1l2l3[i];
    1553                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+1)=0;
    1554                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+2)=0;
    1555                 *(LprimeStokes+num_dof*numgrids2d*2+num_dof*i+3)=0;
    1556                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i)=0;
    1557                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+1)=l1l2l3[i];
    1558                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+2)=0;
    1559                 *(LprimeStokes+num_dof*numgrids2d*3+num_dof*i+3)=0;
    1560                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i)=0;
    1561                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+1)=0;
    1562                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+2)=l1l2l3[i];
    1563                 *(LprimeStokes+num_dof*numgrids2d*4+num_dof*i+3)=0;
    1564                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i)=0;
    1565                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+1)=0;
    1566                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+2)=l1l2l3[i];
    1567                 *(LprimeStokes+num_dof*numgrids2d*5+num_dof*i+3)=0;
    1568                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i)=0;
    1569                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+1)=0;
    1570                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+2)=dh1dh6[2][i];
    1571                 *(LprimeStokes+num_dof*numgrids2d*6+num_dof*i+3)=0;
    1572                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i)=0;
    1573                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+1)=0;
    1574                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+2)=dh1dh6[2][i];
    1575                 *(LprimeStokes+num_dof*numgrids2d*7+num_dof*i+3)=0;
    1576                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i)=0;
    1577                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+1)=0;
    1578                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+2)=dh1dh6[2][i];
    1579                 *(LprimeStokes+num_dof*numgrids2d*8+num_dof*i+3)=0;
    1580                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i)=dh1dh6[2][i];
    1581                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+1)=0;
    1582                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+2)=dh1dh6[0][i];
    1583                 *(LprimeStokes+num_dof*numgrids2d*9+num_dof*i+3)=0;
    1584                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i)=0;
    1585                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+1)=dh1dh6[2][i];
    1586                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+2)=dh1dh6[1][i];
    1587                 *(LprimeStokes+num_dof*numgrids2d*10+num_dof*i+3)=0;
    1588                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i)=0;
    1589                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+1)=0;
    1590                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+2)=0;
    1591                 *(LprimeStokes+num_dof*numgrids2d*11+num_dof*i+3)=l1l2l3[i];
    1592                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i)=0;
    1593                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+1)=0;
    1594                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+2)=0;
    1595                 *(LprimeStokes+num_dof*numgrids2d*12+num_dof*i+3)=l1l2l3[i];
    1596                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i)=0;
    1597                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+1)=0;
    1598                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
    1599                 *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
    1600 
    1601         }
    1602 }
    1603 /*}}}*/
    1604 /*FUNCTION PentaRef::GetJacobian {{{1*/
    1605675void PentaRef::GetJacobian(double* J, double* xyz_list,GaussPenta* gauss){
    1606676

  • issm/trunk/src/c/objects/Elements/PentaRef.h

    r5726 r5743  
    2323
    2424                /*Numerics*/
    25                 void GetNodalFunctionsP1(double* l1l6, double* gauss);
    26                 void GetNodalFunctionsMINI(double* l1l7, double* gauss);
    27                 void GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss);
    28                 void GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss);
    29                 void GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss);
    30                 void GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss);
    31                 void GetJacobian(double* J, double* xyz_list,double* gauss);
    32                 void GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss);
    33                 void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
    34                 void GetBMacAyealPattyn(double* B, double* xyz_list, double* gauss);
    35                 void GetBPattyn(double* B, double* xyz_list, double* gauss);
    36                 void GetBprimePattyn(double* B, double* xyz_list, double* gauss);
    37                 void GetBStokes(double* B, double* xyz_list, double* gauss);
    38                 void GetBprimeStokes(double* B_prime, double* xyz_list, double* gauss);
    39                 void GetBprimeVert(double* B, double* xyz_list, double* gauss);
    40                 void GetBAdvec(double* B_advec, double* xyz_list, double* gauss);
    41                 void GetBArtdiff(double* B_artdiff, double* xyz_list, double* gauss);
    42                 void GetBConduct(double* B_conduct, double* xyz_list, double* gauss);
    43                 void GetBVert(double* B, double* xyz_list, double* gauss);
    44                 void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, double* gauss);
    45                 void GetLStokes(double* LStokes, double* gauss);
    46                 void GetLprimeStokes(double* LprimeStokes, double* xyz_list, double* gauss);
    47                 void GetParameterValue(double* pvalue,double* plist,double* gauss);
    48                 void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    49                 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, double* gauss);
    50                 void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    51 
    5225                void GetNodalFunctionsP1(double* l1l6, GaussPenta* gauss);
    5326                void GetNodalFunctionsMINI(double* l1l7, GaussPenta* gauss);
    … …  
    7548                void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);
    7649                void GetParameterValue(double* pvalue,double* plist, GaussPenta* gauss);
    77                 //void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
     50                void GetParameterValue(double* pvalue,double* plist,GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    7851                void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussPenta* gauss);
    79                 //void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
     52                void GetParameterDerivativeValue(double* pvalues, double* plist,double* xyz_list, GaussTria* gauss){ISSMERROR("only PentaGauss are supported");};
    8053
    8154};

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r5741 r5743  
    807807}
    808808/*}}}*/
    809 /*FUNCTION Tria::GetNodes {{{1*/
    810 void  Tria::GetNodes(void** vpnodes){
    811         int i;
    812         Node** pnodes=NULL;
    813 
    814         /*recover nodes: */
    815         pnodes=(Node**)vpnodes;
    816 
    817         for(i=0;i<3;i++){
    818                 pnodes[i]=nodes[i];
    819         }
    820 }
    821 /*}}}*/
    822809/*FUNCTION Tria::GetNodeIndex {{{1*/
    823810int Tria::GetNodeIndex(Node* node){
    … …  
    55025489}
    55035490/*}}}*/
    5504 /*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype) {{{1*/
    5505 void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype){
    5506 
    5507         /*Output*/
    5508         double value;
    5509 
    5510         /*Intermediaries*/
    5511         const int numnodes=3;
    5512         int       grid1=-1,grid2=-1;
    5513         int       grid3;
    5514         int       i;
    5515         double    gauss_tria[numnodes];
    5516 
    5517         /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */
    5518         ISSMASSERT(nodes);
    5519         for(i=0;i<numnodes;i++){
    5520                 if (node1==nodes[i]) grid1=i;
    5521                 if (node2==nodes[i]) grid2=i;
    5522         }
    5523 
    5524         /*Reverse grid1 and 2 if necessary*/
    5525         if (grid1>grid2){
    5526 
    5527                 /*Reverse grid1 and grid2*/
    5528                 grid3=grid1; grid1=grid2; grid2=grid3;
    5529 
    5530                 /*Change segment gauss point (between -1 and +1)*/
    5531                 gauss_seg=-gauss_seg;
    5532         }
    5533 
    5534         /*Build Triangle Gauss point*/
    5535         if (grid1==0 && grid2==1){
    5536                 /*gauss_seg is between 0 and 1*/
    5537                 gauss_tria[0]=0.5*(1.-gauss_seg);
    5538                 gauss_tria[1]=1.-0.5*(1.-gauss_seg);
    5539                 gauss_tria[2]=0.;
    5540         }
    5541         else if (grid1==0 && grid2==2){
    5542                 /*gauss_seg is between 0 and 2*/
    5543                 gauss_tria[0]=0.5*(1.-gauss_seg);
    5544                 gauss_tria[1]=0.;
    5545                 gauss_tria[2]=1.-0.5*(1.-gauss_seg);
    5546         }
    5547         else if (grid1==1 && grid2==2){
    5548                 /*gauss_seg is between 1 and 2*/
    5549                 gauss_tria[0]=0.;
    5550                 gauss_tria[1]=0.5*(1.-gauss_seg);
    5551                 gauss_tria[2]=1.-0.5*(1.-gauss_seg);
    5552         }
    5553         else{
    5554                 ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes");
    5555         }
    5556 
    5557         /*OK, now we can call input method*/
    5558         this->inputs->GetParameterValue(&value, &gauss_tria[0],enumtype);
    5559 
    5560         /*Assign output pointers:*/
    5561         *pvalue=value;
    5562 }
    5563 /*}}}*/
    5564 /*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in) {{{1*/
    5565 void Tria::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in){
    5566 
    5567         /*Output*/
    5568         double value;
    5569 
    5570         /*Intermediaries*/
    5571         const int numnodes=3;
    5572         int       grid1=-1,grid2=-1;
    5573         int       grid3;
    5574         int       i;
    5575         double    gauss_tria[numnodes];
    5576 
    5577         /*go through 3 nodes (all nodes for tria) and identify 1st and 2nd nodes: */
    5578         ISSMASSERT(nodes);
    5579         for(i=0;i<numnodes;i++){
    5580                 if (node1==nodes[i]) grid1=i;
    5581                 if (node2==nodes[i]) grid2=i;
    5582         }
    5583 
    5584         /*Reverse grid1 and 2 if necessary*/
    5585         if (grid1>grid2){
    5586 
    5587                 /*Reverse grid1 and grid2*/
    5588                 grid3=grid1; grid1=grid2; grid2=grid3;
    5589 
    5590                 /*Change segment gauss point (between -1 and +1)*/
    5591                 gauss_seg=-gauss_seg;
    5592         }
    5593 
    5594         /*Build Triangle Gauss point*/
    5595         if (grid1==0 && grid2==1){
    5596                 /*gauss_seg is between 0 and 1*/
    5597                 gauss_tria[0]=0.5*(1.-gauss_seg);
    5598                 gauss_tria[1]=1.-0.5*(1.-gauss_seg);
    5599                 gauss_tria[2]=0.;
    5600         }
    5601         else if (grid1==0 && grid2==2){
    5602                 /*gauss_seg is between 0 and 2*/
    5603                 gauss_tria[0]=0.5*(1.-gauss_seg);
    5604                 gauss_tria[1]=0.;
    5605                 gauss_tria[2]=1.-0.5*(1.-gauss_seg);
    5606         }
    5607         else if (grid1==1 && grid2==2){
    5608                 /*gauss_seg is between 1 and 2*/
    5609                 gauss_tria[0]=0.;
    5610                 gauss_tria[1]=0.5*(1.-gauss_seg);
    5611                 gauss_tria[2]=1.-0.5*(1.-gauss_seg);
    5612         }
    5613         else{
    5614                 ISSMERROR("The 2 nodes provided are either the same or did not match current Tria nodes");
    5615         }
    5616 
    5617         /*OK, now we can call input method*/
    5618         input_in->GetParameterValue(&value, &gauss_tria[0]);
    5619 
    5620         /*Assign output pointers:*/
    5621         *pvalue=value;
    5622 }
    5623 /*}}}*/
    56245491/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    56255492void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){

  • issm/trunk/src/c/objects/Elements/Tria.h

    r5734 r5743  
    7474                void   CreateKMatrix(Mat Kgg);
    7575                void   CreatePVector(Vec pg);
    76                 void   GetNodes(void** nodes);
    7776                int    GetNodeIndex(Node* node);
    7877                bool   GetOnBed();
    … …  
    153152                void    GetParameterListOnVertices(double* pvalue,int enumtype);
    154153                void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
    155                 void      GetParameterValue(double* pvalue,Node* node,int enumtype);
    156                 void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
    157                 void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);
     154                void    GetParameterValue(double* pvalue,Node* node,int enumtype);
    158155                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    159156                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);

  • issm/trunk/src/c/objects/Elements/TriaRef.cpp

    r5738 r5743  
    5252
    5353/*Reference Element numerics*/
    54 /*FUNCTION TriaRef::GetBMacAyeal {{{1*/
    55 void TriaRef::GetBMacAyeal(double* B, double* xyz_list, double* gauss){
    56         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    57          * For grid i, Bi can be expressed in the actual coordinate system
    58          * by:
    59          *       Bi=[ dh/dx           0    ]
    60          *          [   0           dh/dy  ]
    61          *          [ 1/2*dh/dy  1/2*dh/dx ]
    62          * where h is the interpolation function for grid i.
    63          *
    64          * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
    65          */
    66 
    67         int i;
    68         const int NDOF2=2;
    69         const int numgrids=3;
    70 
    71         double dh1dh3[NDOF2][numgrids];
    72 
    73 
    74         /*Get dh1dh2dh3 in actual coordinate system: */
    75         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
    76 
    77         /*Build B: */
    78         for (i=0;i<numgrids;i++){
    79                 *(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
    80                 *(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
    81                 *(B+NDOF2*numgrids*1+NDOF2*i)=0;
    82                 *(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
    83                 *(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i];
    84                 *(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i];
    85         }
    86 }
    87 /*}}}*/
    8854/*FUNCTION TriaRef::GetBMacAyeal {{{1*/
    8955void TriaRef::GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss){
    … …  
    161127/*}}}*/
    162128/*FUNCTION TriaRef::GetBPrognostic{{{1*/
    163 void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, double* gauss){
    164         /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    165          * For grid i, Bi can be expressed in the actual coordinate system
    166          * by:
    167          *       Bi=[ h ]
    168          *          [ h ]
    169          * where h is the interpolation function for grid i.
    170          *
    171          * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
    172          */
    173 
    174         int i;
    175         const int NDOF1=1;
    176         const int numgrids=3;
    177 
    178         double l1l2l3[numgrids];
    179 
    180 
    181         /*Get dh1dh2dh3 in actual coordinate system: */
    182         GetNodalFunctions(&l1l2l3[0],gauss);
    183 
    184         /*Build B_prog: */
    185         for (i=0;i<numgrids;i++){
    186                 *(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
    187                 *(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
    188         }
    189 }
    190 /*}}}*/
    191 /*FUNCTION TriaRef::GetBPrognostic{{{1*/
    192129void TriaRef::GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss){
    193130        /*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
    … …  
    219156/*}}}*/
    220157/*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/
    221 void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss){
     158void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss){
    222159
    223160        /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    … …  
    253190}
    254191/*}}}*/
    255 /*FUNCTION TriaRef::GetBprimeMacAyeal {{{1*/
    256 void TriaRef::GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss){
    257 
    258         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    259          * For grid i, Bi' can be expressed in the actual coordinate system
    260          * by:
    261          *       Bi_prime=[ 2*dh/dx    dh/dy ]
    262          *                [   dh/dx  2*dh/dy ]
    263          *                [   dh/dy    dh/dx ]
    264          * where h is the interpolation function for grid i.
    265          *
    266          * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
    267          */
    268 
    269         int i;
    270         const int NDOF2=2;
    271         const int numgrids=3;
    272 
    273         /*Same thing in the actual coordinate system: */
    274         double dh1dh3[NDOF2][numgrids];
    275 
    276         /*Get dh1dh2dh3 in actual coordinates system : */
    277         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
    278 
    279         /*Build B': */
    280         for (i=0;i<numgrids;i++){
    281                 *(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i];
    282                 *(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i];
    283                 *(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i];
    284                 *(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i];
    285                 *(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i];
    286                 *(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i];
    287         }
    288 }
    289 /*}}}*/
    290 /*FUNCTION TriaRef::GetBprimePrognostic{{{1*/
    291 void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss){
    292         /*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
    293          * For grid i, Bi' can be expressed in the actual coordinate system
    294          * by:
    295          *       Bi_prime=[ dh/dx ]
    296          *                [ dh/dy ]
    297          * where h is the interpolation function for grid i.
    298          *
    299          * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
    300          */
    301 
    302         int i;
    303         const int NDOF1=1;
    304         const int NDOF2=2;
    305         const int numgrids=3;
    306 
    307         /*Same thing in the actual coordinate system: */
    308         double dh1dh3[NDOF2][numgrids];
    309 
    310         /*Get dh1dh2dh3 in actual coordinates system : */
    311         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
    312 
    313         /*Build B': */
    314         for (i=0;i<numgrids;i++){
    315                 *(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i];
    316                 *(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i];
    317         }
    318 }
    319 /*}}}*/
    320192/*FUNCTION TriaRef::GetBprimePrognostic{{{1*/
    321193void TriaRef::GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss){
    … …  
    348220}
    349221/*}}}*/
    350 /*FUNCTION TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof) {{{1*/
    351 void TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof){
    352         /*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    353          * For grid i, Li can be expressed in the actual coordinate system
    354          * by:
    355          *       numdof=1:
    356          *                 Li=h;
    357          *       numdof=2:
    358          *                 Li=[ h   0 ]
    359          *                    [ 0   h ]
    360          * where h is the interpolation function for grid i.
    361          *
    362          * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
    363          */
    364 
    365         int i;
    366         const int NDOF2=2;
    367         const int numgrids=3;
    368 
    369         double l1l2l3[3];
    370 
    371 
    372         /*Get l1l2l3 in actual coordinate system: */
    373         GetNodalFunctions(l1l2l3, gauss);
    374 
    375         /*Build L: */
    376         if(numdof==1){
    377                 for (i=0;i<numgrids;i++){
    378                         L[i]=l1l2l3[i];
    379                 }
    380         }
    381         else{
    382                 for (i=0;i<numgrids;i++){
    383                         *(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
    384                         *(L+numdof*numgrids*0+numdof*i+1)=0;
    385                         *(L+numdof*numgrids*1+numdof*i)=0;
    386                         *(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
    387                 }
    388         }
    389 }
    390 /*}}}*/
    391 /*FUNCTION TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof) {{{1*/
     222/*FUNCTION TriaRef::GetL{{{1*/
    392223void TriaRef::GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof){
    393224        /*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    … …  
    428259}
    429260/*}}}*/
    430 /*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss) {{{1*/
    431 void TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss){
     261/*FUNCTION TriaRef::GetJacobian{{{1*/
     262void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){
    432263        /*The Jacobian is constant over the element, discard the gaussian points.
    433264         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    … …  
    451282}
    452283/*}}}*/
    453 /*FUNCTION TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/
    454 void TriaRef::GetJacobian(double* J, double* xyz_list,GaussTria* gauss){
    455         /*The Jacobian is constant over the element, discard the gaussian points.
    456          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    457 
    458         const int NDOF2=2;
    459         const int numgrids=3;
    460         double x1,y1,x2,y2,x3,y3;
    461 
    462         x1=*(xyz_list+numgrids*0+0);
    463         y1=*(xyz_list+numgrids*0+1);
    464         x2=*(xyz_list+numgrids*1+0);
    465         y2=*(xyz_list+numgrids*1+1);
    466         x3=*(xyz_list+numgrids*2+0);
    467         y3=*(xyz_list+numgrids*2+1);
    468 
    469 
    470         *(J+NDOF2*0+0)=0.5*(x2-x1);
    471         *(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
    472         *(J+NDOF2*0+1)=0.5*(y2-y1);
    473         *(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
    474 }
    475 /*}}}*/
    476 /*FUNCTION TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss) {{{1*/
     284/*FUNCTION TriaRef::GetSegmentJacobian{{{1*/
    477285void TriaRef::GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss){
    478286        /*The Jacobian is constant over the element, discard the gaussian points.
    … …  
    489297}
    490298/*}}}*/
    491 /*FUNCTION TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/
     299/*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{1*/
    492300void TriaRef::GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss){
    493301        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    … …  
    499307}
    500308/*}}}*/
    501 /*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss) {{{1*/
    502 void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss){
    503 
    504         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    505          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    506         double J[2][2];
    507 
    508         /*Get Jacobian*/
    509         GetJacobian(&J[0][0],xyz_list,gauss);
    510 
    511         /*Get Determinant*/
    512         Matrix2x2Determinant(Jdet,&J[0][0]);
    513         if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
    514 
    515 }
    516 /*}}}*/
    517 /*FUNCTION TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss) {{{1*/
     309/*FUNCTION TriaRef::GetJacobianDeterminant2d{{{1*/
    518310void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss){
    519311        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    … …  
    531323/*}}}*/
    532324/*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/
    533 void TriaRef::GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,double* gauss){
     325void TriaRef::GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,GaussTria* gauss){
    534326        /*The Jacobian determinant is constant over the element, discard the gaussian points.
    535327         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    … …  
    552344}
    553345/*}}}*/
    554 /*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/
    555 void TriaRef::GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,GaussTria* gauss){
    556         /*The Jacobian determinant is constant over the element, discard the gaussian points.
    557          * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    558 
    559         double x1,x2,x3,y1,y2,y3,z1,z2,z3;
    560 
    561         x1=*(xyz_list+3*0+0);
    562         y1=*(xyz_list+3*0+1);
    563         z1=*(xyz_list+3*0+2);
    564         x2=*(xyz_list+3*1+0);
    565         y2=*(xyz_list+3*1+1);
    566         z2=*(xyz_list+3*1+2);
    567         x3=*(xyz_list+3*2+0);
    568         y3=*(xyz_list+3*2+1);
    569         z3=*(xyz_list+3*2+2);
    570 
    571         *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
    572         if(*Jdet<0) ISSMERROR("negative jacobian determinant!");
    573 
    574 }
    575 /*}}}*/
    576 /*FUNCTION TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss) {{{1*/
    577 void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
     346/*FUNCTION TriaRef::GetJacobianInvert{{{1*/
     347void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss){
    578348
    579349        /*Jacobian*/
    … …  
    588358}
    589359/*}}}*/
    590 /*FUNCTION TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss) {{{1*/
    591 void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss){
    592 
    593         /*Jacobian*/
    594         double J[2][2];
    595 
    596         /*Call Jacobian routine to get the jacobian:*/
    597         GetJacobian(&J[0][0], xyz_list, gauss);
    598 
    599         /*Invert Jacobian matrix: */
    600         Matrix2x2Invert(Jinv,&J[0][0]);
    601 
    602 }
    603 /*}}}*/
    604 /*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss) {{{1*/
    605 void TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss){
    606         /*This routine returns the values of the nodal functions  at the gaussian point.*/
    607 
    608         l1l2l3[0]=gauss[0];
    609         l1l2l3[1]=gauss[1];
    610         l1l2l3[2]=gauss[2];
    611 
    612 }
    613 /*}}}*/
    614 /*FUNCTION TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss) {{{1*/
     360/*FUNCTION TriaRef::GetNodalFunctions{{{1*/
    615361void TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss){
    616362        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    622368}
    623369/*}}}*/
    624 /*FUNCTION TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss) {{{1*/
     370/*FUNCTION TriaRef::GetSegmentNodalFunctions{{{1*/
    625371void TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss,int index1,int index2){
    626372        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    … …  
    636382}
    637383/*}}}*/
    638 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss) {{{1*/
    639 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){
     384/*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{1*/
     385void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss){
    640386
    641387        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    … …  
    665411}
    666412/*}}}*/
    667 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss) {{{1*/
    668 void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss){
    669 
    670         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    671          * actual coordinate system): */
    672         int       i;
    673         const int NDOF2    = 2;
    674         const int numgrids = 3;
    675         double    dh1dh3_ref[NDOF2][numgrids];
    676         double    Jinv[NDOF2][NDOF2];
    677 
    678         /*Get derivative values with respect to parametric coordinate system: */
    679         GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss);
    680 
    681         /*Get Jacobian invert: */
    682         GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
    683 
    684         /*Build dh1dh3:
    685          *
    686          * [dhi/dx]= Jinv*[dhi/dr]
    687          * [dhi/dy]       [dhi/ds]
    688          */
    689         for (i=0;i<numgrids;i++){
    690                 dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
    691                 dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
    692         }
    693 
    694 }
    695 /*}}}*/
    696 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss) {{{1*/
    697 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss){
     413/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{1*/
     414void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){
    698415
    699416        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    … …  
    717434}
    718435/*}}}*/
    719 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss) {{{1*/
    720 void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss){
    721 
    722         /*This routine returns the values of the nodal functions derivatives  (with respect to the
    723          * natural coordinate system) at the gaussian point. */
    724 
    725         const int NDOF2=2;
    726         const int numgrids=3;
    727 
    728         /*First nodal function: */
    729         *(dl1dl3+numgrids*0+0)=-0.5;
    730         *(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
    731 
    732         /*Second nodal function: */
    733         *(dl1dl3+numgrids*0+1)=0.5;
    734         *(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
    735 
    736         /*Third nodal function: */
    737         *(dl1dl3+numgrids*0+2)=0;
    738         *(dl1dl3+numgrids*1+2)=1.0/SQRT3;
    739 
    740 }
    741 /*}}}*/
    742 /*FUNCTION TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss) {{{1*/
    743 void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss){
     436/*FUNCTION TriaRef::GetParameterDerivativeValue{{{1*/
     437void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss){
    744438
    745439        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
    … …  
    763457}
    764458/*}}}*/
    765 /*FUNCTION TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss) {{{1*/
    766 void TriaRef::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, GaussTria* gauss){
    767 
    768         /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian
    769          * point specified by gauss_l1l2l3:
    770          *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    771          *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
    772          *
    773          * p is a vector of size 2x1 already allocated.
    774          */
    775 
    776         /*Nodal Derivatives*/
    777         double dh1dh3[2][3]; //nodal derivative functions in actual coordinate system.
    778 
    779         /*Get dh1dh2dh3 in actual coordinate system: */
    780         GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss);
    781 
    782         /*Assign values*/
    783         *(p+0)=plist[0]*dh1dh3[0][0]+plist[1]*dh1dh3[0][1]+plist[2]*dh1dh3[0][2];
    784         *(p+1)=plist[0]*dh1dh3[1][0]+plist[1]*dh1dh3[1][1]+plist[2]*dh1dh3[1][2];
    785 
    786 }
    787 /*}}}*/
    788 /*FUNCTION TriaRef::GetParameterValue(double* p, double* plist, double* gauss){{{1*/
    789 void TriaRef::GetParameterValue(double* p, double* plist, double* gauss){
     459/*FUNCTION TriaRef::GetParameterValue{{{1*/
     460void TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){
    790461
    791462        /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian
    … …  
    802473}
    803474/*}}}*/
    804 /*FUNCTION TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){{{1*/
    805 void TriaRef::GetParameterValue(double* p, double* plist, GaussTria* gauss){
    806 
    807         /*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian
    808          * point specifie by gauss: */
    809 
    810         /*nodal functions annd output: */
    811         double l1l2l3[3];
    812 
    813         /*Get nodal functions*/
    814         GetNodalFunctions(l1l2l3, gauss);
    815 
    816         /*Get parameter*/
    817         *p=l1l2l3[0]*plist[0]+l1l2l3[1]*plist[1]+l1l2l3[2]*plist[2];
    818 }
    819 /*}}}*/

  • issm/trunk/src/c/objects/Elements/TriaRef.h

    r5738 r5743  
    2525
    2626                /*Numerics*/
    27                 void GetBMacAyeal(double* B, double* xyz_list, double* gauss);
    2827                void GetBMacAyeal(double* B, double* xyz_list, GaussTria* gauss);
    29                 void GetBprimeMacAyeal(double* Bprime, double* xyz_list, double* gauss);
    3028                void GetBprimeMacAyeal(double* Bprime, double* xyz_list, GaussTria* gauss);
    31                 void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, double* gauss);
    3229                void GetBprimePrognostic(double* Bprime_prog, double* xyz_list, GaussTria* gauss);
    33                 void GetBPrognostic(double* B_prog, double* xyz_list, double* gauss);
    3430                void GetBPrognostic(double* B_prog, double* xyz_list, GaussTria* gauss);
    35                 void GetL(double* L, double* xyz_list,double* gauss,int numdof);
    3631                void GetL(double* L, double* xyz_list,GaussTria* gauss,int numdof);
    37                 void GetJacobian(double* J, double* xyz_list,double* gauss);
    3832                void GetJacobian(double* J, double* xyz_list,GaussTria* gauss);
    3933                void GetSegmentJacobian(double* J, double* xyz_list,GaussTria* gauss);
    4034                void GetSegmentJacobianDeterminant(double* Jdet, double* xyz_list,GaussTria* gauss);
    41                 void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss);
    4235                void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,GaussTria* gauss);
    43                 void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss);
    4436                void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,GaussTria* gauss);
    45                 void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
    4637                void GetJacobianInvert(double*  Jinv, double* xyz_list,GaussTria* gauss);
    47                 void GetNodalFunctions(double* l1l2l3, double* gauss);
    4838                void GetNodalFunctions(double* l1l2l3,GaussTria* gauss);
    4939                void GetSegmentNodalFunctions(double* l1l2l3,GaussTria* gauss, int index1,int index2);
    5040                void GetSegmentBFlux(double* B,GaussTria* gauss, int index1,int index2);
    5141                void GetSegmentBprimeFlux(double* Bprime,GaussTria* gauss, int index1,int index2);
    52                 void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);
    5342                void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, GaussTria* gauss);
    54                 void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss);
    5543                void GetNodalFunctionsDerivativesReference(double* dl1dl3,GaussTria* gauss);
    56                 void GetParameterValue(double* pp, double* plist, double* gauss);
    5744                void GetParameterValue(double* pp, double* plist, GaussTria* gauss);
    58                 void GetParameterDerivativeValue(double* pp, double* plist,double* xyz_list, double* gauss);
    5945                void GetParameterDerivativeValue(double* pp, double* plist,double* xyz_list, GaussTria* gauss);
    6046

  • issm/trunk/src/c/objects/Inputs/BoolInput.cpp

    r5660 r5743  
    166166void BoolInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");}
    167167/*}}}*/
    168 /*FUNCTION BoolInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
    169 void BoolInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
    170 /*}}}*/
    171168/*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
    172169void BoolInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
    … …  
    174171/*FUNCTION BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
    175172void BoolInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");}
    176 /*}}}*/
    177 /*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
    178 void BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
    179173/*}}}*/
    180174/*FUNCTION BoolInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/

  • issm/trunk/src/c/objects/Inputs/BoolInput.h

    r5736 r5743  
    4747                void GetParameterValue(int* pvalue);
    4848                void GetParameterValue(double* pvalue);
    49                 void GetParameterValue(double* pvalue,double* gauss);
    5049                void GetParameterValue(double* pvalue,GaussTria* gauss);
    5150                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    52                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5351                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5452                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);

  • issm/trunk/src/c/objects/Inputs/DoubleInput.cpp

    r5707 r5743  
    182182}
    183183/*}}}*/
    184 /*FUNCTION DoubleInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
    185 void DoubleInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
    186 /*}}}*/
    187184/*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
    188185void DoubleInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
    … …  
    190187/*FUNCTION DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
    191188void DoubleInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");}
    192 /*}}}*/
    193 /*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
    194 void DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
    195189/*}}}*/
    196190/*FUNCTION DoubleInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/

  • issm/trunk/src/c/objects/Inputs/DoubleInput.h

    r5736 r5743  
    4646                void GetParameterValue(int* pvalue);
    4747                void GetParameterValue(double* pvalue);
    48                 void GetParameterValue(double* pvalue,double* gauss);
    4948                void GetParameterValue(double* pvalue,GaussTria* gauss);
    5049                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    51                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5250                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5351                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);

  • issm/trunk/src/c/objects/Inputs/Input.h

    r5736 r5743  
    2525                virtual void GetParameterValue(int* pvalue)=0;
    2626                virtual void GetParameterValue(double* pvalue)=0;
    27                 virtual void GetParameterValue(double* pvalue,double* gauss)=0;
    2827                virtual void GetParameterValue(double* pvalue,GaussTria* gauss)=0;
    2928                virtual void GetParameterValue(double* pvalue,GaussPenta* gauss)=0;
    30                 virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss)=0;
    3129                virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss)=0;
    3230                virtual void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss)=0;

  • issm/trunk/src/c/objects/Inputs/IntInput.cpp

    r5660 r5743  
    167167}
    168168/*}}}*/
    169 /*FUNCTION IntInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
    170 void IntInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
    171 /*}}}*/
    172169/*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
    173170void IntInput::GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR(" not supported yet!");}
    … …  
    175172/*FUNCTION IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
    176173void IntInput::GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR(" not supported yet!");}
    177 /*}}}*/
    178 /*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
    179 void IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
    180174/*}}}*/
    181175/*FUNCTION IntInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){{{1*/

  • issm/trunk/src/c/objects/Inputs/IntInput.h

    r5736 r5743  
    4747                void GetParameterValue(int* pvalue);
    4848                void GetParameterValue(double* pvalue);
    49                 void GetParameterValue(double* pvalue,double* gauss);
    5049                void GetParameterValue(double* pvalue,GaussTria* gauss);
    5150                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    52                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5351                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5452                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp

    r5736 r5743  
    176176
    177177/*Object functions*/
    178 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
    179 void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
     178/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
     179void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){
    180180
    181181        /*Call PentaRef function*/
    182182        PentaRef::GetParameterValue(pvalue,&values[0],gauss);
    183183
    184 }
    185 /*}}}*/
    186 /*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){{{1*/
    187 void PentaVertexInput::GetParameterValue(double* pvalue,GaussPenta* gauss){
    188 
    189         /*Call PentaRef function*/
    190         PentaRef::GetParameterValue(pvalue,&values[0],gauss);
    191 
    192 }
    193 /*}}}*/
    194 /*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/
    195 void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
    196 
    197         /*Call PentaRef function*/
    198         PentaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
    199184}
    200185/*}}}*/

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.h

    r5736 r5743  
    4747                void GetParameterValue(int* pvalue){ISSMERROR("not implemented yet");};
    4848                void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");};
    49                 void GetParameterValue(double* pvalue,double* gauss);
    5049                void GetParameterValue(double* pvalue,GaussTria* gauss){ISSMERROR("not implemented yet");};
    5150                void GetParameterValue(double* pvalue,GaussPenta* gauss);
    52                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5351                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss){ISSMERROR("not implemented yet");};
    5452                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss);

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp

    r5736 r5743  
    165165
    166166/*Object functions*/
    167 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
    168 void TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){
     167/*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
     168void TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){
    169169
    170170        /*Call TriaRef function*/
    171171        TriaRef::GetParameterValue(pvalue,&values[0],gauss);
    172172
    173 }
    174 /*}}}*/
    175 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){{{1*/
    176 void TriaVertexInput::GetParameterValue(double* pvalue,GaussTria* gauss){
    177 
    178         /*Call TriaRef function*/
    179         TriaRef::GetParameterValue(pvalue,&values[0],gauss);
    180 
    181 }
    182 /*}}}*/
    183 /*FUNCTION TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){{{1*/
    184 void TriaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
    185 
    186         /*Call TriaRef function*/
    187         TriaRef::GetParameterDerivativeValue(p,&values[0],xyz_list,gauss);
    188173}
    189174/*}}}*/

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.h

    r5736 r5743  
    4747                void GetParameterValue(int* pvalue){ISSMERROR("not implemented yet");}
    4848                void GetParameterValue(double* pvalue){ISSMERROR("not implemented yet");}
    49                 void GetParameterValue(double* pvalue,double* gauss);
    5049                void GetParameterValue(double* pvalue,GaussTria* gauss);
    5150                void GetParameterValue(double* pvalue,GaussPenta* gauss){ISSMERROR("not implemented yet");};
    52                 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss);
    5351                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussTria* gauss);
    5452                void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, GaussPenta* gauss){ISSMERROR("not implemented yet");};

  • issm/trunk/src/c/objects/Loads/Friction.cpp

    r5682 r5743  
    5555}
    5656/*}}}*/
    57 /*FUNCTION Friction::GetAlpha2{{{1*/
    58 void Friction::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){
     57/*FUNCTION Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){{{1*/
     58void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
    5959
    6060        /*This routine calculates the basal friction coefficient
    61         alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
     61          alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
    6262
    6363        /*diverse: */
    … …  
    8787        r=drag_q/drag_p;
    8888        s=1./drag_p;
    89                
     89
    9090        //From bed and thickness, compute effective pressure when drag is viscous:
    9191        Neff=gravity*(rho_ice*thickness+rho_water*bed);
    92        
     92
    9393        /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    94         the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    95         for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    96         replace it by Neff=0 (ie, equival it to an ice shelf)*/
     94          the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
     95          for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
     96          replace it by Neff=0 (ie, equival it to an ice shelf)*/
    9797        if (Neff<0)Neff=0;
    9898
    … …  
    109109        }
    110110        else ISSMERROR("element_type %s not supported yet",element_type);
    111        
     111
    112112        alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
    113113
    … …  
    116116}
    117117/*}}}*/
    118 /*FUNCTION Friction::GetAlpha2{{{1*/
    119 void Friction::GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum){
     118/*FUNCTION Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){{{1*/
     119void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
    120120
    121121        /*This routine calculates the basal friction coefficient
    … …  
    177177}
    178178/*}}}*/
    179 /*FUNCTION Friction::GetAlpha2{{{1*/
    180 void Friction::GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum){
    181 
    182         /*This routine calculates the basal friction coefficient
    183           alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
    184 
    185         /*diverse: */
    186         double  r,s;
    187         double  drag_p, drag_q;
    188         double  gravity,rho_ice,rho_water;
    189         double  Neff;
    190         double  thickness,bed;
    191         double  vx,vy,vz,vmag;
    192         double  drag_coefficient;
    193         double  alpha2;
    194 
    195 
    196         /*Recover parameters: */
    197         inputs->GetParameterValue(&drag_p,DragPEnum);
    198         inputs->GetParameterValue(&drag_q,DragQEnum);
    199         inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
    200         inputs->GetParameterValue(&bed, gauss,BedEnum);
    201         inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
    202 
    203         /*Get material parameters: */
    204         gravity=matpar->GetG();
    205         rho_ice=matpar->GetRhoIce();
    206         rho_water=matpar->GetRhoWater();
    207 
    208         //compute r and q coefficients: */
    209         r=drag_q/drag_p;
    210         s=1./drag_p;
    211 
    212         //From bed and thickness, compute effective pressure when drag is viscous:
    213         Neff=gravity*(rho_ice*thickness+rho_water*bed);
    214 
    215         /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    216           the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    217           for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    218           replace it by Neff=0 (ie, equival it to an ice shelf)*/
    219         if (Neff<0)Neff=0;
    220 
    221         if(strcmp(element_type,"2d")==0){
    222                 inputs->GetParameterValue(&vx, gauss,vxenum);
    223                 inputs->GetParameterValue(&vy, gauss,vyenum);
    224                 vmag=sqrt(pow(vx,2)+pow(vy,2));
    225         }
    226         else if (strcmp(element_type,"3d")==0){
    227                 inputs->GetParameterValue(&vx, gauss,vxenum);
    228                 inputs->GetParameterValue(&vy, gauss,vyenum);
    229                 inputs->GetParameterValue(&vz, gauss,vzenum);
    230                 vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
    231         }
    232         else ISSMERROR("element_type %s not supported yet",element_type);
    233 
    234         alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
    235 
    236         /*Assign output pointers:*/
    237         *palpha2=alpha2;
    238 }
    239 /*}}}*/
    240179/*FUNCTION Friction::GetAlphaComplement {{{1*/
    241 void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){
    242        
     180void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){
     181
    243182        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
    244183         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    … …  
    274213        r=drag_q/drag_p;
    275214        s=1./drag_p;
    276                
     215
    277216        //From bed and thickness, compute effective pressure when drag is viscous:
    278217        Neff=gravity*(rho_ice*thickness+rho_water*bed);
    … …  
    288227        inputs->GetParameterValue(&vy, gauss,vyenum);
    289228        vmag=sqrt(pow(vx,2)+pow(vy,2));
    290        
     229
    291230        alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
    292231
    … …  
    296235}
    297236/*}}}*/
    298 /*FUNCTION Friction::GetAlphaComplement {{{1*/
    299 void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){
    300 
    301         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
    302          * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    303          * alpha_complement= Neff ^r * vel ^s*/
    304 
    305         /*diverse: */
    306         int     i;
    307         double  Neff;
    308         double  r,s;
    309         double  vx;
    310         double  vy;
    311         double  vmag;
    312         double  drag_p,drag_q;
    313         double  drag_coefficient;
    314         double  bed,thickness;
    315         double  gravity,rho_ice,rho_water;
    316         double  alpha_complement;
    317 
    318         /*Recover parameters: */
    319         inputs->GetParameterValue(&drag_p,DragPEnum);
    320         inputs->GetParameterValue(&drag_q,DragQEnum);
    321         inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
    322         inputs->GetParameterValue(&bed, gauss,BedEnum);
    323         inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
    324 
    325         /*Get material parameters: */
    326         gravity=matpar->GetG();
    327         rho_ice=matpar->GetRhoIce();
    328         rho_water=matpar->GetRhoWater();
    329 
    330 
    331         //compute r and q coefficients: */
    332         r=drag_q/drag_p;
    333         s=1./drag_p;
    334 
    335         //From bed and thickness, compute effective pressure when drag is viscous:
    336         Neff=gravity*(rho_ice*thickness+rho_water*bed);
    337 
    338         /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    339           the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    340           for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    341           replace it by Neff=0 (ie, equival it to an ice shelf)*/
    342         if (Neff<0)Neff=0;
    343 
    344         //We need the velocity magnitude to evaluate the basal stress:
    345         inputs->GetParameterValue(&vx, gauss,vxenum);
    346         inputs->GetParameterValue(&vy, gauss,vyenum);
    347         vmag=sqrt(pow(vx,2)+pow(vy,2));
    348 
    349         alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
    350 
    351         /*Assign output pointers:*/
    352         *palpha_complement=alpha_complement;
    353 
    354 }
    355 /*}}}*/

  • issm/trunk/src/c/objects/Loads/Friction.h

    r5682 r5743  
    2727       
    2828                void  Echo(void);
    29                 void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
    3029                void  GetAlpha2(double* palpha2, GaussTria* gauss,int vxenum,int vyenum,int vzenum);
    3130                void  GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    32                 void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
    3331                void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);
    3432

  • issm/trunk/src/c/objects/Loads/Icefront.cpp

    r5742 r5743  
    549549        /*Objects: */
    550550        double  pe_g[numdofs]={0.0};
    551         Node**  element_nodes=NULL;
    552551        int    element_type;
    553552               
    … …  
    579578       
    580579        //Identify which grids are comprised in the quad:
    581         if(element_type==PentaEnum)element_nodes=(Node**)xmalloc(6*sizeof(Node*));
    582         element->GetNodes((void**)element_nodes);
    583 
    584         grid1=-1; grid2=-1; grid3=-1; grid4=-1;
    585         for(i=0;i<6;i++){
    586                 if (nodes[0]==element_nodes[i])grid1=i;
    587                 if (nodes[1]==element_nodes[i])grid2=i;
    588                 if (nodes[2]==element_nodes[i])grid3=i;
    589                 if (nodes[3]==element_nodes[i])grid4=i;
    590         }
    591        
    592         if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){
    593                 ISSMERROR("could not find element grids corresponding to quad icefront!");
    594         }
     580        grid1=element->GetNodeIndex(nodes[0]);
     581        grid2=element->GetNodeIndex(nodes[1]);
     582        grid3=element->GetNodeIndex(nodes[2]);
     583        grid4=element->GetNodeIndex(nodes[3]);
    595584
    596585        /*Build new xyz, bed and thickness lists for grids 1 to 4: */
    … …  
    653642
    654643        /*Free ressources:*/
    655         xfree((void**)&element_nodes);
    656644        xfree((void**)&doflist);
    657645
    … …  
    678666        /*Objects: */
    679667        double  pe_g[numdofs]={0.0};
    680         Node**  element_nodes=NULL;
    681668        Penta*   penta=NULL;
    682669
    … …  
    709696       
    710697        //Identify which grids are comprised in the quad:
    711         element_nodes=(Node**)xmalloc(6*sizeof(Node*));
    712         element->GetNodes((void**)element_nodes);
    713 
    714         grid1=-1; grid2=-1; grid3=-1; grid4=-1;
    715         for(i=0;i<6;i++){
    716                 if (nodes[0]==element_nodes[i])grid1=i;
    717                 if (nodes[1]==element_nodes[i])grid2=i;
    718                 if (nodes[2]==element_nodes[i])grid3=i;
    719                 if (nodes[3]==element_nodes[i])grid4=i;
    720         }
    721        
    722         if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){
    723                 ISSMERROR("could not find element grids corresponding to quad icefront!");
    724         }
     698        grid1=element->GetNodeIndex(nodes[0]);
     699        grid2=element->GetNodeIndex(nodes[1]);
     700        grid3=element->GetNodeIndex(nodes[2]);
     701        grid4=element->GetNodeIndex(nodes[3]);
    725702
    726703        /*Build new xyz, bed and thickness lists for grids 1 to 4: */
    … …  
    784761
    785762        /*Free ressources:*/
    786         xfree((void**)&element_nodes);
    787763        xfree((void**)&doflist);
    788764
    … …  
    1019995                                complete_list[j][2]=0;
    1020996                        }
    1021                         tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],l1l2l3_tria);
     997                        tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss);
    1022998                }
    1023999
    … …  
    12611237                                complete_list[j][2]=0;
    12621238                        }
    1263                         tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],l1l2l3_tria);
     1239                        tria->GetJacobianDeterminant2d(&J[i],&complete_list[0][0],gauss);
    12641240                }
    12651241
    … …  
    13101286
    13111287                } //for(i=0;i<4;i++)
    1312         } //for(ig=0;ig<num_gauss;ig++)
     1288        }
    13131289
    13141290        delete tria;
Note: See TracChangeset for help on using the changeset viewer.