Changeset 16838


Ignore:
Timestamp:
11/20/13 12:09:18 (11 years ago)
Author:
seroussi
Message:

CHG: starting to move Front vectors in analysis

Location:
issm/trunk-jpl/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r16831 r16838  
    11101110ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
    11111111
     1112        /*If no front, return NULL*/
     1113        if(!element->IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
     1114
     1115        /*Intermediaries*/
     1116        IssmDouble  rho_ice,rho_water,gravity;
     1117        IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure;
     1118        IssmDouble  surface_under_water,base_under_water,pressure;
     1119        IssmDouble *xyz_list = NULL;
     1120        IssmDouble *xyz_list_front = NULL;
     1121        IssmDouble  normal[2];
     1122
     1123        /*Fetch number of nodes for this finite element*/
     1124        int numnodes = element->GetNumberOfNodes();
     1125
     1126        /*Initialize Element vector and other vectors*/
     1127        ElementVector* pe    = element->NewElementVector(SSAApproximationEnum);
     1128        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
     1129
     1130        /*Retrieve all inputs and parameters*/
     1131        Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
     1132        Input* bed_input      =element->GetInput(BedEnum);       _assert_(bed_input);
     1133        rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum);
     1134        rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
     1135        gravity   = element->GetMaterialParameter(ConstantsGEnum);
     1136        element->GetVerticesCoordinates(&xyz_list);
     1137        element->ZeroLevelsetCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
     1138        element->NormalBase(&normal[0],xyz_list_front);
     1139
     1140        /*Start looping on Gaussian points*/
     1141        Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
     1142        for(int ig=gauss->begin();ig<gauss->end();ig++){
     1143
     1144                gauss->GaussPoint(ig);
     1145                thickness_input->GetInputValue(&thickness,gauss);
     1146                bed_input->GetInputValue(&bed,gauss);
     1147                element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
     1148                element->NodalFunctions(basis,gauss);
     1149
     1150                surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
     1151                base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
     1152                water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water);
     1153                ice_pressure   = 1.0/2.0*gravity*rho_ice*thickness*thickness;
     1154                pressure = ice_pressure + water_pressure;
     1155
     1156                for (int i=0;i<numnodes;i++){
     1157                        pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
     1158                        pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
     1159                }
     1160        }
     1161
     1162        /*Transform coordinate system*/
     1163        element->TransformLoadVectorCoord(pe,XYEnum);
     1164
     1165        /*Clean up and return*/
     1166        xDelete<IssmDouble>(xyz_list);
     1167        xDelete<IssmDouble>(xyz_list_front);
     1168        xDelete<IssmDouble>(basis);
     1169        delete gauss;
     1170        return pe;
    11121171        return NULL;
    11131172}/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r16833 r16838  
    8282                virtual void    GetDofListPressure(int** pdoflist,int setenum)=0;
    8383                virtual void   JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
     84                virtual void   JacobianDeterminantSurface(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss)=0;
    8485                virtual void   JacobianDeterminantBase(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0;
    8586                virtual void   JacobianDeterminantTop(IssmDouble* Jdet,IssmDouble* xyz_list_base,Gauss* gauss)=0;
    … …  
    130131                virtual Gauss* NewGauss(void)=0;
    131132                virtual Gauss* NewGauss(int order)=0;
     133      virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order)=0;
    132134                virtual Gauss* NewGaussBase(int order)=0;
    133135                virtual Gauss* NewGaussTop(int order)=0;
    … …  
    152154                virtual int    PressureInterpolation()=0;
    153155                virtual bool   IsZeroLevelset(int levelset_enum)=0;
     156                virtual void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum)=0;
    154157
    155158                #ifdef _HAVE_RESPONSES_

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r16836 r16838  
    884884/*}}}*/
    885885/*FUNCTION Penta::GetAreaCoordinates{{{*/
    886 void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[4][3],IssmDouble xyz_list[6][3],int numpoints){
     886void Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){
    887887        /*Computeportion of the element that is grounded*/
    888888
    … …  
    891891        IssmDouble  xyz_bis[3][3];
    892892
    893         area_init=fabs(xyz_list[1][0]*xyz_list[2][1] - xyz_list[1][1]*xyz_list[2][0] + xyz_list[0][0]*xyz_list[1][1] - xyz_list[0][1]*xyz_list[1][0] + xyz_list[2][0]*xyz_list[0][1] - xyz_list[2][1]*xyz_list[0][0])/2.;
     893        area_init=fabs(xyz_list[1*3+0]*xyz_list[2*3+1] - xyz_list[1*3+1]*xyz_list[2*3+0] + xyz_list[0*3+0]*xyz_list[1*3+1] - xyz_list[0*3+1]*xyz_list[1*3+0] + xyz_list[2*3+0]*xyz_list[0*3+1] - xyz_list[2*3+1]*xyz_list[0*3+0])/2.;
    894894
    895895        /*Initialize xyz_list with original xyz_list of triangle coordinates*/
    896896        for(j=0;j<3;j++){
    897897                for(k=0;k<3;k++){
    898                         xyz_bis[j][k]=xyz_list[j][k];
     898                        xyz_bis[j][k]=xyz_list[j*3+k];
    899899                }
    900900        }
    … …  
    903903                        for(k=0;k<3;k++){
    904904                                /*Change appropriate line*/
    905                                 xyz_bis[j][k]=xyz_zero[i][k];
     905                                xyz_bis[j][k]=xyz_zero[i*3+k];
    906906                        }
    907907
    … …  
    913913                        for(k=0;k<3;k++){
    914914                                /*Reinitialize xyz_list with original coordinates*/
    915                                 xyz_bis[j][k]=xyz_list[j][k];
     915                                xyz_bis[j][k]=xyz_list[j*3+k];
    916916                        }
    917917                }
    … …  
    15121512/*}}}*/
    15131513/*FUNCTION Penta::GetQuadNormal {{{*/
    1514 void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){
     1514void Penta:: GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list){
    15151515
    15161516        /*Build unit outward pointing vector*/
    … …  
    15191519        IssmDouble norm;
    15201520
    1521         AB[0]=xyz_list[1][0] - xyz_list[0][0];
    1522         AB[1]=xyz_list[1][1] - xyz_list[0][1];
    1523         AB[2]=xyz_list[1][2] - xyz_list[0][2];
    1524         AC[0]=xyz_list[2][0] - xyz_list[0][0];
    1525         AC[1]=xyz_list[2][1] - xyz_list[0][1];
    1526         AC[2]=xyz_list[2][2] - xyz_list[0][2];
     1521        AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
     1522        AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
     1523        AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
     1524        AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
     1525        AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
     1526        AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
    15271527
    15281528        cross(normal,AB,AC);
    … …  
    16481648}
    16491649/*}}}*/
    1650 /*FUNCTION Penta::GetZeroLevelsetCoordinates{{{*/
    1651 void Penta::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum){
     1650/*FUNCTION Penta::ZeroLevelsetCoordinates{{{*/
     1651void Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){
    16521652        /*Computeportion of the element that is grounded*/
    16531653
    … …  
    16551655        IssmDouble  s1,s2;
    16561656        IssmDouble  levelset[NUMVERTICES];
     1657        IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
    16571658
    16581659        /*Recover parameters and values*/
    … …  
    16661667                if(levelset[2]>0.) normal_orientation=1;
    16671668                /*New point 1*/
    1668                 xyz_zero[3*normal_orientation+0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
    1669                 xyz_zero[3*normal_orientation+1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
    1670                 xyz_zero[3*normal_orientation+2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
     1669                xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
     1670                xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
     1671                xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
    16711672
    16721673                /*New point 0*/
    1673                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
    1674                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
    1675                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
     1674                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
     1675                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
     1676                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
    16761677
    16771678                /*New point 3*/
    1678                 xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5][0]+s1*(xyz_list[4][0]-xyz_list[5][0]);
    1679                 xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5][1]+s1*(xyz_list[4][1]-xyz_list[5][1]);
    1680                 xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5][2]+s1*(xyz_list[4][2]-xyz_list[5][2]);
     1679                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5*3+0]+s1*(xyz_list[4*3+0]-xyz_list[5*3+0]);
     1680                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5*3+1]+s1*(xyz_list[4*3+1]-xyz_list[5*3+1]);
     1681                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5*3+2]+s1*(xyz_list[4*3+2]-xyz_list[5*3+2]);
    16811682
    16821683                /*New point 4*/
    1683                 xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5][0]+s2*(xyz_list[3][0]-xyz_list[5][0]);
    1684                 xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5][1]+s2*(xyz_list[3][1]-xyz_list[5][1]);
    1685                 xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5][2]+s2*(xyz_list[3][2]-xyz_list[5][2]);
     1684                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5*3+0]+s2*(xyz_list[3*3+0]-xyz_list[5*3+0]);
     1685                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5*3+1]+s2*(xyz_list[3*3+1]-xyz_list[5*3+1]);
     1686                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5*3+2]+s2*(xyz_list[3*3+2]-xyz_list[5*3+2]);
    16861687        }
    16871688        else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
    … …  
    16921693                if(levelset[0]>0.) normal_orientation=1;
    16931694                /*New point 1*/
    1694                 xyz_zero[3*normal_orientation+0]=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
    1695                 xyz_zero[3*normal_orientation+1]=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]);
    1696                 xyz_zero[3*normal_orientation+2]=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]);
     1695                xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
     1696                xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
     1697                xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
    16971698
    16981699                /*New point 2*/
    1699                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]);
    1700                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]);
    1701                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]);
     1700                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
     1701                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
     1702                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
    17021703
    17031704                /*New point 3*/
    1704                 xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3][0]+s1*(xyz_list[5][0]-xyz_list[3][0]);
    1705                 xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3][1]+s1*(xyz_list[5][1]-xyz_list[3][1]);
    1706                 xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3][2]+s1*(xyz_list[5][2]-xyz_list[3][2]);
     1705                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3*3+0]+s1*(xyz_list[5*3+0]-xyz_list[3*3+0]);
     1706                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3*3+1]+s1*(xyz_list[5*3+1]-xyz_list[3*3+1]);
     1707                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3*3+2]+s1*(xyz_list[5*3+2]-xyz_list[3*3+2]);
    17071708
    17081709                /*New point 4*/
    1709                 xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3][0]+s2*(xyz_list[4][0]-xyz_list[3][0]);
    1710                 xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3][1]+s2*(xyz_list[4][1]-xyz_list[3][1]);
    1711                 xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3][2]+s2*(xyz_list[4][2]-xyz_list[3][2]);
     1710                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3*3+0]+s2*(xyz_list[4*3+0]-xyz_list[3*3+0]);
     1711                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3*3+1]+s2*(xyz_list[4*3+1]-xyz_list[3*3+1]);
     1712                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3*3+2]+s2*(xyz_list[4*3+2]-xyz_list[3*3+2]);
    17121713        }
    17131714        else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
    … …  
    17181719                if(levelset[1]>0.) normal_orientation=1;
    17191720                /*New point 0*/
    1720                 xyz_zero[3*normal_orientation+0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
    1721                 xyz_zero[3*normal_orientation+1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
    1722                 xyz_zero[3*normal_orientation+2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
     1721                xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
     1722                xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
     1723                xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
    17231724
    17241725                /*New point 2*/
    1725                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
    1726                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
    1727                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
     1726                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
     1727                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
     1728                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
    17281729
    17291730                /*New point 3*/
    1730                 xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4][0]+s1*(xyz_list[3][0]-xyz_list[4][0]);
    1731                 xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4][1]+s1*(xyz_list[3][1]-xyz_list[4][1]);
    1732                 xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4][2]+s1*(xyz_list[3][2]-xyz_list[4][2]);
     1731                xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4*3+0]+s1*(xyz_list[3*3+0]-xyz_list[4*3+0]);
     1732                xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4*3+1]+s1*(xyz_list[3*3+1]-xyz_list[4*3+1]);
     1733                xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4*3+2]+s1*(xyz_list[3*3+2]-xyz_list[4*3+2]);
    17331734
    17341735                /*New point 4*/
    1735                 xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4][0]+s2*(xyz_list[5][0]-xyz_list[4][0]);
    1736                 xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4][1]+s2*(xyz_list[5][1]-xyz_list[4][1]);
    1737                 xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4][2]+s2*(xyz_list[5][2]-xyz_list[4][2]);
     1736                xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4*3+0]+s2*(xyz_list[5*3+0]-xyz_list[4*3+0]);
     1737                xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4*3+1]+s2*(xyz_list[5*3+1]-xyz_list[4*3+1]);
     1738                xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4*3+2]+s2*(xyz_list[5*3+2]-xyz_list[4*3+2]);
    17381739        }
    17391740        else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
    1740                 xyz_zero[3*0+0]=xyz_list[0][0];
    1741                 xyz_zero[3*0+1]=xyz_list[0][1];
    1742                 xyz_zero[3*0+2]=xyz_list[0][2];
     1741                xyz_zero[3*0+0]=xyz_list[0*3+0];
     1742                xyz_zero[3*0+1]=xyz_list[0*3+1];
     1743                xyz_zero[3*0+2]=xyz_list[0*3+2];
    17431744
    17441745                /*New point 2*/
    1745                 xyz_zero[3*1+0]=xyz_list[1][0];
    1746                 xyz_zero[3*1+1]=xyz_list[1][1];
    1747                 xyz_zero[3*1+2]=xyz_list[1][2];
     1746                xyz_zero[3*1+0]=xyz_list[1*3+0];
     1747                xyz_zero[3*1+1]=xyz_list[1*3+1];
     1748                xyz_zero[3*1+2]=xyz_list[1*3+2];
    17481749
    17491750                /*New point 3*/
    1750                 xyz_zero[3*2+0]=xyz_list[4][0];
    1751                 xyz_zero[3*2+1]=xyz_list[4][1];
    1752                 xyz_zero[3*2+2]=xyz_list[4][2];
     1751                xyz_zero[3*2+0]=xyz_list[4*3+0];
     1752                xyz_zero[3*2+1]=xyz_list[4*3+1];
     1753                xyz_zero[3*2+2]=xyz_list[4*3+2];
    17531754
    17541755                /*New point 4*/
    1755                 xyz_zero[3*3+0]=xyz_list[3][0];
    1756                 xyz_zero[3*3+1]=xyz_list[3][1];
    1757                 xyz_zero[3*3+2]=xyz_list[3][2];
     1756                xyz_zero[3*3+0]=xyz_list[3*3+0];
     1757                xyz_zero[3*3+1]=xyz_list[3*3+1];
     1758                xyz_zero[3*3+2]=xyz_list[3*3+2];
    17581759        }
    17591760        else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1760                 xyz_zero[3*0+0]=xyz_list[2][0];
    1761                 xyz_zero[3*0+1]=xyz_list[2][1];
    1762                 xyz_zero[3*0+2]=xyz_list[2][2];
     1761                xyz_zero[3*0+0]=xyz_list[2*3+0];
     1762                xyz_zero[3*0+1]=xyz_list[2*3+1];
     1763                xyz_zero[3*0+2]=xyz_list[2*3+2];
    17631764
    17641765                /*New point 2*/
    1765                 xyz_zero[3*1+0]=xyz_list[0][0];
    1766                 xyz_zero[3*1+1]=xyz_list[0][1];
    1767                 xyz_zero[3*1+2]=xyz_list[0][2];
     1766                xyz_zero[3*1+0]=xyz_list[0*3+0];
     1767                xyz_zero[3*1+1]=xyz_list[0*3+1];
     1768                xyz_zero[3*1+2]=xyz_list[0*3+2];
    17681769
    17691770                /*New point 3*/
    1770                 xyz_zero[3*2+0]=xyz_list[3][0];
    1771                 xyz_zero[3*2+1]=xyz_list[3][1];
    1772                 xyz_zero[3*2+2]=xyz_list[3][2];
     1771                xyz_zero[3*2+0]=xyz_list[3*3+0];
     1772                xyz_zero[3*2+1]=xyz_list[3*3+1];
     1773                xyz_zero[3*2+2]=xyz_list[3*3+2];
    17731774
    17741775                /*New point 4*/
    1775                 xyz_zero[3*3+0]=xyz_list[5][0];
    1776                 xyz_zero[3*3+1]=xyz_list[5][1];
    1777                 xyz_zero[3*3+2]=xyz_list[5][2];
     1776                xyz_zero[3*3+0]=xyz_list[5*3+0];
     1777                xyz_zero[3*3+1]=xyz_list[5*3+1];
     1778                xyz_zero[3*3+2]=xyz_list[5*3+2];
    17781779        }
    17791780        else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1780                 xyz_zero[3*0+0]=xyz_list[1][0];
    1781                 xyz_zero[3*0+1]=xyz_list[1][1];
    1782                 xyz_zero[3*0+2]=xyz_list[1][2];
     1781                xyz_zero[3*0+0]=xyz_list[1*3+0];
     1782                xyz_zero[3*0+1]=xyz_list[1*3+1];
     1783                xyz_zero[3*0+2]=xyz_list[1*3+2];
    17831784
    17841785                /*New point 2*/
    1785                 xyz_zero[3*1+0]=xyz_list[2][0];
    1786                 xyz_zero[3*1+1]=xyz_list[2][1];
    1787                 xyz_zero[3*1+2]=xyz_list[2][2];
     1786                xyz_zero[3*1+0]=xyz_list[2*3+0];
     1787                xyz_zero[3*1+1]=xyz_list[2*3+1];
     1788                xyz_zero[3*1+2]=xyz_list[2*3+2];
    17881789
    17891790                /*New point 3*/
    1790                 xyz_zero[3*2+0]=xyz_list[5][0];
    1791                 xyz_zero[3*2+1]=xyz_list[5][1];
    1792                 xyz_zero[3*2+2]=xyz_list[5][2];
     1791                xyz_zero[3*2+0]=xyz_list[5*3+0];
     1792                xyz_zero[3*2+1]=xyz_list[5*3+1];
     1793                xyz_zero[3*2+2]=xyz_list[5*3+2];
    17931794
    17941795                /*New point 4*/
    1795                 xyz_zero[3*3+0]=xyz_list[4][0];
    1796                 xyz_zero[3*3+1]=xyz_list[4][1];
    1797                 xyz_zero[3*3+2]=xyz_list[4][2];
     1796                xyz_zero[3*3+0]=xyz_list[4*3+0];
     1797                xyz_zero[3*3+1]=xyz_list[4*3+1];
     1798                xyz_zero[3*3+2]=xyz_list[4*3+2];
    17981799        }
    17991800        else _error_("Case not covered");
     1801
     1802        /*Assign output pointer*/
     1803        *pxyz_zero= xyz_zero;
    18001804}
    18011805/*}}}*/
    … …  
    87648768        IssmDouble  surface_under_water,base_under_water,pressure;
    87658769        IssmDouble  xyz_list[NUMVERTICES][3];
    8766         IssmDouble  xyz_list_front[4][3];
     8770        IssmDouble* xyz_list_front = NULL;
    87678771        IssmDouble  area_coordinates[4][3];
    87688772        IssmDouble  normal[3];
    … …  
    87828786        rho_ice  =matpar->GetRhoIce();
    87838787        gravity  =matpar->GetG();
    8784         GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
    8785         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
     8788        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
     8789        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
    87868790        GetQuadNormal(&normal[0],xyz_list_front);
    87878791
    … …  
    88548858        IssmDouble  surface_under_water,base_under_water,pressure;
    88558859        IssmDouble  xyz_list[NUMVERTICES][3];
    8856         IssmDouble  xyz_list_front[4][3];
     8860        IssmDouble* xyz_list_front = NULL;
    88578861        IssmDouble  area_coordinates[4][3];
    88588862        IssmDouble  normal[3];
    … …  
    88698873
    88708874        /*Initialize Element matrix and vectors*/
     8875        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    88718876        ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    88728877        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    88738878
    88748879        /*Retrieve all inputs and parameters*/
    8875         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    88768880        rho_water=matpar->GetRhoWater();
    88778881        rho_ice  =matpar->GetRhoIce();
    88788882        gravity  =matpar->GetG();
    8879         GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
    8880         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
     8883        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
     8884        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
    88818885        GetQuadNormal(&normal[0],xyz_list_front);
    88828886
    … …  
    89108914        /*Clean up and return*/
    89118915        xDelete<int>(cs_list);
     8916        xDelete<IssmDouble>(xyz_list_front);
    89128917        xDelete<IssmDouble>(vbasis);
    89138918        delete gauss;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r16833 r16838  
    114114                int    PressureInterpolation();
    115115                bool   IsZeroLevelset(int levelset_enum);
     116                void   ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    116117
    117118                void   ResultInterpolation(int* pinterpolation,int output_enum);
    … …  
    211212                IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
    212213                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    213                 void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints);
     214                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
    214215
    215216                void             GetVertexPidList(int* doflist);
    … …  
    231232                void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
    232233                void             GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
    233                 void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
     234                void           GetQuadNormal(IssmDouble* normal,IssmDouble* xyz_list);
    234235                void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
    235236                void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    236237                Penta*         GetUpperElement(void);
    237                 void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
    238238                Penta*         GetLowerElement(void);
    239239                void           InputChangeName(int input_enum, int enum_type_old);
    … …  
    247247                bool           IsNodeOnShelfFromFlags(IssmDouble* flags);
    248248                void           JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
     249                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    249250                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    250251                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
    … …  
    252253                Gauss*         NewGauss(void);
    253254                Gauss*         NewGauss(int order);
     255      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
    254256                Gauss*         NewGaussBase(int order);
    255257                Gauss*         NewGaussTop(int order);

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

    r16818 r16838  
    19451945/*}}}*/
    19461946/*FUNCTION PentaRef::GetQuadJacobianDeterminant{{{*/
    1947 void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble xyz_list[4][3],Gauss* gauss){
     1947void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble* xyz_list,Gauss* gauss){
    19481948        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    19491949
    19501950        IssmDouble x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
    19511951
    1952         x1=xyz_list[0][0];
    1953         y1=xyz_list[0][1];
    1954         z1=xyz_list[0][2];
    1955         x2=xyz_list[1][0];
    1956         y2=xyz_list[1][1];
    1957         z2=xyz_list[1][2];
    1958         x3=xyz_list[2][0];
    1959         y3=xyz_list[2][1];
    1960         z3=xyz_list[2][2];
    1961         x4=xyz_list[3][0];
    1962         y4=xyz_list[3][1];
    1963         z4=xyz_list[3][2];
     1952        x1=xyz_list[0*3+0];
     1953        y1=xyz_list[0*3+1];
     1954        z1=xyz_list[0*3+2];
     1955        x2=xyz_list[1*3+0];
     1956        y2=xyz_list[1*3+1];
     1957        z2=xyz_list[1*3+2];
     1958        x3=xyz_list[2*3+0];
     1959        y3=xyz_list[2*3+1];
     1960        z3=xyz_list[2*3+2];
     1961        x4=xyz_list[3*3+0];
     1962        y4=xyz_list[3*3+1];
     1963        z4=xyz_list[3*3+2];
    19641964
    19651965        /*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/

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

    r16818 r16838  
    3535                void GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,Gauss* gauss);
    3636                void GetNodalFunctionsMINIDerivativesReference(IssmDouble* dl1dl7,Gauss* gauss);
    37                 void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble xyz_list[4][3],Gauss* gauss);
     37                void GetQuadJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
    3838                void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
    3939                void GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r16833 r16838  
    114114                bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
    115115                void        JacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
     116                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
    116117                void        JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    117118                void        JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    … …  
    142143                Gauss*      NewGauss(void){_error_("not implemented yet");};
    143144                Gauss*      NewGauss(int order){_error_("not implemented yet");};
     145      Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
    144146                Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
    145147                Gauss*      NewGaussTop(int order){_error_("not implemented yet");};
    … …  
    151153                void        ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
    152154                bool        IsZeroLevelset(int levelset_enum){_error_("not implemented");};
     155                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){_error_("not implemented");};
    153156
    154157                #ifdef _HAVE_THERMAL_

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r16834 r16838  
    879879/*}}}*/
    880880/*FUNCTION Tria::GetAreaCoordinates{{{*/
    881 void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints){
     881void Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){
    882882        /*Computeportion of the element that is grounded*/
    883883
    … …  
    891891        for(j=0;j<3;j++){
    892892                for(k=0;k<3;k++){
    893                         xyz_bis[j][k]=xyz_list[j][k];
     893                        xyz_bis[j][k]=xyz_list[j*3+k];
    894894                }
    895895        }
    … …  
    898898                        for(k=0;k<3;k++){
    899899                                /*Change appropriate line*/
    900                                 xyz_bis[j][k]=xyz_zero[i][k];
     900                                xyz_bis[j][k]=xyz_zero[i*3+k];
    901901                        }
    902902
    … …  
    908908                        for(k=0;k<3;k++){
    909909                                /*Reinitialize xyz_list with original coordinates*/
    910                                 xyz_bis[j][k]=xyz_list[j][k];
     910                                xyz_bis[j][k]=xyz_list[j*3+k];
    911911                        }
    912912                }
    … …  
    11841184}
    11851185/*}}}*/
    1186 /*FUNCTION Tria::GetZeroLevelsetCoordinates{{{*/
    1187 void Tria::GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum){
    1188         /*Computeportion of the element that is grounded*/
     1186/*FUNCTION Tria::ZeroLevelsetCoordinates{{{*/
     1187void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){
    11891188
    11901189        int         normal_orientation=0;
    … …  
    11931192
    11941193        /*Recover parameters and values*/
     1194        IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
    11951195        GetInputListOnVertices(&levelset[0],levelsetenum);
    11961196
    … …  
    12021202                if(levelset[2]>0.) normal_orientation=1;
    12031203                /*New point 1*/
    1204                 xyz_zero[3*normal_orientation+0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]);
    1205                 xyz_zero[3*normal_orientation+1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]);
    1206                 xyz_zero[3*normal_orientation+2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]);
     1204                xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
     1205                xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
     1206                xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
    12071207
    12081208                /*New point 0*/
    1209                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]);
    1210                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]);
    1211                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]);
     1209                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
     1210                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
     1211                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
    12121212        }
    12131213        else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
    … …  
    12181218                if(levelset[0]>0.) normal_orientation=1;
    12191219                /*New point 1*/
    1220                 xyz_zero[3*normal_orientation+0]=xyz_list[0][0]+s1*(xyz_list[2][0]-xyz_list[0][0]);
    1221                 xyz_zero[3*normal_orientation+1]=xyz_list[0][1]+s1*(xyz_list[2][1]-xyz_list[0][1]);
    1222                 xyz_zero[3*normal_orientation+2]=xyz_list[0][2]+s1*(xyz_list[2][2]-xyz_list[0][2]);
     1220                xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
     1221                xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
     1222                xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
    12231223
    12241224                /*New point 2*/
    1225                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0][0]+s2*(xyz_list[1][0]-xyz_list[0][0]);
    1226                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0][1]+s2*(xyz_list[1][1]-xyz_list[0][1]);
    1227                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0][2]+s2*(xyz_list[1][2]-xyz_list[0][2]);
     1225                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
     1226                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
     1227                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
    12281228        }
    12291229        else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
    … …  
    12341234                if(levelset[1]>0.) normal_orientation=1;
    12351235                /*New point 0*/
    1236                 xyz_zero[3*normal_orientation+0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]);
    1237                 xyz_zero[3*normal_orientation+1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]);
    1238                 xyz_zero[3*normal_orientation+2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]);
     1236                xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
     1237                xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
     1238                xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
    12391239
    12401240                /*New point 2*/
    1241                 xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]);
    1242                 xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]);
    1243                 xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]);
     1241                xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
     1242                xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
     1243                xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
    12441244        }
    12451245        else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
    1246                 xyz_zero[3*0+0]=xyz_list[0][0];
    1247                 xyz_zero[3*0+1]=xyz_list[0][1];
    1248                 xyz_zero[3*0+2]=xyz_list[0][2];
     1246                xyz_zero[3*0+0]=xyz_list[0*3+0];
     1247                xyz_zero[3*0+1]=xyz_list[0*3+1];
     1248                xyz_zero[3*0+2]=xyz_list[0*3+2];
    12491249
    12501250                /*New point 2*/
    1251                 xyz_zero[3*1+0]=xyz_list[1][0];
    1252                 xyz_zero[3*1+1]=xyz_list[1][1];
    1253                 xyz_zero[3*1+2]=xyz_list[1][2];
     1251                xyz_zero[3*1+0]=xyz_list[1*3+0];
     1252                xyz_zero[3*1+1]=xyz_list[1*3+1];
     1253                xyz_zero[3*1+2]=xyz_list[1*3+2];
    12541254        }
    12551255        else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1256                 xyz_zero[3*0+0]=xyz_list[2][0];
    1257                 xyz_zero[3*0+1]=xyz_list[2][1];
    1258                 xyz_zero[3*0+2]=xyz_list[2][2];
     1256                xyz_zero[3*0+0]=xyz_list[2*3+0];
     1257                xyz_zero[3*0+1]=xyz_list[2*3+1];
     1258                xyz_zero[3*0+2]=xyz_list[2*3+2];
    12591259
    12601260                /*New point 2*/
    1261                 xyz_zero[3*1+0]=xyz_list[0][0];
    1262                 xyz_zero[3*1+1]=xyz_list[0][1];
    1263                 xyz_zero[3*1+2]=xyz_list[0][2];
     1261                xyz_zero[3*1+0]=xyz_list[0*3+0];
     1262                xyz_zero[3*1+1]=xyz_list[0*3+1];
     1263                xyz_zero[3*1+2]=xyz_list[0*3+2];
    12641264        }
    12651265        else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
    1266                 xyz_zero[3*0+0]=xyz_list[1][0];
    1267                 xyz_zero[3*0+1]=xyz_list[1][1];
    1268                 xyz_zero[3*0+2]=xyz_list[1][2];
     1266                xyz_zero[3*0+0]=xyz_list[1*3+0];
     1267                xyz_zero[3*0+1]=xyz_list[1*3+1];
     1268                xyz_zero[3*0+2]=xyz_list[1*3+2];
    12691269
    12701270                /*New point 2*/
    1271                 xyz_zero[3*1+0]=xyz_list[2][0];
    1272                 xyz_zero[3*1+1]=xyz_list[2][1];
    1273                 xyz_zero[3*1+2]=xyz_list[2][2];
     1271                xyz_zero[3*1+0]=xyz_list[2*3+0];
     1272                xyz_zero[3*1+1]=xyz_list[2*3+1];
     1273                xyz_zero[3*1+2]=xyz_list[2*3+2];
    12741274        }
    12751275        else _error_("Case not covered");
     1276
     1277        /*Assign output pointer*/
     1278        *pxyz_zero= xyz_zero;
    12761279}
    12771280/*}}}*/
    … …  
    20142017}
    20152018/*}}}*/
     2019/*FUNCTION Tria::JacobianDeterminantSurface{{{*/
     2020void Tria::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){
     2021
     2022        _assert_(gauss->Enum()==GaussTriaEnum);
     2023        this->GetSegmentJacobianDeterminant(pJdet,xyz_list,(GaussTria*)gauss);
     2024
     2025}
     2026/*}}}*/
    20162027/*FUNCTION Tria::HasEdgeOnBed {{{*/
    20172028bool Tria::HasEdgeOnBed(){
    … …  
    22072218Gauss* Tria::NewGauss(int order){
    22082219        return new GaussTria(order);
     2220}
     2221/*}}}*/
     2222/*FUNCTION Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){{{*/
     2223Gauss* Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){
     2224
     2225        IssmDouble  area_coordinates[2][3];
     2226
     2227        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
     2228
     2229        return new GaussTria(area_coordinates,order);
    22092230}
    22102231/*}}}*/
    … …  
    39303951        IssmDouble  Jdet,water_pressure,air_pressure,pressure;
    39313952        IssmDouble  xyz_list[NUMVERTICES][3];
    3932         IssmDouble  xyz_list_front[2][3];
     3953        IssmDouble* xyz_list_front = NULL;
    39333954        IssmDouble  area_coordinates[2][3];
    39343955        IssmDouble  normal[2];
    … …  
    39533974        rho_ice = matpar->GetRhoIce();
    39543975        gravity   = matpar->GetG();
    3955         GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
    3956         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
    3957         NormalBase(&normal[0],&xyz_list_front[0][0]);
     3976        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
     3977        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
     3978        NormalBase(&normal[0],xyz_list_front);
    39583979
    39593980        /*Start looping on Gaussian points*/
    3960         IssmDouble ymax=max(xyz_list_front[0][1],xyz_list_front[1][1]);
    3961         IssmDouble ymin=min(xyz_list_front[0][1],xyz_list_front[1][1]);
     3981        IssmDouble ymax=max(xyz_list_front[0*3+1],xyz_list_front[1*3+1]);
     3982        IssmDouble ymin=min(xyz_list_front[0*3+1],xyz_list_front[1*3+1]);
    39623983        if(ymax>0. && ymin<0.) gauss=new GaussTria(area_coordinates,30); //refined in vertical because of the sea level discontinuity
    39633984        else                   gauss=new GaussTria(area_coordinates,3);
    … …  
    39683989                y_g=GetYcoord(gauss);
    39693990                GetNodalFunctionsVelocity(vbasis,gauss);
    3970                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
     3991                GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss);
    39713992
    39723993                water_pressure = rho_water*gravity*min(0.,y_g); //0 if the gaussian point is above water level
    … …  
    39804001        }
    39814002
    3982 
    39834003        /*Transform coordinate system*/
    39844004        ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    … …  
    39864006        /*Clean up and return*/
    39874007        xDelete<int>(cs_list);
     4008        xDelete<IssmDouble>(xyz_list_front);
    39884009        xDelete<IssmDouble>(vbasis);
    39894010        delete gauss;
    … …  
    42014222        IssmDouble  surface_under_water,base_under_water,pressure;
    42024223        IssmDouble  xyz_list[NUMVERTICES][3];
    4203         IssmDouble  xyz_list_front[2][3];
     4224        IssmDouble  *xyz_list_front = NULL;
    42044225        IssmDouble  area_coordinates[2][3];
    42054226        IssmDouble  normal[2];
    … …  
    42194240        rho_ice   = matpar->GetRhoIce();
    42204241        gravity   = matpar->GetG();
    4221         GetZeroLevelsetCoordinates(&xyz_list_front[0][0],xyz_list,MaskIceLevelsetEnum);
    4222         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
    4223         NormalBase(&normal[0],&xyz_list_front[0][0]);
     4242        ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
     4243        GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
     4244        NormalBase(&normal[0],xyz_list_front);
    42244245
    42254246        /*Start looping on Gaussian points*/
    … …  
    42304251                thickness_input->GetInputValue(&thickness,gauss);
    42314252                bed_input->GetInputValue(&bed,gauss);
    4232                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_front[0][0],gauss);
     4253                GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss);
    42334254                GetNodalFunctions(basis,gauss);
    42344255
    … …  
    42514272        /*Clean up and return*/
    42524273        xDelete<IssmDouble>(basis);
     4274        xDelete<IssmDouble>(xyz_list_front);
    42534275        delete gauss;
    42544276        return pe;

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16833 r16838  
    135135                void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    136136                IssmDouble  TimeAdapt();
     137                void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
    137138                bool        IsZeroLevelset(int levelset_enum);
    138139
    … …  
    248249                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    249250                IssmDouble     GetArea(void);
    250                 void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints);
     251                void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
    251252                int            GetElementType(void);
    252253
    … …  
    259260                void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
    260261                IssmDouble     GetMaterialParameter(int enum_in);
    261                 void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[3][3],int levelsetenum);
    262262                Input*         GetInput(int inputenum);
    263263                void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
    … …  
    277277                bool             IsInput(int name);
    278278                void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
     279                void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
    279280                void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    280281                void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");};
    … …  
    282283                Gauss*         NewGauss(void);
    283284                Gauss*         NewGauss(int order);
     285      Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
    284286                Gauss*         NewGaussBase(int order){_error_("not implemented yet");};
    285287                Gauss*         NewGaussTop(int order){_error_("not implemented yet");};
Note: See TracChangeset for help on using the changeset viewer.