source: issm/oecreview/Archive/16554-17801/ISSM-17547-17548.diff

Last change on this file was 17802, checked in by Mathieu Morlighem, 11 years ago

Added archives
File size: 7.9 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17547)
4+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17548)
5@@ -4257,6 +4257,7 @@
6 IssmDouble dvx[3],dvy[3],dvz[3],B,n;
7 IssmDouble *xyz_list = NULL;
8 IssmDouble Jdet,r;
9+ Gauss* gauss = NULL;
10
11 parameters->FindParam(&meshtype,MeshTypeEnum);
12 switch(meshtype){
13@@ -4308,7 +4309,7 @@
14 sigmapyz_input=element->GetInput(DeviatoricStressyzEnum); _assert_(sigmapyz_input);
15 }
16
17- Gauss* gauss=element->NewGauss(5);
18+ gauss=element->NewGauss(5);
19 for(int ig=gauss->begin();ig<gauss->end();ig++){
20 gauss->GaussPoint(ig);
21 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
22@@ -4373,30 +4374,128 @@
23 Ke,1);
24
25 /*Create Right hand sides*/
26- for(int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += sigmapxx*tbasis[ii]*gauss->weight*Jdet;
27- for(int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += sigmapyy*tbasis[ii]*gauss->weight*Jdet;
28- for(int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += sigmapxy*tbasis[ii]*gauss->weight*Jdet;
29+ for(int ii=0;ii<tnumnodes;ii++) pe_xx[ii] += (r*epsxx+sigmapxx)*tbasis[ii]*gauss->weight*Jdet;
30+ for(int ii=0;ii<tnumnodes;ii++) pe_yy[ii] += (r*epsyy+sigmapyy)*tbasis[ii]*gauss->weight*Jdet;
31+ for(int ii=0;ii<tnumnodes;ii++) pe_xy[ii] += (r*epsxy+sigmapxy)*tbasis[ii]*gauss->weight*Jdet;
32 if(dim==3){
33- for(int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += sigmapzz*tbasis[ii]*gauss->weight*Jdet;
34- for(int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += sigmapxz*tbasis[ii]*gauss->weight*Jdet;
35- for(int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += sigmapyz*tbasis[ii]*gauss->weight*Jdet;
36+ for(int ii=0;ii<tnumnodes;ii++) pe_zz[ii] += (r*epszz+sigmapzz)*tbasis[ii]*gauss->weight*Jdet;
37+ for(int ii=0;ii<tnumnodes;ii++) pe_xz[ii] += (r*epsxz+sigmapxz)*tbasis[ii]*gauss->weight*Jdet;
38+ for(int ii=0;ii<tnumnodes;ii++) pe_yz[ii] += (r*epsyz+sigmapyz)*tbasis[ii]*gauss->weight*Jdet;
39 }
40 }
41
42 /*Solve the systems*/
43- _error_("S");
44+ IssmDouble* d_xx = xNew<IssmDouble>(tnumnodes);
45+ IssmDouble* d_yy = xNew<IssmDouble>(tnumnodes);
46+ IssmDouble* d_xy = xNew<IssmDouble>(tnumnodes);
47+ IssmDouble* d_zz = NULL;
48+ IssmDouble* d_xz = NULL;
49+ IssmDouble* d_yz = NULL;
50+ if(dim==2){
51+ _assert_(tnumnodes==3);
52+ Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
53+ Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
54+ Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
55+ element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
56+ element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
57+ element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
58+ }
59+ else{
60+ _assert_(tnumnodes==4);
61+ Matrix3x3Solve(&d_xx[0],Ke,pe_xx);
62+ Matrix3x3Solve(&d_yy[0],Ke,pe_yy);
63+ Matrix3x3Solve(&d_xy[0],Ke,pe_xy);
64+ Matrix3x3Solve(&d_zz[0],Ke,pe_zz);
65+ Matrix3x3Solve(&d_xz[0],Ke,pe_xz);
66+ Matrix3x3Solve(&d_yz[0],Ke,pe_yz);
67+ element->AddInput(StrainRatexxEnum,d_xx,P1DGEnum);
68+ element->AddInput(StrainRateyyEnum,d_yy,P1DGEnum);
69+ element->AddInput(StrainRatexyEnum,d_xy,P1DGEnum);
70+ element->AddInput(StrainRatezzEnum,d_zz,P1DGEnum);
71+ element->AddInput(StrainRatexzEnum,d_xz,P1DGEnum);
72+ element->AddInput(StrainRateyzEnum,d_yz,P1DGEnum);
73+ }
74
75+ /*Update tau accordingly*/
76+ IssmDouble* tau_xx = xNew<IssmDouble>(tnumnodes);
77+ IssmDouble* tau_yy = xNew<IssmDouble>(tnumnodes);
78+ IssmDouble* tau_xy = xNew<IssmDouble>(tnumnodes);
79+ IssmDouble* tau_zz = NULL;
80+ IssmDouble* tau_xz = NULL;
81+ IssmDouble* tau_yz = NULL;
82+ Input* epsxx_input=element->GetInput(StrainRatexxEnum); _assert_(epsxx_input);
83+ Input* epsyy_input=element->GetInput(StrainRateyyEnum); _assert_(epsyy_input);
84+ Input* epsxy_input=element->GetInput(StrainRatexyEnum); _assert_(epsxy_input);
85+ Input* epszz_input=NULL; Input* epsxz_input=NULL; Input* epsyz_input=NULL;
86+ if(dim==3){
87+ epszz_input=element->GetInput(StrainRatezzEnum); _assert_(epszz_input);
88+ epsxz_input=element->GetInput(StrainRatexzEnum); _assert_(epsxz_input);
89+ epsyz_input=element->GetInput(StrainRateyzEnum); _assert_(epsyz_input);
90+ }
91+ Gauss* gauss = element->NewGauss();
92+ for(int ig=0;ig<tnumnodes;ig++){
93+ gauss->GaussNode(P1DGEnum,ig);
94+
95+ /*Get D(u)*/
96+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
97+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
98+ if(dim==3){
99+ vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
100+ }
101+ epsxx = dvx[0];
102+ epsyy = dvy[1];
103+ epsxy = 0.5*(dvx[1] + dvy[0]);
104+ if(dim==3){
105+ epszz = dvz[2];
106+ epsxz = 0.5*(dvx[2] + dvz[0]);
107+ epsyz = 0.5*(dvy[2] + dvz[1]);
108+ }
109+
110+ /*Get tau^(n-1) from inputs*/
111+ sigmapxx_input->GetInputValue(&sigmapxx,gauss);
112+ sigmapyy_input->GetInputValue(&sigmapyy,gauss);
113+ sigmapxy_input->GetInputValue(&sigmapxy,gauss);
114+ if(dim==3){
115+ sigmapzz_input->GetInputValue(&sigmapzz,gauss);
116+ sigmapxz_input->GetInputValue(&sigmapxz,gauss);
117+ sigmapyz_input->GetInputValue(&sigmapyz,gauss);
118+ }
119+
120+ /*Get d and update tau accordingly*/
121+ tau_xx[ig] = sigmapxx + r*(epsxx - d_xx[ig]);
122+ tau_yy[ig] = sigmapyy + r*(epsyy - d_yy[ig]);
123+ tau_xy[ig] = sigmapxy + r*(epsxy - d_xy[ig]);
124+ if(dim==3){
125+ tau_zz[ig] = sigmapzz + r*(epszz - d_zz[ig]);
126+ tau_xz[ig] = sigmapxz + r*(epsxz - d_xz[ig]);
127+ tau_yz[ig] = sigmapyz + r*(epsyz - d_yz[ig]);
128+ }
129+ }
130+
131+ /*Add inputs*/
132+ element->AddInput(StrainRatexxEnum,tau_xx,P1DGEnum);
133+ element->AddInput(StrainRateyyEnum,tau_yy,P1DGEnum);
134+ element->AddInput(StrainRatexyEnum,tau_xy,P1DGEnum);
135+ if(dim==3){
136+ element->AddInput(StrainRatezzEnum,tau_zz,P1DGEnum);
137+ element->AddInput(StrainRatexzEnum,tau_xz,P1DGEnum);
138+ element->AddInput(StrainRateyzEnum,tau_yz,P1DGEnum);
139+ }
140+
141 /*Clean up and */
142+ xDelete<IssmDouble>(xyz_list);
143 xDelete<IssmDouble>(tbasis);
144 xDelete<IssmDouble>(Ke);
145- xDelete<IssmDouble>(pe_xx);
146- xDelete<IssmDouble>(pe_yy);
147- xDelete<IssmDouble>(pe_zz);
148- xDelete<IssmDouble>(pe_xy);
149- xDelete<IssmDouble>(pe_xz);
150- xDelete<IssmDouble>(pe_yz);
151+ xDelete<IssmDouble>(pe_xx); xDelete<IssmDouble>(d_xx); xDelete<IssmDouble>(tau_xx);
152+ xDelete<IssmDouble>(pe_yy); xDelete<IssmDouble>(d_yy); xDelete<IssmDouble>(tau_yy);
153+ xDelete<IssmDouble>(pe_zz); xDelete<IssmDouble>(d_zz); xDelete<IssmDouble>(tau_zz);
154+ xDelete<IssmDouble>(pe_xy); xDelete<IssmDouble>(d_xy); xDelete<IssmDouble>(tau_xy);
155+ xDelete<IssmDouble>(pe_xz); xDelete<IssmDouble>(d_xz); xDelete<IssmDouble>(tau_xz);
156+ xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz); xDelete<IssmDouble>(tau_yz);
157 }
158
159+ delete gauss;
160+
161 }/*}}}*/
162
163 /*Coupling (Tiling)*/
164Index: ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
165===================================================================
166--- ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 17547)
167+++ ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 17548)
168@@ -366,3 +366,16 @@
169 Ainv[8]=(A[0]*A[4]-A[1]*A[3])/det; /* = (a*e-b*d)/det */
170
171 }/*}}}*/
172+/*FUNCTION Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B) {{{*/
173+void Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B){
174+
175+ IssmDouble Ainv[3][3];
176+
177+ Matrix3x3Invert(&Ainv[0][0],A);
178+ for(int i=0;i<3;i++) X[i]=Ainv[i][0]*B[0] + Ainv[i][1]*B[1] + Ainv[i][2]*B[2];
179+
180+}/*}}}*/
181+/*FUNCTION Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble *B) {{{*/
182+void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble *B){
183+ _error_("STOP");
184+}/*}}}*/
185Index: ../trunk-jpl/src/c/shared/Matrix/matrix.h
186===================================================================
187--- ../trunk-jpl/src/c/shared/Matrix/matrix.h (revision 17547)
188+++ ../trunk-jpl/src/c/shared/Matrix/matrix.h (revision 17548)
189@@ -14,5 +14,7 @@
190 void Matrix2x2Determinant(IssmDouble* Adet,IssmDouble* A);
191 void Matrix3x3Invert(IssmDouble* Ainv, IssmDouble* A);
192 void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A);
193+void Matrix3x3Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B);
194+void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B);
195
196 #endif //ifndef _MATRIXUTILS_H_
Note: See TracBrowser for help on using the repository browser.