Changeset 25723
- Timestamp:
- 10/29/20 14:30:20 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r25722 r25723 425 425 IssmDouble xyz_front[2][3]; 426 426 427 IssmDouble *xyz_list = NULL;428 this->GetVerticesCoordinates(&xyz_list);427 IssmDouble xyz_list[NUMVERTICES][3]; 428 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 429 429 430 430 /*Recover parameters and values*/ … … 446 446 pt1 = 1; pt2 = 0; 447 447 } 448 xyz_front[pt2][0]=xyz_list[ 3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);449 xyz_front[pt2][1]=xyz_list[ 3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);450 xyz_front[pt2][2]=xyz_list[ 3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);451 xyz_front[pt1][0]=xyz_list[ 3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);452 xyz_front[pt1][1]=xyz_list[ 3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);453 xyz_front[pt1][2]=xyz_list[ 3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);448 xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]); 449 xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]); 450 xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]); 451 xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]); 452 xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]); 453 xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]); 454 454 } 455 455 else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 … … 462 462 } 463 463 464 xyz_front[pt1][0]=xyz_list[ 3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);465 xyz_front[pt1][1]=xyz_list[ 3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);466 xyz_front[pt1][2]=xyz_list[ 3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);467 xyz_front[pt2][0]=xyz_list[ 3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);468 xyz_front[pt2][1]=xyz_list[ 3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);469 xyz_front[pt2][2]=xyz_list[ 3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);464 xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]); 465 xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]); 466 xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]); 467 xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]); 468 xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]); 469 xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]); 470 470 } 471 471 else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 … … 478 478 } 479 479 480 xyz_front[pt2][0]=xyz_list[ 3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);481 xyz_front[pt2][1]=xyz_list[ 3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);482 xyz_front[pt2][2]=xyz_list[ 3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);483 xyz_front[pt1][0]=xyz_list[ 3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);484 xyz_front[pt1][1]=xyz_list[ 3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);485 xyz_front[pt1][2]=xyz_list[ 3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);480 xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]); 481 xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]); 482 xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]); 483 xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]); 484 xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]); 485 xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]); 486 486 } 487 487 else{ … … 511 511 512 512 /*Start looping on Gaussian points*/ 513 Gauss* gauss=this->NewGaussBase( xyz_list,&xyz_front[0][0],3);513 Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3); 514 514 while(gauss->next()){ 515 515 thickness_input->GetInputValue(&thickness,gauss); … … 545 545 IssmDouble xyz_front[2][3]; 546 546 547 IssmDouble *xyz_list = NULL;548 this->GetVerticesCoordinates(&xyz_list);547 IssmDouble xyz_list[NUMVERTICES][3]; 548 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 549 549 550 550 /*Recover parameters and values*/ … … 566 566 pt1 = 1; pt2 = 0; 567 567 } 568 xyz_front[pt2][0]=xyz_list[ 3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);569 xyz_front[pt2][1]=xyz_list[ 3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);570 xyz_front[pt2][2]=xyz_list[ 3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);571 xyz_front[pt1][0]=xyz_list[ 3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);572 xyz_front[pt1][1]=xyz_list[ 3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);573 xyz_front[pt1][2]=xyz_list[ 3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);568 xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]); 569 xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]); 570 xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]); 571 xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]); 572 xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]); 573 xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]); 574 574 } 575 575 else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 … … 582 582 } 583 583 584 xyz_front[pt1][0]=xyz_list[ 3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);585 xyz_front[pt1][1]=xyz_list[ 3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);586 xyz_front[pt1][2]=xyz_list[ 3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);587 xyz_front[pt2][0]=xyz_list[ 3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);588 xyz_front[pt2][1]=xyz_list[ 3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);589 xyz_front[pt2][2]=xyz_list[ 3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);584 xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]); 585 xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]); 586 xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]); 587 xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]); 588 xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]); 589 xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]); 590 590 } 591 591 else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 … … 598 598 } 599 599 600 xyz_front[pt2][0]=xyz_list[ 3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);601 xyz_front[pt2][1]=xyz_list[ 3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);602 xyz_front[pt2][2]=xyz_list[ 3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);603 xyz_front[pt1][0]=xyz_list[ 3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);604 xyz_front[pt1][1]=xyz_list[ 3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);605 xyz_front[pt1][2]=xyz_list[ 3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);600 xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]); 601 xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]); 602 xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]); 603 xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]); 604 xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]); 605 xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]); 606 606 } 607 607 else{ … … 637 637 638 638 /*Start looping on Gaussian points*/ 639 Gauss* gauss=this->NewGaussBase( xyz_list,&xyz_front[0][0],3);639 Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3); 640 640 while(gauss->next()){ 641 641 thickness_input->GetInputValue(&thickness,gauss); … … 1115 1115 1116 1116 /*Intermediaries*/ 1117 IssmDouble* xyz_list = NULL;1118 1117 IssmDouble bed_normal[3],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES]; 1119 1118 IssmDouble water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES]; … … 1131 1130 1132 1131 /* Get node coordinates and dof list: */ 1133 GetVerticesCoordinates(&xyz_list); 1132 IssmDouble xyz_list[NUMVERTICES][3]; 1133 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 1134 1134 1135 1135 /*Retrieve all inputs we will be needing: */ … … 1144 1144 1145 1145 /*Compute strain rate viscosity and pressure: */ 1146 this->StrainRateFS(&epsilon[0], xyz_list,gauss,vx_input,vy_input,vz_input);1147 this->material->ViscosityFS(&viscosity,3, xyz_list,gauss,vx_input,vy_input,vz_input);1146 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 1147 this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 1148 1148 /*FIXME: this is for Hongju only*/ 1149 1149 //pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]); … … 1164 1164 /*If was grounded*/ 1165 1165 if (phi[i]>=0.){ 1166 NormalBase(&bed_normal[0], xyz_list);1166 NormalBase(&bed_normal[0],&xyz_list[0][0]); 1167 1167 sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1] + sigmazz[i]*bed_normal[2]*bed_normal[2]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]+2*sigmaxz[i]*bed_normal[0]*bed_normal[2]+2*sigmayz[i]*bed_normal[1]*bed_normal[2]); 1168 1168 water_pressure[i]=-gravity*rho_water*base[i]; … … 1183 1183 /*clean up*/ 1184 1184 delete gauss; 1185 xDelete<IssmDouble>(xyz_list);1186 1185 } 1187 1186 /*}}}*/ … … 1896 1895 IssmDouble xyz_front[2][3]; 1897 1896 1898 IssmDouble *xyz_list = NULL;1899 this->GetVerticesCoordinates(&xyz_list);1897 IssmDouble xyz_list[NUMVERTICES][3]; 1898 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 1900 1899 1901 1900 /*Recover parameters and values*/ … … 1917 1916 pt1 = 1; pt2 = 0; 1918 1917 } 1919 xyz_front[pt2][0]=xyz_list[ 3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);1920 xyz_front[pt2][1]=xyz_list[ 3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);1921 xyz_front[pt2][2]=xyz_list[ 3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);1922 xyz_front[pt1][0]=xyz_list[ 3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);1923 xyz_front[pt1][1]=xyz_list[ 3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);1924 xyz_front[pt1][2]=xyz_list[ 3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);1918 xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]); 1919 xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]); 1920 xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]); 1921 xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]); 1922 xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]); 1923 xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]); 1925 1924 } 1926 1925 else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 … … 1933 1932 } 1934 1933 1935 xyz_front[pt1][0]=xyz_list[ 3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);1936 xyz_front[pt1][1]=xyz_list[ 3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);1937 xyz_front[pt1][2]=xyz_list[ 3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);1938 xyz_front[pt2][0]=xyz_list[ 3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);1939 xyz_front[pt2][1]=xyz_list[ 3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);1940 xyz_front[pt2][2]=xyz_list[ 3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);1934 xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]); 1935 xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]); 1936 xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]); 1937 xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]); 1938 xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]); 1939 xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]); 1941 1940 } 1942 1941 else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 … … 1949 1948 } 1950 1949 1951 xyz_front[pt2][0]=xyz_list[ 3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);1952 xyz_front[pt2][1]=xyz_list[ 3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);1953 xyz_front[pt2][2]=xyz_list[ 3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);1954 xyz_front[pt1][0]=xyz_list[ 3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);1955 xyz_front[pt1][1]=xyz_list[ 3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);1956 xyz_front[pt1][2]=xyz_list[ 3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);1950 xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]); 1951 xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]); 1952 xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]); 1953 xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]); 1954 xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]); 1955 xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]); 1957 1956 } 1958 1957 else{ … … 1983 1982 1984 1983 /*Start looping on Gaussian points*/ 1985 Gauss* gauss=this->NewGaussBase( xyz_list,&xyz_front[0][0],3);1984 Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3); 1986 1985 while(gauss->next()){ 1987 1986 thickness_input->GetInputValue(&thickness,gauss); … … 1995 1994 /*Clean up and return*/ 1996 1995 delete gauss; 1997 xDelete<IssmDouble>(xyz_list);1998 1996 return flux; 1999 2000 1997 } 2001 1998 /*}}}*/ … … 2014 2011 IssmDouble xyz_front[2][3]; 2015 2012 2016 IssmDouble *xyz_list = NULL;2017 this->GetVerticesCoordinates(&xyz_list);2013 IssmDouble xyz_list[NUMVERTICES][3]; 2014 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 2018 2015 2019 2016 /*Recover parameters and values*/ … … 2035 2032 pt1 = 1; pt2 = 0; 2036 2033 } 2037 xyz_front[pt2][0]=xyz_list[ 3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);2038 xyz_front[pt2][1]=xyz_list[ 3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);2039 xyz_front[pt2][2]=xyz_list[ 3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);2040 xyz_front[pt1][0]=xyz_list[ 3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);2041 xyz_front[pt1][1]=xyz_list[ 3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);2042 xyz_front[pt1][2]=xyz_list[ 3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);2034 xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]); 2035 xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]); 2036 xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]); 2037 xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]); 2038 xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]); 2039 xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]); 2043 2040 } 2044 2041 else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 … … 2051 2048 } 2052 2049 2053 xyz_front[pt1][0]=xyz_list[ 3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);2054 xyz_front[pt1][1]=xyz_list[ 3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);2055 xyz_front[pt1][2]=xyz_list[ 3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);2056 xyz_front[pt2][0]=xyz_list[ 3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);2057 xyz_front[pt2][1]=xyz_list[ 3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);2058 xyz_front[pt2][2]=xyz_list[ 3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);2050 xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]); 2051 xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]); 2052 xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]); 2053 xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]); 2054 xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]); 2055 xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]); 2059 2056 } 2060 2057 else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 … … 2067 2064 } 2068 2065 2069 xyz_front[pt2][0]=xyz_list[ 3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);2070 xyz_front[pt2][1]=xyz_list[ 3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);2071 xyz_front[pt2][2]=xyz_list[ 3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);2072 xyz_front[pt1][0]=xyz_list[ 3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);2073 xyz_front[pt1][1]=xyz_list[ 3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);2074 xyz_front[pt1][2]=xyz_list[ 3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);2066 xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]); 2067 xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]); 2068 xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]); 2069 xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]); 2070 xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]); 2071 xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]); 2075 2072 } 2076 2073 else{ … … 2103 2100 2104 2101 /*Start looping on Gaussian points*/ 2105 Gauss* gauss=this->NewGaussBase( xyz_list,&xyz_front[0][0],3);2102 Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3); 2106 2103 while(gauss->next()){ 2107 2104 thickness_input->GetInputValue(&thickness,gauss); … … 2114 2111 2115 2112 /*Clean up and return*/ 2116 xDelete<IssmDouble>(xyz_list);2117 2113 delete gauss; 2118 2114 return flux; 2119 2120 2115 } 2121 2116 /*}}}*/ … … 3418 3413 void Penta::StrainRateparallel(){/*{{{*/ 3419 3414 3420 IssmDouble *xyz_list = NULL;3421 3415 IssmDouble epsilon[6]; 3422 3416 GaussPenta* gauss=NULL; … … 3428 3422 3429 3423 /* Get node coordinates and dof list: */ 3430 this->GetVerticesCoordinates(&xyz_list); 3424 IssmDouble xyz_list[NUMVERTICES][3]; 3425 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 3431 3426 3432 3427 /*Retrieve all inputs we will need*/ … … 3446 3441 3447 3442 /*Compute strain rate and viscosity: */ 3448 this->StrainRateFS(&epsilon[0], xyz_list,gauss,vx_input,vy_input,vz_input);3443 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 3449 3444 strainxx=epsilon[0]; 3450 3445 strainyy=epsilon[1]; … … 3460 3455 /*Clean up and return*/ 3461 3456 delete gauss; 3462 xDelete<IssmDouble>(xyz_list);3463 3457 } 3464 3458 /*}}}*/ 3465 3459 void Penta::StrainRateperpendicular(){/*{{{*/ 3466 3460 3467 IssmDouble *xyz_list = NULL;3468 3461 IssmDouble epsilon[6]; 3469 3462 GaussPenta* gauss=NULL; … … 3475 3468 3476 3469 /* Get node coordinates and dof list: */ 3477 this->GetVerticesCoordinates(&xyz_list); 3470 IssmDouble xyz_list[NUMVERTICES][3]; 3471 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 3478 3472 3479 3473 /*Retrieve all inputs we will need*/ … … 3493 3487 3494 3488 /*Compute strain rate and viscosity: */ 3495 this->StrainRateFS(&epsilon[0], xyz_list,gauss,vx_input,vy_input,vz_input);3489 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); 3496 3490 strainxx=epsilon[0]; 3497 3491 strainyy=epsilon[1]; … … 3507 3501 /*Clean up and return*/ 3508 3502 delete gauss; 3509 xDelete<IssmDouble>(xyz_list);3510 3503 } 3511 3504 /*}}}*/ … … 3720 3713 IssmDouble xyz_front[2][3]; 3721 3714 3722 IssmDouble *xyz_list = NULL;3723 this->GetVerticesCoordinates(&xyz_list);3715 IssmDouble xyz_list[NUMVERTICES][3]; 3716 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 3724 3717 3725 3718 /*Recover parameters and values*/ … … 3741 3734 pt1 = 1; pt2 = 0; 3742 3735 } 3743 xyz_front[pt2][0]=xyz_list[ 3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);3744 xyz_front[pt2][1]=xyz_list[ 3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);3745 xyz_front[pt2][2]=xyz_list[ 3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);3746 xyz_front[pt1][0]=xyz_list[ 3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);3747 xyz_front[pt1][1]=xyz_list[ 3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);3748 xyz_front[pt1][2]=xyz_list[ 3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);3736 xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]); 3737 xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]); 3738 xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]); 3739 xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]); 3740 xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]); 3741 xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]); 3749 3742 } 3750 3743 else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 … … 3757 3750 } 3758 3751 3759 xyz_front[pt1][0]=xyz_list[ 3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);3760 xyz_front[pt1][1]=xyz_list[ 3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);3761 xyz_front[pt1][2]=xyz_list[ 3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);3762 xyz_front[pt2][0]=xyz_list[ 3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);3763 xyz_front[pt2][1]=xyz_list[ 3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);3764 xyz_front[pt2][2]=xyz_list[ 3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);3752 xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]); 3753 xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]); 3754 xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]); 3755 xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]); 3756 xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]); 3757 xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]); 3765 3758 } 3766 3759 else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 … … 3773 3766 } 3774 3767 3775 xyz_front[pt2][0]=xyz_list[ 3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);3776 xyz_front[pt2][1]=xyz_list[ 3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);3777 xyz_front[pt2][2]=xyz_list[ 3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);3778 xyz_front[pt1][0]=xyz_list[ 3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);3779 xyz_front[pt1][1]=xyz_list[ 3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);3780 xyz_front[pt1][2]=xyz_list[ 3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);3768 xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]); 3769 xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]); 3770 xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]); 3771 xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]); 3772 xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]); 3773 xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]); 3781 3774 } 3782 3775 else{ … … 3806 3799 3807 3800 /*Start looping on Gaussian points*/ 3808 Gauss* gauss=this->NewGaussBase( xyz_list,&xyz_front[0][0],3);3801 Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3); 3809 3802 while(gauss->next()){ 3810 3803 thickness_input->GetInputValue(&thickness,gauss); … … 3819 3812 delete gauss; 3820 3813 return flux; 3821 3822 3814 } 3823 3815 /*}}}*/ … … 3836 3828 IssmDouble xyz_front[2][3]; 3837 3829 3838 IssmDouble *xyz_list = NULL;3839 this->GetVerticesCoordinates(&xyz_list);3830 IssmDouble xyz_list[NUMVERTICES][3]; 3831 ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); 3840 3832 3841 3833 /*Recover parameters and values*/ … … 3857 3849 pt1 = 1; pt2 = 0; 3858 3850 } 3859 xyz_front[pt2][0]=xyz_list[ 3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);3860 xyz_front[pt2][1]=xyz_list[ 3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);3861 xyz_front[pt2][2]=xyz_list[ 3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);3862 xyz_front[pt1][0]=xyz_list[ 3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);3863 xyz_front[pt1][1]=xyz_list[ 3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);3864 xyz_front[pt1][2]=xyz_list[ 3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);3851 xyz_front[pt2][0]=xyz_list[2][0]+s1*(xyz_list[1][0]-xyz_list[2][0]); 3852 xyz_front[pt2][1]=xyz_list[2][1]+s1*(xyz_list[1][1]-xyz_list[2][1]); 3853 xyz_front[pt2][2]=xyz_list[2][2]+s1*(xyz_list[1][2]-xyz_list[2][2]); 3854 xyz_front[pt1][0]=xyz_list[2][0]+s2*(xyz_list[0][0]-xyz_list[2][0]); 3855 xyz_front[pt1][1]=xyz_list[2][1]+s2*(xyz_list[0][1]-xyz_list[2][1]); 3856 xyz_front[pt1][2]=xyz_list[2][2]+s2*(xyz_list[0][2]-xyz_list[2][2]); 3865 3857 } 3866 3858 else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 … … 3873 3865 } 3874 3866 3875 xyz_front[pt1][0]=xyz_list[ 3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);3876 xyz_front[pt1][1]=xyz_list[ 3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);3877 xyz_front[pt1][2]=xyz_list[ 3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);3878 xyz_front[pt2][0]=xyz_list[ 3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);3879 xyz_front[pt2][1]=xyz_list[ 3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);3880 xyz_front[pt2][2]=xyz_list[ 3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+2]);3867 xyz_front[pt1][0]=xyz_list[0][0]+s1*(xyz_list[1][0]-xyz_list[0][0]); 3868 xyz_front[pt1][1]=xyz_list[0][1]+s1*(xyz_list[1][1]-xyz_list[0][1]); 3869 xyz_front[pt1][2]=xyz_list[0][2]+s1*(xyz_list[1][2]-xyz_list[0][2]); 3870 xyz_front[pt2][0]=xyz_list[0][0]+s2*(xyz_list[2][0]-xyz_list[0][0]); 3871 xyz_front[pt2][1]=xyz_list[0][1]+s2*(xyz_list[2][1]-xyz_list[0][1]); 3872 xyz_front[pt2][2]=xyz_list[0][2]+s2*(xyz_list[2][2]-xyz_list[0][2]); 3881 3873 } 3882 3874 else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 … … 3889 3881 } 3890 3882 3891 xyz_front[pt2][0]=xyz_list[ 3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);3892 xyz_front[pt2][1]=xyz_list[ 3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);3893 xyz_front[pt2][2]=xyz_list[ 3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);3894 xyz_front[pt1][0]=xyz_list[ 3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);3895 xyz_front[pt1][1]=xyz_list[ 3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);3896 xyz_front[pt1][2]=xyz_list[ 3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);3883 xyz_front[pt2][0]=xyz_list[1][0]+s1*(xyz_list[0][0]-xyz_list[1][0]); 3884 xyz_front[pt2][1]=xyz_list[1][1]+s1*(xyz_list[0][1]-xyz_list[1][1]); 3885 xyz_front[pt2][2]=xyz_list[1][2]+s1*(xyz_list[0][2]-xyz_list[1][2]); 3886 xyz_front[pt1][0]=xyz_list[1][0]+s2*(xyz_list[2][0]-xyz_list[1][0]); 3887 xyz_front[pt1][1]=xyz_list[1][1]+s2*(xyz_list[2][1]-xyz_list[1][1]); 3888 xyz_front[pt1][2]=xyz_list[1][2]+s2*(xyz_list[2][2]-xyz_list[1][2]); 3897 3889 } 3898 3890 else{ … … 3931 3923 3932 3924 /*Start looping on Gaussian points*/ 3933 Gauss* gauss=this->NewGaussBase( xyz_list,&xyz_front[0][0],3);3925 Gauss* gauss=this->NewGaussBase(&xyz_list[0][0],&xyz_front[0][0],3); 3934 3926 while(gauss->next()){ 3935 3927 thickness_input->GetInputValue(&thickness,gauss); … … 3950 3942 delete gauss; 3951 3943 return flux; 3952 3953 3944 } 3954 3945 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.