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

CHG: starting to move Front vectors in analysis

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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;
Note: See TracChangeset for help on using the changeset viewer.