Changeset 9217
- Timestamp:
- 08/09/11 11:30:31 (14 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r9206 r9217 783 783 ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum); 784 784 ElementMatrix* Ke2=new ElementMatrix(this->nodes ,NUMVERTICES,this->parameters,PattynApproximationEnum); 785 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);785 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); 786 786 delete Ke1; delete Ke2; 787 787 … … 1451 1451 double B[5][numdof]; 1452 1452 double Bprime[5][numdof]; 1453 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.1454 1453 Tria* tria=NULL; 1455 1454 GaussPenta *gauss=NULL; … … 1489 1488 &D[0][0],5,5,0, 1490 1489 &Bprime[0][0],5,numdof,0, 1491 &Ke_gg_gaussian[0][0],0); 1492 1493 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j]; 1490 &Ke->values[0],1); 1494 1491 } 1495 1492 … … 1517 1514 double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag 1518 1515 double DL_scalar; 1519 double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag1520 1516 Friction *friction = NULL; 1521 1517 GaussPenta *gauss=NULL; … … 1565 1561 &DL[0][0],2,2,0, 1566 1562 &L[0][0],2,numdof,0, 1567 &Ke_gg_gaussian[0][0],0); 1568 1569 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j]; 1563 &Ke->values[0],1); 1570 1564 } 1571 1565 … … 1780 1774 double Bprime[NDOF1][NUMVERTICES]; 1781 1775 double DL_scalar; 1782 double Ke_gg[numdof][numdof]={0.0};1783 1776 GaussPenta *gauss=NULL; 1784 1777 … … 1804 1797 &DL_scalar,1,1,0, 1805 1798 &Bprime[0][0],1,NUMVERTICES,0, 1806 &Ke_gg[0][0],0); 1807 1808 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg[i][j]; 1799 &Ke->values[0],1); 1809 1800 } 1810 1801 … … 1830 1821 double basis[NUMVERTICES]; 1831 1822 GaussPenta *gauss=NULL; 1832 double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.1833 1823 1834 1824 /*Initialize Element matrix*/ … … 1849 1839 GetNodalFunctionsP1(&basis[0], gauss); 1850 1840 1851 /**********************Do not forget the sign**********************************/ 1852 DL_scalar=- gauss->weight*Jdet2d*surface_normal[2]; 1853 /******************************************************************************/ 1841 DL_scalar= - gauss->weight*Jdet2d*surface_normal[2]; 1854 1842 1855 1843 TripleMultiply( basis,1,numdof,1, 1856 1844 &DL_scalar,1,1,0, 1857 1845 basis,1,numdof,0, 1858 &Ke_g[0][0],0); 1859 1860 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j]; 1846 &Ke->values[0],1); 1861 1847 } 1862 1848 … … 1908 1894 double D[3][3]; 1909 1895 double K[2][2]={0.0}; 1910 double Ke_gaussian_conduct[numdof][numdof];1911 double Ke_gaussian_advec[numdof][numdof];1912 double Ke_gaussian_artdiff[numdof][numdof];1913 double Ke_gaussian_transient[numdof][numdof];1914 1896 Tria* tria=NULL; 1915 1897 GaussPenta *gauss=NULL; … … 1965 1947 &D[0][0],3,3,0, 1966 1948 &B_conduct[0][0],3,numdof,0, 1967 &Ke _gaussian_conduct[0][0],0);1949 &Ke->values[0],1); 1968 1950 1969 1951 /*Advection: */ … … 1989 1971 &D[0][0],3,3,0, 1990 1972 &Bprime_advec[0][0],3,numdof,0, 1991 &Ke _gaussian_advec[0][0],0);1973 &Ke->values[0],1); 1992 1974 1993 1975 /*Transient: */ … … 2001 1983 &D_scalar_trans,1,1,0, 2002 1984 &L[0],1,numdof,0, 2003 &Ke_gaussian_transient[0][0],0); 2004 } 2005 else{ 2006 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_transient[i][j]=0; 1985 &Ke->values[0],1); 2007 1986 } 2008 1987 … … 2021 2000 &K[0][0],2,2,0, 2022 2001 &B_artdiff[0][0],2,numdof,0, 2023 &Ke _gaussian_artdiff[0][0],0);2002 &Ke->values[0],1); 2024 2003 } 2025 2004 else if(artdiff==2){ … … 2031 2010 for(i=0;i<numdof;i++){ 2032 2011 for(j=0;j<numdof;j++){ 2033 Ke _gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*2012 Ke->values[i*numdof+j]=tau_parameter*D_scalar_advec* 2034 2013 ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]); 2035 2014 } … … 2038 2017 for(i=0;i<numdof;i++){ 2039 2018 for(j=0;j<numdof;j++){ 2040 Ke _gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);2019 Ke->values[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]); 2041 2020 } 2042 2021 } 2043 2022 } 2044 2023 } 2045 else{2046 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_artdiff[i][j]=0;2047 }2048 2049 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];2050 2024 } 2051 2025 … … 2070 2044 double basis[NUMVERTICES]; 2071 2045 double D_scalar; 2072 double Ke_gaussian[numdof][numdof]={0.0};2073 2046 GaussPenta *gauss=NULL; 2074 2047 … … 2102 2075 &D_scalar,1,1,0, 2103 2076 &basis[0],1,numdof,0, 2104 &Ke_gaussian[0][0],0); 2105 2106 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j]; 2077 &Ke->values[0],1); 2107 2078 } 2108 2079 … … 2198 2169 double D[3][3]; 2199 2170 double K[2][2]={0.0}; 2200 double Ke_gaussian_conduct[numdof][numdof];2201 double Ke_gaussian_advec[numdof][numdof];2202 double Ke_gaussian_artdiff[numdof][numdof];2203 double Ke_gaussian_transient[numdof][numdof];2204 2171 Tria* tria=NULL; 2205 2172 GaussPenta *gauss=NULL; … … 2248 2215 &D[0][0],3,3,0, 2249 2216 &B_conduct[0][0],3,numdof,0, 2250 &Ke _gaussian_conduct[0][0],0);2217 &Ke->values[0],1); 2251 2218 2252 2219 /*Advection: */ … … 2272 2239 &D[0][0],3,3,0, 2273 2240 &Bprime_advec[0][0],3,numdof,0, 2274 &Ke _gaussian_advec[0][0],0);2241 &Ke->values[0],1); 2275 2242 2276 2243 /*Transient: */ … … 2284 2251 &D_scalar_trans,1,1,0, 2285 2252 &L[0],1,numdof,0, 2286 &Ke_gaussian_transient[0][0],0); 2287 } 2288 else{ 2289 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_transient[i][j]=0; 2253 &Ke->values[0],1); 2290 2254 } 2291 2255 … … 2304 2268 &K[0][0],2,2,0, 2305 2269 &B_artdiff[0][0],2,numdof,0, 2306 &Ke _gaussian_artdiff[0][0],0);2270 &Ke->values[0],1); 2307 2271 } 2308 2272 else if(artdiff==2){ … … 2314 2278 for(i=0;i<numdof;i++){ 2315 2279 for(j=0;j<numdof;j++){ 2316 Ke _gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*2280 Ke->values[i*numdof+j]=tau_parameter*D_scalar_advec* 2317 2281 ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]); 2318 2282 } … … 2321 2285 for(i=0;i<numdof;i++){ 2322 2286 for(j=0;j<numdof;j++){ 2323 Ke _gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);2287 Ke->values[i*numdof+j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]); 2324 2288 } 2325 2289 } 2326 2290 } 2327 2291 } 2328 else{2329 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gaussian_artdiff[i][j]=0;2330 }2331 2332 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j];2333 2292 } 2334 2293 … … 2354 2313 double basis[NUMVERTICES]; 2355 2314 double D_scalar; 2356 double Ke_gaussian[numdof][numdof]={0.0};2357 2315 GaussPenta *gauss=NULL; 2358 2316 … … 2386 2344 &D_scalar,1,1,0, 2387 2345 &basis[0],1,numdof,0, 2388 &Ke_gaussian[0][0],0); 2389 2390 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gaussian[i][j]; 2346 &Ke->values[0],1); 2391 2347 } 2392 2348 -
issm/trunk/src/c/objects/Elements/Tria.cpp
r9206 r9217 484 484 double DLprime[2][2] = {0.0}; 485 485 double DL_scalar; 486 double Ke_gg_gaussian[numdof][numdof] = {0.0};487 double Ke_gg_thickness1[numdof][numdof] = {0.0};488 double Ke_gg_thickness2[numdof][numdof] = {0.0};489 486 GaussTria *gauss = NULL; 490 487 … … 548 545 &DL[0][0],2,2,0, 549 546 &B[0][0],2,numdof,0, 550 &Ke _gg_thickness1[0][0],0);547 &Ke->values[0],1); 551 548 552 549 TripleMultiply( &B[0][0],2,numdof,1, 553 550 &DLprime[0][0],2,2,0, 554 551 &Bprime[0][0],2,numdof,0, 555 &Ke_gg_thickness2[0][0],0); 556 557 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness1[i][j]; 558 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_thickness2[i][j]; 552 &Ke->values[0],1); 559 553 560 554 if(artdiff){ … … 565 559 &KDL[0][0],2,2,0, 566 560 &Bprime[0][0],2,numdof,0, 567 &Ke_gg_gaussian[0][0],0); 568 569 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j]; 561 &Ke->values[0],1); 570 562 } 571 563 } … … 590 582 double DL[2][2]={0.0}; 591 583 double DL_scalar; 592 double Ke_gg[numdof][numdof]={0.0};593 584 GaussTria *gauss=NULL; 594 585 … … 623 614 &DL[0][0],2,2,0, 624 615 &Bprime[0][0],2,numdof,0, 625 &Ke_gg[0][0],0); 626 627 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg[i][j]; 616 &Ke->values[0],1); 628 617 } 629 618 … … 655 644 double DLprime[2][2]={0.0}; 656 645 double DL_scalar; 657 double Ke_gg_gaussian[numdof][numdof] = {0.0};658 double Ke_gg_velocities[numdof][numdof] = {0.0};659 646 GaussTria *gauss=NULL; 660 647 … … 734 721 &DLprime[0][0],2,2,0, 735 722 &Bprime[0][0],2,numdof,0, 736 &Ke_gg_velocities[0][0],0); 737 738 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_velocities[i][j]; 723 &Ke->values[0],1); 739 724 740 725 if(artdiff){ … … 745 730 &KDL[0][0],2,2,0, 746 731 &Bprime[0][0],2,numdof,0, 747 &Ke_gg_gaussian[0][0],0); 748 749 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j]; 732 &Ke->values[0],1); 750 733 } 751 734 } … … 786 769 double D[3][3] = {0.0}; 787 770 double D_scalar; 788 double Ke_g[numdof][numdof];789 771 GaussTria *gauss = NULL; 790 772 … … 824 806 &D[0][0],3,3,0, 825 807 &Bprime[0][0],3,numdof,0, 826 &Ke_g[0][0],0); 827 828 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j]; 808 &Ke->values[0],1); 829 809 } 830 810 … … 849 829 double L[2][numdof]; 850 830 double DL[2][2] = {{ 0,0 },{0,0}}; 851 double Ke_g[numdof][numdof];852 831 double DL_scalar; 853 832 double slope[2] = {0.0,0.0}; … … 894 873 &DL[0][0],2,2,0, 895 874 &L[0][0],2,numdof,0, 896 &Ke_g[0][0],0); 897 898 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j]; 875 &Ke->values[0],1); 899 876 } 900 877 … … 985 962 double dLy[NUMVERTICES]; 986 963 double K[2]; 987 double Ke_gg_gaussian1[numdof][numdof] ={0.0};988 double Ke_gg_gaussian2[numdof][numdof] ={0.0};989 double Ke_gg_gaussian3[numdof][numdof] ={0.0};990 964 double hydro_p,hydro_q,hydro_CR,hydro_n,rho_water,g,gamma; 991 965 double DL_scalar1; … … 1019 993 GetNodalFunctions(&L[0],gauss); 1020 994 GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss); 1021 //printf("NUMVERTICES=%d, Jdettria=%f\n",NUMVERTICES,Jdettria); //BK1022 995 for(i=0;i<NUMVERTICES;i++)dLx[i]=dL[0][i]; 1023 996 for(i=0;i<NUMVERTICES;i++)dLy[i]=dL[1][i]; 1024 /* for(i=0;i<NUMVERTICES;i++)1025 { dLx[i]=dL[0][i];1026 printf("dLx[%d]=%f, L[%d]=%f\n",i,dLx[i]),i,L[i];}1027 for(i=0;i<NUMVERTICES;i++)1028 { dLy[i]=dL[1][i];1029 printf("dLy[%d]=%f\n",i,dLy[i]);} */1030 997 1031 998 /*Get K parameter: */ 1032 999 GetHydrologyK(&K[0],&xyz_list[0][0],gauss); 1033 //printf("K[0]=%g,K[1]=%g\n",K[0],K[1]);1034 1000 1035 1001 if(dt){ … … 1038 1004 &DL_scalar1,1,1,0, 1039 1005 &L[0],1,numdof,0, 1040 &Ke _gg_gaussian1[0][0],0);1006 &Ke->values[0],1); 1041 1007 } 1042 1008 1043 1009 if(dt){ 1044 DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; //don't forget the -1045 DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; //don't forget the -1010 DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; 1011 DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; 1046 1012 1047 1013 } 1048 1014 else{ 1049 DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; //don't forget the -1050 DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; //don't forget the -1015 DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; 1016 DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; 1051 1017 } 1052 1018 1053 TripleMultiply( &dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian2[0][0],0); 1054 TripleMultiply( &dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian3[0][0],0); 1055 1056 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian1[i][j]; 1057 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian2[i][j]; 1058 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian3[i][j]; 1059 1019 TripleMultiply(&dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1); 1020 TripleMultiply(&dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1); 1060 1021 } 1061 1022 … … 1210 1171 double DLprime[2][2]={0.0}; 1211 1172 double DL_scalar; 1212 double Ke_gg1[numdof][numdof]={0.0};1213 double Ke_gg2[numdof][numdof]={0.0};1214 1173 GaussTria *gauss=NULL; 1215 1174 … … 1249 1208 &DL_scalar,1,1,0, 1250 1209 &L[0],1,numdof,0, 1251 &Ke _gg1[0][0],0);1210 &Ke->values[0],1); 1252 1211 1253 1212 /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/ … … 1263 1222 &DLprime[0][0],2,2,0, 1264 1223 &Bprime[0][0],2,numdof,0, 1265 &Ke_gg2[0][0],0); 1266 1267 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg1[i][j]; 1268 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg2[i][j]; 1224 &Ke->values[0],1); 1269 1225 } 1270 1226 … … 1285 1241 double xyz_list[NUMVERTICES][3]; 1286 1242 double L[1][3]; 1287 double Ke_g[numdof][numdof];1288 1243 GaussTria *gauss = NULL; 1289 1244 … … 1307 1262 &DL_scalar,1,1,0, 1308 1263 &L[0][0],1,3,0, 1309 &Ke_g[0][0],0); 1310 1311 for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_g[i][j]; 1264 &Ke->values[0],1); 1312 1265 } 1313 1266 … … 2868 2821 K[1]=pow(w,hydro_p)*pow((rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn/w * (dsdy - dbdy ) * surface_slope,hydro_q); 2869 2822 2870 //printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w); //bk2871 //if (w<0) {printf("negative w!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n");2872 //printf("K[0]=%g,K[1]=%g,w=%g\n",K[0],K[1],w);}2873 2823 } 2874 2824 /*}}}*/
Note:
See TracChangeset
for help on using the changeset viewer.