source: issm/oecreview/Archive/15392-16133/ISSM-15474-15475.diff@ 16134

Last change on this file since 16134 was 16134, checked in by Mathieu Morlighem, 12 years ago

Added Archive/15392-16133
File size: 4.7 KB
RevLine 
[16134]1Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15474)
4+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15475)
5@@ -6917,7 +6917,7 @@
6 GaussPenta *gauss=NULL;
7
8 /*Stabilization*/
9- bool stabilization = false;
10+ bool stabilization = true;
11 IssmDouble dbasis[3][6];
12 IssmDouble dmu[3];
13 IssmDouble mu;
14@@ -6943,6 +6943,20 @@
15 rigidity=material->GetB();
16 diameter=MinEdgeLength(xyz_list);
17
18+ if(stabilization){
19+ gauss=new GaussPenta();
20+ for(int iv=0;iv<6;iv++){
21+ gauss->GaussVertex(iv);
22+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
23+ for(i=0;i<6;i++){
24+ for(j=0;j<3;j++){
25+ dnodalbasis[i][iv][j] = dbasis[j][i];
26+ }
27+ }
28+ }
29+ delete gauss;
30+ }
31+
32 /* Start looping on the number of gaussian points: */
33 gauss=new GaussPenta(3,2);
34 for(int ig=gauss->begin();ig<gauss->end();ig++){
35@@ -6971,14 +6985,7 @@
36 if(stabilization){
37 GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
38 dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
39- mu = viscosity;
40- for(i=0;i<6;i++){
41- for(p=0;p<6;p++){
42- for(j=0;j<3;j++){
43- dnodalbasis[i][p][j] = dbasis[j][i];
44- }
45- }
46- }
47+ mu = 2.*viscosity;
48 for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
49 SU[p][i][j]=0.;
50 SW[p][i][j]=0.;
51@@ -7017,34 +7024,8 @@
52 }
53
54 }
55-
56-
57- //viscosity = 1000*rigidity;
58- //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
59- //if(id==24){
60- //printf("tau %g\n",1./3.*pow(diameter,2)/(8.*2.*viscosity));
61- //}
62- //GetBConduct(&B_stab[0][0],&xyz_list[0][0],gauss);
63-
64- //D_stab[0][0]=D_scalar_stab; D_stab[0][1]=0; D_stab[0][2]=0;
65- //D_stab[1][0]=0; D_stab[1][1]=D_scalar_stab; D_stab[1][2]=0;
66- //D_stab[2][0]=0; D_stab[2][1]=0; D_stab[2][2]=D_scalar_stab;
67-
68- //TripleMultiply(&B_stab[0][0],3,6,1,
69- // &D_stab[0][0],3,3,0,
70- // &B_stab[0][0],3,6,0,
71- // &Ke_temp_stab[0][0],0);
72-
73- //for(i=0;i<NUMVERTICES;i++) for(j=0;j<NUMVERTICES;j++) Ke->values[numdof*(i*4+3)+j*4+3]+=Ke_temp_stab[i][j];
74-
75 }
76
77-// if(id==24){
78-// _error_("");
79-// }
80-
81- //Ke->Echo();
82- //_error_("stop");
83 /*Transform Coordinate System*/
84 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
85 if(id==24){
86@@ -7943,8 +7924,22 @@
87 /*Find minimal length*/
88 diameter=MinEdgeLength(xyz_list);
89
90+ if(stabilization){
91+ gauss=new GaussPenta();
92+ for(int iv=0;iv<6;iv++){
93+ gauss->GaussVertex(iv);
94+ GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
95+ for(i=0;i<6;i++){
96+ for(j=0;j<3;j++){
97+ dnodalbasis[i][iv][j] = dbasis[j][i];
98+ }
99+ }
100+ }
101+ delete gauss;
102+ }
103+
104 /* Start looping on the number of gaussian points: */
105- gauss=new GaussPenta(2,3);
106+ gauss=new GaussPenta(3,2);
107 for(int ig=gauss->begin();ig<gauss->end();ig++){
108
109 gauss->GaussPoint(ig);
110@@ -7969,15 +7964,7 @@
111 if(stabilization){
112 GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
113 dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
114- mu = viscosity;
115- for(i=0;i<6;i++){
116- for(p=0;p<6;p++){
117- for(j=0;j<3;j++){
118- dnodalbasis[i][p][j] = dbasis[j][i];
119- }
120- }
121- }
122- //dNodalBasisdx(1:n,p,:) = dBasisdx(1:n,:)
123+ mu = 2.*viscosity;
124 for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
125 SW[p][i][j]=0.;
126 }
127@@ -7988,10 +7975,8 @@
128 SW[p][i][i] += -dmu[j]*dbasis[j][p];
129 SW[p][j][i] += -dmu[j]*dbasis[i][p];
130 for(ii=0;ii<6;ii++){
131- //SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
132- //SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
133- SW[p][i][i] += -mu*dnodalbasis[p][ii][j];
134- SW[p][j][i] += -mu*dnodalbasis[p][ii][i];
135+ SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
136+ SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
137 }
138 }
139 }
140@@ -8006,26 +7991,11 @@
141 }
142
143 }
144- ///*Add stabilization*/
145- //D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
146- //GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0],gauss);
147-
148- //for(i=0;i<NUMVERTICES;i++){
149- // Pe_gaussian[i*NDOF4+3]+=-rho_ice*gravity*D_scalar_stab*dh1dh6[2][i];
150- // Pe_gaussian[i*NDOF4+3]+=forcex*D_scalar_stab*dh1dh6[0][i];
151- // Pe_gaussian[i*NDOF4+3]+=forcey*D_scalar_stab*dh1dh6[1][i];
152- // Pe_gaussian[i*NDOF4+3]+=forcez*D_scalar_stab*dh1dh6[2][i];
153- //}
154 }
155
156 /*Transform coordinate system*/
157 TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
158
159- if(id==24){
160- printarray(pe->values,24,1);
161- _error_("STOP");
162- }
163-
164 /*Clean up and return*/
165 delete gauss;
166 return pe;
Note: See TracBrowser for help on using the repository browser.