source: issm/oecreview/Archive/11431-11450/ISSM-11447-11448.diff@ 11515

Last change on this file since 11515 was 11515, checked in by Eric.Larour, 13 years ago

Oecreview up to 11509
File size: 14.7 KB
RevLine 
[11515]1Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.cpp
2===================================================================
3--- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11447)
4+++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 11448)
5@@ -335,9 +335,12 @@
6 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
7 switch(analysis_type){
8 #ifdef _HAVE_DIAGNOSTIC_
9- case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
10+ case DiagnosticHorizAnalysisEnum:
11 Ke=CreateKMatrixDiagnosticMacAyeal();
12 break;
13+ case AdjointHorizAnalysisEnum:
14+ Ke=CreateKMatrixAdjointMacAyeal();
15+ break;
16 case DiagnosticHutterAnalysisEnum:
17 Ke=CreateKMatrixDiagnosticHutter();
18 break;
19@@ -3062,7 +3065,7 @@
20 return pe;
21 }
22 /*}}}*/
23-/*FUNCTION Tria::CreateJacobianDiagnosticPattyn{{{1*/
24+/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{1*/
25 ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
26
27 /*Constants*/
28@@ -3863,7 +3866,7 @@
29 rheologyb_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
30
31 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
32- //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
33+ Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
34 }
35
36 /*Clean up and return*/
37@@ -4735,7 +4738,7 @@
38 drag_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
39
40 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
41- //Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
42+ Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
43 }
44
45 /*Clean up and return*/
46@@ -4765,6 +4768,76 @@
47 return Ke;
48 }
49 /*}}}*/
50+/*FUNCTION Tria::CreateKMatrixAdjointMacAyeal{{{1*/
51+ElementMatrix* Tria::CreateKMatrixAdjointMacAyeal(void){
52+
53+ /*Constants*/
54+ const int numdof=NDOF2*NUMVERTICES;
55+
56+ /*Intermediaries */
57+ int i,j,ig;
58+ bool incomplete_adjoint;
59+ double xyz_list[NUMVERTICES][3];
60+ double Jdet,thickness;
61+ double eps1dotdphii,eps1dotdphij;
62+ double eps2dotdphii,eps2dotdphij;
63+ double mu_prime;
64+ double epsilon[3];/* epsilon=[exx,eyy,exy];*/
65+ double eps1[2],eps2[2];
66+ double phi[NUMVERTICES];
67+ double dphi[2][NUMVERTICES];
68+ GaussTria *gauss=NULL;
69+
70+ /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
71+ parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
72+ ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
73+ if(incomplete_adjoint) return Ke;
74+
75+ /*Retrieve all inputs and parameters*/
76+ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
77+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
78+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
79+ Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
80+
81+ /* Start looping on the number of gaussian points: */
82+ gauss=new GaussTria(2);
83+ for (ig=gauss->begin();ig<gauss->end();ig++){
84+
85+ gauss->GaussPoint(ig);
86+
87+ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
88+ GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
89+
90+ thickness_input->GetInputValue(&thickness, gauss);
91+ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
92+ matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
93+ eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
94+ eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
95+
96+ for(i=0;i<3;i++){
97+ for(j=0;j<3;j++){
98+ eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
99+ eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
100+ eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
101+ eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
102+
103+ Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
104+ Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
105+ Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
106+ Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
107+ }
108+ }
109+ }
110+
111+ /*Transform Coordinate System*/
112+ TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
113+
114+ /*Clean up and return*/
115+ delete gauss;
116+ //Ke->Transpose();
117+ return Ke;
118+}
119+/*}}}*/
120 /*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/
121 void Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){
122
123Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.h
124===================================================================
125--- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11447)
126+++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Tria.h (revision 11448)
127@@ -216,6 +216,7 @@
128
129 #ifdef _HAVE_CONTROL_
130 ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
131+ ElementMatrix* CreateKMatrixAdjointMacAyeal(void);
132 ElementVector* CreatePVectorAdjointHoriz(void);
133 ElementVector* CreatePVectorAdjointStokes(void);
134 ElementVector* CreatePVectorAdjointBalancethickness(void);
135Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.cpp
136===================================================================
137--- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11447)
138+++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 11448)
139@@ -584,9 +584,12 @@
140 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
141 switch(analysis_type){
142 #ifdef _HAVE_DIAGNOSTIC_
143- case DiagnosticHorizAnalysisEnum: case AdjointHorizAnalysisEnum:
144+ case DiagnosticHorizAnalysisEnum:
145 Ke=CreateKMatrixDiagnosticHoriz(); De=CreateDVectorDiagnosticHoriz();
146 break;
147+ case AdjointHorizAnalysisEnum:
148+ Ke=CreateKMatrixAdjointHoriz();
149+ break;
150 case DiagnosticHutterAnalysisEnum:
151 Ke=CreateKMatrixDiagnosticHutter();
152 break;
153@@ -4442,6 +4445,196 @@
154 ((ControlInput*)input)->SetGradient(grad_input);
155
156 }/*}}}*/
157+/*FUNCTION Penta::CreateKMatrixAdjointHoriz{{{1*/
158+ElementMatrix* Penta::CreateKMatrixAdjointHoriz(void){
159+
160+ int approximation;
161+ inputs->GetInputValue(&approximation,ApproximationEnum);
162+
163+ switch(approximation){
164+ case MacAyealApproximationEnum:
165+ return CreateKMatrixAdjointMacAyeal2d();
166+ case PattynApproximationEnum:
167+ return CreateKMatrixAdjointPattyn();
168+ case StokesApproximationEnum:
169+ return CreateKMatrixAdjointStokes();
170+ case NoneApproximationEnum:
171+ return NULL;
172+ default:
173+ _error_("Approximation %s not supported yet",EnumToStringx(approximation));
174+ }
175+}
176+/*}}}*/
177+/*FUNCTION Penta::CreateKMatrixAdjointMacAyeal2d{{{1*/
178+ElementMatrix* Penta::CreateKMatrixAdjointMacAyeal2d(void){
179+
180+ /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
181+ bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
182+ the stiffness matrix. */
183+ if (!IsOnBed()) return NULL;
184+
185+ /*Depth Averaging B*/
186+ this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
187+
188+ /*Call Tria function*/
189+ Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
190+ ElementMatrix* Ke=tria->CreateKMatrixAdjointMacAyeal();
191+ delete tria->matice; delete tria;
192+
193+ /*Delete B averaged*/
194+ this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
195+
196+ /*clean up and return*/
197+ return Ke;
198+}
199+/*}}}*/
200+/*FUNCTION Penta::CreateKMatrixAdjointPattyn{{{1*/
201+ElementMatrix* Penta::CreateKMatrixAdjointPattyn(void){
202+
203+ /*Constants*/
204+ const int numdof=NDOF2*NUMVERTICES;
205+
206+ /*Intermediaries */
207+ int i,j,ig;
208+ bool incomplete_adjoint;
209+ double xyz_list[NUMVERTICES][3];
210+ double Jdet;
211+ double eps1dotdphii,eps1dotdphij;
212+ double eps2dotdphii,eps2dotdphij;
213+ double mu_prime;
214+ double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
215+ double eps1[3],eps2[3];
216+ double phi[NUMVERTICES];
217+ double dphi[3][NUMVERTICES];
218+ GaussPenta *gauss=NULL;
219+
220+ /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
221+ parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
222+ ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
223+ if(incomplete_adjoint) return Ke;
224+
225+ /*Retrieve all inputs and parameters*/
226+ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
227+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
228+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
229+
230+ /* Start looping on the number of gaussian points: */
231+ gauss=new GaussPenta(5,5);
232+ for (ig=gauss->begin();ig<gauss->end();ig++){
233+
234+ gauss->GaussPoint(ig);
235+
236+ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
237+ GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
238+
239+ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
240+ matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
241+ eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
242+ eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
243+ eps1[2]=epsilon[3]; eps2[2]=epsilon[4];
244+
245+ for(i=0;i<6;i++){
246+ for(j=0;j<6;j++){
247+ eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
248+ eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
249+ eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
250+ eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
251+
252+ Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
253+ Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
254+ Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
255+ Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
256+ }
257+ }
258+ }
259+
260+ /*Transform Coordinate System*/
261+ TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
262+
263+ /*Clean up and return*/
264+ delete gauss;
265+ return Ke;
266+}
267+/*}}}*/
268+/*FUNCTION Penta::CreateKMatrixAdjointStokes{{{1*/
269+ElementMatrix* Penta::CreateKMatrixAdjointStokes(void){
270+
271+ /*Constants*/
272+ const int numdof=NDOF4*NUMVERTICES;
273+
274+ /*Intermediaries */
275+ int i,j,ig;
276+ bool incomplete_adjoint;
277+ double xyz_list[NUMVERTICES][3];
278+ double Jdet;
279+ double eps1dotdphii,eps1dotdphij;
280+ double eps2dotdphii,eps2dotdphij;
281+ double eps3dotdphii,eps3dotdphij;
282+ double mu_prime;
283+ double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
284+ double eps1[3],eps2[3],eps3[3];
285+ double phi[NUMVERTICES];
286+ double dphi[3][NUMVERTICES];
287+ GaussPenta *gauss=NULL;
288+
289+ /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
290+ parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
291+ ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
292+ if(incomplete_adjoint) return Ke;
293+
294+ /*Retrieve all inputs and parameters*/
295+ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
296+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
297+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
298+ Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
299+
300+ /* Start looping on the number of gaussian points: */
301+ gauss=new GaussPenta(5,5);
302+ for (ig=gauss->begin();ig<gauss->end();ig++){
303+
304+ gauss->GaussPoint(ig);
305+
306+ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
307+ GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
308+
309+ this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
310+ matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
311+ eps1[0]=epsilon[0]; eps2[0]=epsilon[2]; eps3[0]=epsilon[3];
312+ eps1[1]=epsilon[2]; eps2[1]=epsilon[1]; eps3[1]=epsilon[4];
313+ eps1[2]=epsilon[3]; eps2[2]=epsilon[4]; eps3[2]= -epsilon[0] -epsilon[1];
314+
315+ for(i=0;i<6;i++){
316+ for(j=0;j<6;j++){
317+ eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
318+ eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
319+ eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
320+ eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
321+ eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
322+ eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
323+
324+ Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
325+ Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
326+ Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
327+
328+ Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
329+ Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
330+ Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
331+
332+ Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
333+ Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
334+ Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
335+ }
336+ }
337+ }
338+
339+ /*Transform Coordinate System*/
340+ TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
341+
342+ /*Clean up and return*/
343+ delete gauss;
344+ return Ke;
345+}
346+/*}}}*/
347 /*FUNCTION Penta::CreatePVectorAdjointHoriz{{{1*/
348 ElementVector* Penta::CreatePVectorAdjointHoriz(void){
349
350Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.h
351===================================================================
352--- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11447)
353+++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Elements/Penta.h (revision 11448)
354@@ -215,6 +215,7 @@
355 ElementMatrix* CreateKMatrixCouplingMacAyealStokesFriction(void);
356 ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
357 ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
358+ ElementMatrix* CreateKMatrixAdjointHoriz(void);
359 ElementVector* CreateDVectorDiagnosticHoriz(void);
360 ElementVector* CreateDVectorDiagnosticStokes(void);
361 ElementMatrix* CreateKMatrixDiagnosticHutter(void);
362@@ -274,6 +275,9 @@
363
364 #ifdef _HAVE_CONTROL_
365 ElementVector* CreatePVectorAdjointHoriz(void);
366+ ElementMatrix* CreateKMatrixAdjointMacAyeal2d(void);
367+ ElementMatrix* CreateKMatrixAdjointPattyn(void);
368+ ElementMatrix* CreateKMatrixAdjointStokes(void);
369 ElementVector* CreatePVectorAdjointMacAyeal(void);
370 ElementVector* CreatePVectorAdjointPattyn(void);
371 ElementVector* CreatePVectorAdjointStokes(void);
Note: See TracBrowser for help on using the repository browser.