source: issm/oecreview/Archive/17984-18295/ISSM-18059-18060.diff@ 18296

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

Added 17984-18295
File size: 39.9 KB
RevLine 
[18296]1Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18059)
4+++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 18060)
5@@ -922,7 +922,9 @@
6 element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
7
8 /*Check that control_type is supported*/
9- if(control_type!=MaterialsRheologyBbarEnum || control_type!=FrictionCoefficientEnum || control_type!=DamageDbarEnum){
10+ if(control_type!=MaterialsRheologyBbarEnum &&
11+ control_type!=FrictionCoefficientEnum &&
12+ control_type!=DamageDbarEnum){
13 _error_("Control "<<EnumToStringx(control_type)<<" not supported");
14 }
15
16@@ -974,7 +976,7 @@
17 }/*}}}*/
18 void AdjointHorizAnalysis::GradientJDragGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
19
20- /*return if floating*/
21+ /*return if floating (gradient is 0)*/
22 if(element->IsFloating()) return;
23
24 /*Intermediaries*/
25@@ -1011,7 +1013,7 @@
26
27 /*Initialize some vectors*/
28 IssmDouble* dbasis = xNew<IssmDouble>(2*numvertices);
29- IssmDouble* ge = xNew<IssmDouble>(numvertices);
30+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
31 int* vertexpidlist = xNew<int>(numvertices);
32
33 /*Retrieve all inputs we will be needing: */
34@@ -1055,16 +1057,309 @@
35
36 }/*}}}*/
37 void AdjointHorizAnalysis::GradientJBGradient(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
38- _error_("not implemented yet");
39+
40+ /*Intermediaries*/
41+ int domaintype,dim;
42+ Element* basalelement;
43+
44+ /*Get basal element*/
45+ element->FindParam(&domaintype,DomainTypeEnum);
46+ switch(domaintype){
47+ case Domain2DhorizontalEnum:
48+ basalelement = element;
49+ dim = 2;
50+ break;
51+ case Domain2DverticalEnum:
52+ if(!element->IsOnBase()) return;
53+ basalelement = element->SpawnBasalElement();
54+ dim = 1;
55+ break;
56+ case Domain3DEnum:
57+ if(!element->IsOnBase()) return;
58+ basalelement = element->SpawnBasalElement();
59+ dim = 2;
60+ break;
61+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
62+ }
63+
64+ /*Intermediaries*/
65+ IssmDouble Jdet,weight;
66+ IssmDouble dk[3];
67+ IssmDouble *xyz_list= NULL;
68+
69+ /*Fetch number of vertices for this finite element*/
70+ int numvertices = basalelement->GetNumberOfVertices();
71+
72+ /*Initialize some vectors*/
73+ IssmDouble* dbasis = xNew<IssmDouble>(2*numvertices);
74+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
75+ int* vertexpidlist = xNew<int>(numvertices);
76+
77+ /*Retrieve all inputs we will be needing: */
78+ basalelement->GetVerticesCoordinates(&xyz_list);
79+ basalelement->GradientIndexing(&vertexpidlist[0],control_index);
80+ Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
81+ Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
82+
83+ /* Start looping on the number of gaussian points: */
84+ Gauss* gauss=basalelement->NewGauss(2);
85+ for(int ig=gauss->begin();ig<gauss->end();ig++){
86+ gauss->GaussPoint(ig);
87+
88+ basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
89+ basalelement->NodalFunctionsP1Derivatives(dbasis,xyz_list,gauss);
90+ weights_input->GetInputValue(&weight,gauss,RheologyBbarAbsGradientEnum);
91+
92+ /*Build alpha_complement_list: */
93+ rheologyb_input->GetInputDerivativeValue(&dk[0],xyz_list,gauss);
94+
95+ /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
96+ for(int i=0;i<numvertices;i++){
97+ if(dim==2){
98+ ge[i]+=-weight*Jdet*gauss->weight*(dbasis[0*numvertices+i]*dk[0]+dbasis[1*numvertices+i]*dk[1]);
99+ }
100+ else{
101+ ge[i]+=-weight*Jdet*gauss->weight*dbasis[0*numvertices+i]*dk[0];
102+ }
103+ _assert_(!xIsNan<IssmDouble>(ge[i]));
104+ }
105+ }
106+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
107+
108+ /*Clean up and return*/
109+ xDelete<IssmDouble>(xyz_list);
110+ xDelete<IssmDouble>(dbasis);
111+ xDelete<IssmDouble>(ge);
112+ xDelete<int>(vertexpidlist);
113+ delete gauss;
114+ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
115 }/*}}}*/
116 void AdjointHorizAnalysis::GradientJDragSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
117- _error_("not implemented yet");
118+
119+ /*return if floating (gradient is 0)*/
120+ if(element->IsFloating()) return;
121+
122+ /*Intermediaries*/
123+ int domaintype,dim;
124+ Element* basalelement;
125+
126+ /*Get basal element*/
127+ element->FindParam(&domaintype,DomainTypeEnum);
128+ switch(domaintype){
129+ case Domain2DhorizontalEnum:
130+ basalelement = element;
131+ dim = 2;
132+ break;
133+ case Domain2DverticalEnum:
134+ if(!element->IsOnBase()) return;
135+ basalelement = element->SpawnBasalElement();
136+ dim = 1;
137+ break;
138+ case Domain3DEnum:
139+ if(!element->IsOnBase()) return;
140+ basalelement = element->SpawnBasalElement();
141+ dim = 2;
142+ break;
143+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
144+ }
145+
146+ /*Intermediaries*/
147+ IssmDouble Jdet,weight;
148+ IssmDouble drag,dalpha2dk;
149+ IssmDouble vx,vy,lambda,mu;
150+ IssmDouble *xyz_list= NULL;
151+
152+ /*Fetch number of vertices for this finite element*/
153+ int numvertices = basalelement->GetNumberOfVertices();
154+
155+ /*Initialize some vectors*/
156+ IssmDouble* basis = xNew<IssmDouble>(numvertices);
157+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
158+ int* vertexpidlist = xNew<int>(numvertices);
159+
160+ /*Build friction element, needed later: */
161+ Friction* friction=new Friction(basalelement,dim);
162+
163+ /*Retrieve all inputs we will be needing: */
164+ basalelement->GetVerticesCoordinates(&xyz_list);
165+ basalelement->GradientIndexing(&vertexpidlist[0],control_index);
166+ Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input);
167+ Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input);
168+ Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input);
169+ Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input);
170+ Input* dragcoeff_input = basalelement->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input);
171+
172+ /* Start looping on the number of gaussian points: */
173+ Gauss* gauss=basalelement->NewGauss(4);
174+ for(int ig=gauss->begin();ig<gauss->end();ig++){
175+ gauss->GaussPoint(ig);
176+
177+ adjointx_input->GetInputValue(&lambda, gauss);
178+ adjointy_input->GetInputValue(&mu, gauss);
179+ vx_input->GetInputValue(&vx,gauss);
180+ vy_input->GetInputValue(&vy,gauss);
181+ dragcoeff_input->GetInputValue(&drag, gauss);
182+
183+ friction->GetAlphaComplement(&dalpha2dk,gauss);
184+
185+ basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
186+ basalelement->NodalFunctionsP1(basis,gauss);
187+
188+ /*Build gradient vector (actually -dJ/dD): */
189+ for(int i=0;i<numvertices;i++){
190+ ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
191+ _assert_(!xIsNan<IssmDouble>(ge[i]));
192+ }
193+ }
194+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
195+
196+ /*Clean up and return*/
197+ xDelete<IssmDouble>(xyz_list);
198+ xDelete<IssmDouble>(basis);
199+ xDelete<IssmDouble>(ge);
200+ xDelete<int>(vertexpidlist);
201+ delete gauss;
202+ delete friction;
203+ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
204 }/*}}}*/
205 void AdjointHorizAnalysis::GradientJDragHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
206- _error_("not implemented yet");
207+
208+ /*return if floating or not on bed (gradient is 0)*/
209+ if(element->IsFloating()) return;
210+ if(!element->IsOnBase()) return;
211+
212+ /*Intermediaries*/
213+ int dim=3;
214+ IssmDouble Jdet,weight;
215+ IssmDouble drag,dalpha2dk;
216+ IssmDouble vx,vy,lambda,mu;
217+ IssmDouble *xyz_list_base= NULL;
218+
219+ /*Fetch number of vertices for this finite element*/
220+ int numvertices = element->GetNumberOfVertices();
221+
222+ /*Initialize some vectors*/
223+ IssmDouble* basis = xNew<IssmDouble>(numvertices);
224+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
225+ int* vertexpidlist = xNew<int>(numvertices);
226+
227+ /*Build friction element, needed later: */
228+ Friction* friction=new Friction(element,dim);
229+
230+ /*Retrieve all inputs we will be needing: */
231+ element->GetVerticesCoordinatesBase(&xyz_list_base);
232+ element->GradientIndexing(&vertexpidlist[0],control_index);
233+ Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
234+ Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
235+ Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input);
236+ Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input);
237+ Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input);
238+
239+ /* Start looping on the number of gaussian points: */
240+ Gauss* gauss=element->NewGaussBase(4);
241+ for(int ig=gauss->begin();ig<gauss->end();ig++){
242+ gauss->GaussPoint(ig);
243+
244+ adjointx_input->GetInputValue(&lambda, gauss);
245+ adjointy_input->GetInputValue(&mu, gauss);
246+ vx_input->GetInputValue(&vx,gauss);
247+ vy_input->GetInputValue(&vy,gauss);
248+ dragcoeff_input->GetInputValue(&drag, gauss);
249+
250+ friction->GetAlphaComplement(&dalpha2dk,gauss);
251+
252+ element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
253+ element->NodalFunctionsP1(basis,gauss);
254+
255+ /*Build gradient vector (actually -dJ/dD): */
256+ for(int i=0;i<numvertices;i++){
257+ ge[i]+=-2.*drag*dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
258+ _assert_(!xIsNan<IssmDouble>(ge[i]));
259+ }
260+ }
261+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
262+
263+ /*Clean up and return*/
264+ xDelete<IssmDouble>(xyz_list_base);
265+ xDelete<IssmDouble>(basis);
266+ xDelete<IssmDouble>(ge);
267+ xDelete<int>(vertexpidlist);
268+ delete gauss;
269+ delete friction;
270 }/*}}}*/
271 void AdjointHorizAnalysis::GradientJDragFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
272- _error_("not implemented yet");
273+
274+ /*return if floating or not on bed (gradient is 0)*/
275+ if(element->IsFloating()) return;
276+ if(!element->IsOnBase()) return;
277+
278+ /*Intermediaries*/
279+ int dim=3;
280+ IssmDouble Jdet,weight;
281+ IssmDouble drag,dalpha2dk,normal[3];
282+ IssmDouble vx,vy,vz,lambda,mu,xi;
283+ IssmDouble *xyz_list_base= NULL;
284+
285+ /*Fetch number of vertices for this finite element*/
286+ int numvertices = element->GetNumberOfVertices();
287+
288+ /*Initialize some vectors*/
289+ IssmDouble* basis = xNew<IssmDouble>(numvertices);
290+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
291+ int* vertexpidlist = xNew<int>(numvertices);
292+
293+ /*Build friction element, needed later: */
294+ Friction* friction=new Friction(element,dim);
295+
296+ /*Retrieve all inputs we will be needing: */
297+ element->GetVerticesCoordinatesBase(&xyz_list_base);
298+ element->GradientIndexing(&vertexpidlist[0],control_index);
299+ Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
300+ Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
301+ Input* vz_input = element->GetInput(VzEnum); _assert_(vy_input);
302+ Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input);
303+ Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input);
304+ Input* adjointz_input = element->GetInput(AdjointzEnum); _assert_(adjointz_input);
305+ Input* dragcoeff_input = element->GetInput(FrictionCoefficientEnum); _assert_(dragcoeff_input);
306+
307+ /* Start looping on the number of gaussian points: */
308+ Gauss* gauss=element->NewGaussBase(4);
309+ for(int ig=gauss->begin();ig<gauss->end();ig++){
310+ gauss->GaussPoint(ig);
311+
312+ adjointx_input->GetInputValue(&lambda, gauss);
313+ adjointy_input->GetInputValue(&mu, gauss);
314+ adjointz_input->GetInputValue(&xi ,gauss);
315+ vx_input->GetInputValue(&vx,gauss);
316+ vy_input->GetInputValue(&vy,gauss);
317+ vz_input->GetInputValue(&vz,gauss);
318+ dragcoeff_input->GetInputValue(&drag, gauss);
319+
320+ friction->GetAlphaComplement(&dalpha2dk,gauss);
321+ element->NormalBase(&normal[0],xyz_list_base);
322+
323+ element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
324+ element->NodalFunctionsP1(basis,gauss);
325+
326+ /*Build gradient vector (actually -dJ/dk): */
327+ for(int i=0;i<numvertices;i++){
328+ ge[i]+=(
329+ -lambda*(2*drag*dalpha2dk*(vx - vz*normal[0]*normal[2]))
330+ -mu *(2*drag*dalpha2dk*(vy - vz*normal[1]*normal[2]))
331+ -xi *(2*drag*dalpha2dk*(-vx*normal[0]*normal[2]-vy*normal[1]*normal[2]))
332+ )*Jdet*gauss->weight*basis[i];
333+ _assert_(!xIsNan<IssmDouble>(ge[i]));
334+ }
335+ }
336+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
337+
338+ /*Clean up and return*/
339+ xDelete<IssmDouble>(xyz_list_base);
340+ xDelete<IssmDouble>(basis);
341+ xDelete<IssmDouble>(ge);
342+ xDelete<int>(vertexpidlist);
343+ delete gauss;
344+ delete friction;
345 }/*}}}*/
346 void AdjointHorizAnalysis::GradientJBbarSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
347
348@@ -1103,25 +1398,24 @@
349
350 /*Initialize some vectors*/
351 IssmDouble* basis = xNew<IssmDouble>(numvertices);
352- IssmDouble* ge = xNew<IssmDouble>(numvertices);
353+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
354 int* vertexpidlist = xNew<int>(numvertices);
355
356 /*Retrieve all inputs we will be needing: */
357 basalelement->GetVerticesCoordinates(&xyz_list);
358 basalelement->GradientIndexing(&vertexpidlist[0],control_index);
359- Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
360- Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
361- Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
362- Input* adjointx_input = element->GetInput(AdjointxEnum); _assert_(adjointx_input);
363- Input* adjointy_input = element->GetInput(AdjointyEnum); _assert_(adjointy_input);
364- Input* rheologyb_input = element->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
365+ Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input);
366+ Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input);
367+ Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input);
368+ Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input);
369+ Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input);
370+ Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
371
372 /* Start looping on the number of gaussian points: */
373 Gauss* gauss=basalelement->NewGauss(4);
374 for(int ig=gauss->begin();ig<gauss->end();ig++){
375 gauss->GaussPoint(ig);
376
377-
378 thickness_input->GetInputValue(&thickness,gauss);
379 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
380 vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
381@@ -1152,13 +1446,96 @@
382 if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
383 }/*}}}*/
384 void AdjointHorizAnalysis::GradientJBbarHO(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
385- _error_("not implemented yet");
386+
387+ /*WARNING: We use SSA as an estimate for now*/
388+ this->GradientJBbarSSA(element,gradient,control_index);
389 }/*}}}*/
390 void AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
391- _error_("not implemented yet");
392+ /*WARNING: We use SSA as an estimate for now*/
393+ this->GradientJBbarSSA(element,gradient,control_index);
394 }/*}}}*/
395 void AdjointHorizAnalysis::GradientJDSSA(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
396- _error_("not implemented yet");
397+
398+ /*Intermediaries*/
399+ int domaintype,dim;
400+ Element* basalelement;
401+
402+ /*Get basal element*/
403+ element->FindParam(&domaintype,DomainTypeEnum);
404+ switch(domaintype){
405+ case Domain2DhorizontalEnum:
406+ basalelement = element;
407+ dim = 2;
408+ break;
409+ case Domain2DverticalEnum:
410+ if(!element->IsOnBase()) return;
411+ basalelement = element->SpawnBasalElement();
412+ dim = 1;
413+ break;
414+ case Domain3DEnum:
415+ if(!element->IsOnBase()) return;
416+ basalelement = element->SpawnBasalElement();
417+ dim = 2;
418+ break;
419+ default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
420+ }
421+
422+ /*Intermediaries*/
423+ IssmDouble Jdet,weight;
424+ IssmDouble thickness,dmudD;
425+ IssmDouble dvx[3],dvy[3],dadjx[3],dadjy[3];
426+ IssmDouble *xyz_list= NULL;
427+
428+ /*Fetch number of vertices for this finite element*/
429+ int numvertices = basalelement->GetNumberOfVertices();
430+
431+ /*Initialize some vectors*/
432+ IssmDouble* basis = xNew<IssmDouble>(numvertices);
433+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
434+ int* vertexpidlist = xNew<int>(numvertices);
435+
436+ /*Retrieve all inputs we will be needing: */
437+ basalelement->GetVerticesCoordinates(&xyz_list);
438+ basalelement->GradientIndexing(&vertexpidlist[0],control_index);
439+ Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input);
440+ Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input);
441+ Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input);
442+ Input* adjointx_input = basalelement->GetInput(AdjointxEnum); _assert_(adjointx_input);
443+ Input* adjointy_input = basalelement->GetInput(AdjointyEnum); _assert_(adjointy_input);
444+ Input* rheologyb_input = basalelement->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
445+
446+ /* Start looping on the number of gaussian points: */
447+ Gauss* gauss=basalelement->NewGauss(4);
448+ for(int ig=gauss->begin();ig<gauss->end();ig++){
449+ gauss->GaussPoint(ig);
450+
451+ thickness_input->GetInputValue(&thickness,gauss);
452+ vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
453+ vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
454+ adjointx_input->GetInputDerivativeValue(&dadjx[0],xyz_list,gauss);
455+ adjointy_input->GetInputDerivativeValue(&dadjy[0],xyz_list,gauss);
456+
457+ basalelement->dViscositydDSSA(&dmudD,dim,xyz_list,gauss,vx_input,vy_input);
458+
459+ basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
460+ basalelement->NodalFunctionsP1(basis,gauss);
461+
462+ for(int i=0;i<numvertices;i++){
463+ ge[i]+=-dmudD*thickness*(
464+ (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
465+ )*Jdet*gauss->weight*basis[i];
466+ _assert_(!xIsNan<IssmDouble>(ge[i]));
467+ }
468+ }
469+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
470+
471+ /*Clean up and return*/
472+ xDelete<IssmDouble>(xyz_list);
473+ xDelete<IssmDouble>(basis);
474+ xDelete<IssmDouble>(ge);
475+ xDelete<int>(vertexpidlist);
476+ delete gauss;
477+ if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
478 }/*}}}*/
479 void AdjointHorizAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
480 int approximation;
481Index: ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp
482===================================================================
483--- ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 18059)
484+++ ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 18060)
485@@ -150,8 +150,171 @@
486 _error_("not implemented yet");
487 }/*}}}*/
488 void AdjointBalancethicknessAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
489- _error_("Not implemented yet");
490+ /*The gradient of the cost function is calculated in 2 parts.
491+ *
492+ * dJ \partial J \partial lambda^T(KU-F)
493+ * -- = ---------- + ------------------------
494+ * dk \partial k \parial k
495+ *
496+ * */
497+
498+ /*If on water, grad = 0: */
499+ if(!element->IsIceInElement()) return;
500+
501+ /*Get list of cost functions*/
502+ int *responses = NULL;
503+ int num_responses,resp;
504+ element->FindParam(&num_responses,InversionNumCostFunctionsEnum);
505+ element->FindParam(&responses,NULL,InversionCostFunctionsEnum);
506+
507+ /*Check that control_type is supported*/
508+ if(control_type!=VxEnum &&
509+ control_type!=VyEnum &&
510+ control_type!=BalancethicknessThickeningRateEnum){
511+ _error_("Control "<<EnumToStringx(control_type)<<" not supported");
512+ }
513+
514+ /*Deal with first part (partial derivative a J with respect to k)*/
515+ for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
516+ case ThicknessAbsMisfitEnum: /*Nothing, \partial J/\partial k = 0*/ break;
517+ case ThicknessAbsGradientEnum: /*Nothing, \partial J/\partial k = 0*/ break;
518+ case ThicknessAlongGradientEnum: /*Nothing, \partial J/\partial k = 0*/ break;
519+ case ThicknessAcrossGradientEnum: /*Nothing, \partial J/\partial k = 0*/ break;
520+ default: _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
521+ }
522+
523+ /*Deal with second term*/
524+ switch(control_type){
525+ case BalancethicknessThickeningRateEnum: GradientJDhDt(element,gradient,control_index); break;
526+ case VxEnum: GradientJVx( element,gradient,control_index); break;
527+ case VyEnum: GradientJVy( element,gradient,control_index); break;
528+ default: _error_("control type not supported yet: " << EnumToStringx(control_type));
529+ }
530+
531+ /*Clean up and return*/
532+ xDelete<int>(responses);
533+
534 }/*}}}*/
535+void AdjointBalancethicknessAnalysis::GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
536+
537+ /*Intermediaries*/
538+ IssmDouble Jdet,weight;
539+ IssmDouble thickness,Dlambda[3],dp[3];
540+ IssmDouble *xyz_list= NULL;
541+
542+ /*Fetch number of vertices for this finite element*/
543+ int numvertices = element->GetNumberOfVertices();
544+
545+ /*Initialize some vectors*/
546+ IssmDouble* basis = xNew<IssmDouble>(numvertices);
547+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
548+ int* vertexpidlist = xNew<int>(numvertices);
549+
550+ /*Retrieve all inputs we will be needing: */
551+ element->GetVerticesCoordinates(&xyz_list);
552+ element->GradientIndexing(&vertexpidlist[0],control_index);
553+ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
554+ Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input);
555+
556+ /* Start looping on the number of gaussian points: */
557+ Gauss* gauss=element->NewGauss(4);
558+ for(int ig=gauss->begin();ig<gauss->end();ig++){
559+ gauss->GaussPoint(ig);
560+
561+ adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
562+ thickness_input->GetInputValue(&thickness, gauss);
563+ thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
564+
565+ element->JacobianDeterminant(&Jdet,xyz_list,gauss);
566+ element->NodalFunctionsP1(basis,gauss);
567+
568+ /*Build gradient vector (actually -dJ/dD): */
569+ for(int i=0;i<numvertices;i++){
570+ ge[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
571+ _assert_(!xIsNan<IssmDouble>(ge[i]));
572+ }
573+ }
574+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
575+
576+ /*Clean up and return*/
577+ xDelete<IssmDouble>(xyz_list);
578+ xDelete<IssmDouble>(basis);
579+ xDelete<IssmDouble>(ge);
580+ xDelete<int>(vertexpidlist);
581+ delete gauss;
582+}/*}}}*/
583+void AdjointBalancethicknessAnalysis::GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
584+
585+ /*Intermediaries*/
586+ IssmDouble Jdet,weight;
587+ IssmDouble thickness,Dlambda[3],dp[3];
588+ IssmDouble *xyz_list= NULL;
589+
590+ /*Fetch number of vertices for this finite element*/
591+ int numvertices = element->GetNumberOfVertices();
592+
593+ /*Initialize some vectors*/
594+ IssmDouble* basis = xNew<IssmDouble>(numvertices);
595+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
596+ int* vertexpidlist = xNew<int>(numvertices);
597+
598+ /*Retrieve all inputs we will be needing: */
599+ element->GetVerticesCoordinates(&xyz_list);
600+ element->GradientIndexing(&vertexpidlist[0],control_index);
601+ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input);
602+ Input* adjoint_input = element->GetInput(AdjointEnum); _assert_(adjoint_input);
603+
604+ /* Start looping on the number of gaussian points: */
605+ Gauss* gauss=element->NewGauss(4);
606+ for(int ig=gauss->begin();ig<gauss->end();ig++){
607+ gauss->GaussPoint(ig);
608+
609+ adjoint_input->GetInputDerivativeValue(&Dlambda[0],xyz_list,gauss);
610+ thickness_input->GetInputValue(&thickness, gauss);
611+ thickness_input->GetInputDerivativeValue(&dp[0],xyz_list,gauss);
612+
613+ element->JacobianDeterminant(&Jdet,xyz_list,gauss);
614+ element->NodalFunctionsP1(basis,gauss);
615+
616+ /*Build gradient vector (actually -dJ/dvy): */
617+ for(int i=0;i<numvertices;i++){
618+ ge[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
619+ _assert_(!xIsNan<IssmDouble>(ge[i]));
620+ }
621+ }
622+ gradient->SetValues(numvertices,vertexpidlist,ge,ADD_VAL);
623+
624+ /*Clean up and return*/
625+ xDelete<IssmDouble>(xyz_list);
626+ xDelete<IssmDouble>(basis);
627+ xDelete<IssmDouble>(ge);
628+ xDelete<int>(vertexpidlist);
629+ delete gauss;
630+}/*}}}*/
631+void AdjointBalancethicknessAnalysis::GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
632+
633+ /*Fetch number of vertices for this finite element*/
634+ int numvertices = element->GetNumberOfVertices();
635+
636+ /*Initialize some vectors*/
637+ IssmDouble* ge = xNewZeroInit<IssmDouble>(numvertices);
638+ IssmDouble* lambda = xNew<IssmDouble>(numvertices);
639+ int* vertexpidlist = xNew<int>(numvertices);
640+
641+ /*Retrieve all inputs we will be needing: */
642+ element->GradientIndexing(&vertexpidlist[0],control_index);
643+ element->GetInputListOnVertices(lambda,AdjointEnum);
644+ for(int i=0;i<numvertices;i++){
645+ ge[i]=-lambda[i];
646+ _assert_(!xIsNan<IssmDouble>(ge[i]));
647+ }
648+ gradient->SetValues(numvertices,vertexpidlist,ge,INS_VAL);
649+
650+ /*Clean up and return*/
651+ xDelete<IssmDouble>(ge);
652+ xDelete<IssmDouble>(lambda);
653+ xDelete<int>(vertexpidlist);
654+}/*}}}*/
655 void AdjointBalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
656
657 int domaintype;
658Index: ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h
659===================================================================
660--- ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 18059)
661+++ ../trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.h (revision 18060)
662@@ -27,6 +27,9 @@
663 ElementVector* CreatePVector(Element* element);
664 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
665 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
666+ void GradientJVx(Element* element,Vector<IssmDouble>* gradient,int control_index);
667+ void GradientJVy(Element* element,Vector<IssmDouble>* gradient,int control_index);
668+ void GradientJDhDt(Element* element,Vector<IssmDouble>* gradient,int control_index);
669 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
670 void UpdateConstraints(FemModel* femmodel);
671 };
672Index: ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
673===================================================================
674--- ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18059)
675+++ ../trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp (revision 18060)
676@@ -8,8 +8,8 @@
677
678 void Gradjx(Vector<IssmDouble>** pgradient,IssmDouble** pnorm_list, Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters){
679
680- int i,j,numberofvertices;
681- int num_controls;
682+ int numberofvertices;
683+ int num_controls,analysisenum;
684 IssmDouble norm_inf;
685 IssmDouble *norm_list = NULL;
686 int *control_type = NULL;
687@@ -21,38 +21,39 @@
688 parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
689 numberofvertices=vertices->NumberOfVertices();
690
691+ /*Get current analysis*/
692+ parameters->FindParam(&analysisenum,AnalysisTypeEnum);
693+ Analysis* analysis = EnumToAnalysis(analysisenum);
694+
695 /*Allocate gradient_list */
696 gradient_list = xNew<Vector<IssmDouble>*>(num_controls);
697 norm_list = xNew<IssmDouble>(num_controls);
698- for(i=0;i<num_controls;i++){
699+ for(int i=0;i<num_controls;i++){
700 gradient_list[i]=new Vector<IssmDouble>(num_controls*numberofvertices);
701 }
702 gradient=new Vector<IssmDouble>(num_controls*numberofvertices);
703
704 /*Compute all gradient_list*/
705- for(i=0;i<num_controls;i++){
706-
707- for(j=0;j<elements->Size();j++){
708-
709+ for(int i=0;i<num_controls;i++){
710+ for(int j=0;j<elements->Size();j++){
711 Element* element=(Element*)elements->GetObjectByOffset(j);
712- element->Gradj(gradient_list[i],control_type[i],i);
713+ //element->Gradj(gradient_list[i],control_type[i],i);
714+ analysis->GradientJ(gradient_list[i],element,control_type[i],i);
715 }
716-
717 gradient_list[i]->Assemble();
718-
719 norm_list[i]=gradient_list[i]->Norm(NORM_INF);
720 }
721
722 /*Add all gradient_list together*/
723- for(i=0;i<num_controls;i++){
724+ for(int i=0;i<num_controls;i++){
725 gradient->AXPY(gradient_list[i],1.0);
726 delete gradient_list[i];
727 }
728
729 /*Check that gradient is clean*/
730 norm_inf=gradient->Norm(NORM_INF);
731- if(norm_inf<=0) _error_("||∂J/∂α||∞ = 0 gradient norm is zero");
732- if(xIsNan<IssmDouble>(norm_inf))_error_("||∂J/∂α||∞ = NaN gradient norm is NaN");
733+ if(norm_inf<=0) _error_("||dJ/dk|| = 0 gradient norm is zero");
734+ if(xIsNan<IssmDouble>(norm_inf))_error_("||dJ/dk|| = NaN gradient norm is NaN");
735
736 /*Clean-up and assign output pointer*/
737 if(pnorm_list){
738Index: ../trunk-jpl/src/c/modules/Gradjx/Gradjx.h
739===================================================================
740--- ../trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 18059)
741+++ ../trunk-jpl/src/c/modules/Gradjx/Gradjx.h (revision 18060)
742@@ -6,6 +6,7 @@
743 #define _GRADJX_H
744
745 #include "../../classes/classes.h"
746+#include "../../analyses/analyses.h"
747
748 /* local prototypes: */
749 void Gradjx(Vector<IssmDouble>** pgrad_g,IssmDouble** pgrad_norm,Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials,Parameters* parameters);
750Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp
751===================================================================
752--- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18059)
753+++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 18060)
754@@ -284,21 +284,10 @@
755 }
756 /*}}}*/
757 /*FUNCTION Matice::GetViscosity_B {{{*/
758-void Matice::GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){
759- /*From a string tensor and a material object, return viscosity, using Glen's flow law.
760- (1-D)
761- viscosity= -----------------------
762- 2 eps_eff ^[(n-1)/n]
763+void Matice::GetViscosity_B(IssmDouble* pdmudB,IssmDouble eps_eff){
764
765- where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
766- and n the flow law exponent.
767-
768- If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
769- return 10^14, initial viscosity.
770- */
771-
772 /*output: */
773- IssmDouble viscosity;
774+ IssmDouble dmudB;
775
776 /*Intermediary: */
777 IssmDouble D=0.,n;
778@@ -310,27 +299,44 @@
779 _assert_(D>=0. && D<1.);
780 }
781
782- if (n==1.){
783- /*Linear Viscous behavior (Newtonian fluid) viscosity=B/2: */
784- viscosity=(1.-D)/2.;
785+ if(n==1.){
786+ /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
787+ dmudB=(1.-D)/2.;
788 }
789 else{
790+ if(eps_eff==0.) dmudB = 0.;
791+ else dmudB = (1.-D)/(2.*pow(eps_eff,(n-1.)/n));
792+ }
793
794- /*if no strain rate, return maximum viscosity*/
795- if(eps_eff==0.){
796- viscosity = 1.e+14/2.;
797- }
798+ /*Return: */
799+ *pdmudB=dmudB;
800+}
801+/*}}}*/
802+/*FUNCTION Matice::GetViscosity_D {{{*/
803+void Matice::GetViscosity_D(IssmDouble* pdmudD,IssmDouble eps_eff){
804
805- else{
806- viscosity=(1.-D)/(2.*pow(eps_eff,(n-1.)/n));
807- }
808+ /*output: */
809+ IssmDouble dmudD;
810+
811+ /*Intermediary: */
812+ IssmDouble n,B;
813+
814+ /*Get B and n*/
815+ n=GetN(); _assert_(n>0.);
816+ B=GetBbar();
817+ _assert_(this->isdamaged);
818+
819+ if(n==1.){
820+ /*Linear Viscous behavior (Newtonian fluid) dmudB=B/2: */
821+ dmudD=-B/2.;
822 }
823+ else{
824+ if(eps_eff==0.) dmudD = 0.;
825+ else dmudD = -B/(2.*pow(eps_eff,(n-1.)/n));
826+ }
827
828- /*Checks in debugging mode*/
829- if(viscosity<=0) _error_("Negative viscosity");
830-
831 /*Return: */
832- *pviscosity=viscosity;
833+ *pdmudD=dmudD;
834 }
835 /*}}}*/
836 /*FUNCTION Matice::GetViscosityBar {{{*/
837Index: ../trunk-jpl/src/c/classes/Materials/Material.h
838===================================================================
839--- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 18059)
840+++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 18060)
841@@ -26,6 +26,7 @@
842 virtual void Configure(Elements* elements)=0;
843 virtual void GetViscosity(IssmDouble* pviscosity,IssmDouble epseff)=0;
844 virtual void GetViscosity_B(IssmDouble* pviscosity,IssmDouble epseff)=0;
845+ virtual void GetViscosity_D(IssmDouble* pviscosity,IssmDouble epseff)=0;
846 virtual void GetViscosityBar(IssmDouble* pviscosity,IssmDouble epseff)=0;
847 virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
848 virtual void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
849Index: ../trunk-jpl/src/c/classes/Materials/Matice.h
850===================================================================
851--- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 18059)
852+++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 18060)
853@@ -55,6 +55,7 @@
854 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
855 void GetViscosity(IssmDouble* pviscosity, IssmDouble eps_eff);
856 void GetViscosity_B(IssmDouble* pviscosity, IssmDouble eps_eff);
857+ void GetViscosity_D(IssmDouble* pviscosity, IssmDouble eps_eff);
858 void GetViscosityBar(IssmDouble* pviscosity, IssmDouble eps_eff);
859 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
860 void GetViscosityDComplement(IssmDouble*, IssmDouble*);
861Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
862===================================================================
863--- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18059)
864+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 18060)
865@@ -76,6 +76,7 @@
866 void Configure(Elements* elements);
867 void GetViscosity(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
868 void GetViscosity_B(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
869+ void GetViscosity_D(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
870 void GetViscosityBar(IssmDouble* pviscosity,IssmDouble eps_eff){_error_("not supported");};
871 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
872 void GetViscosityDComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
873Index: ../trunk-jpl/src/c/classes/Elements/Element.h
874===================================================================
875--- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18059)
876+++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 18060)
877@@ -60,6 +60,7 @@
878 void DeepEcho();
879 void DeleteMaterials(void);
880 void dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
881+ void dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
882 IssmDouble Divergence(void);
883 void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);
884 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);
885Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
886===================================================================
887--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18059)
888+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18060)
889@@ -1456,6 +1456,13 @@
890
891 }
892 /*}}}*/
893+void Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
894+
895+ _assert_(gauss->Enum()==GaussTriaEnum);
896+ this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
897+
898+}
899+/*}}}*/
900 void Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
901
902 _assert_(gauss->Enum()==GaussTriaEnum);
903@@ -1463,6 +1470,13 @@
904
905 }
906 /*}}}*/
907+void Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
908+
909+ _assert_(gauss->Enum()==GaussTriaEnum);
910+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
911+
912+}
913+/*}}}*/
914 void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
915
916 _assert_(gauss->Enum()==GaussTriaEnum);
917Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
918===================================================================
919--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18059)
920+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18060)
921@@ -190,9 +190,9 @@
922 Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
923 Gauss* NewGaussTop(int order);
924 void NodalFunctions(IssmDouble* basis,Gauss* gauss);
925- void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
926+ void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
927 void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
928- void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
929+ void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
930 void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
931 void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
932 void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
933Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
934===================================================================
935--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18059)
936+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 18060)
937@@ -2999,6 +2999,7 @@
938 } /*}}}*/
939 void Penta::GradjDragHO(Vector<IssmDouble>* gradient,int control_index){/*{{{*/
940
941+ printf("-------------- file: Penta.cpp line: %i\n",__LINE__);
942 int i,j;
943 int analysis_type;
944 int vertexpidlist[NUMVERTICES];
945Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
946===================================================================
947--- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18059)
948+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 18060)
949@@ -239,6 +239,33 @@
950
951 }
952 /*}}}*/
953+void Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
954+
955+ /*Intermediaries*/
956+ IssmDouble dmudB;
957+ IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy]; */
958+ IssmDouble epsilon1d; /* epsilon=[exx]; */
959+ IssmDouble eps_eff;
960+
961+ if(dim==2){
962+ /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
963+ this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
964+ eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
965+ }
966+ else{
967+ /* eps_eff^2 = 1/2 exx^2*/
968+ this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
969+ eps_eff = sqrt(epsilon1d*epsilon1d/2.);
970+ }
971+
972+ /*Get viscosity*/
973+ material->GetViscosity_D(&dmudB,eps_eff);
974+
975+ /*Assign output pointer*/
976+ *pdmudB=dmudB;
977+
978+}
979+/*}}}*/
980 void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
981 matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
982 }/*}}}*/
983Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp
984===================================================================
985--- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18059)
986+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 18060)
987@@ -3,7 +3,7 @@
988 */
989
990 /*Headers:*/
991-/*{{{*//*{{{*/
992+/*{{{*/
993 #ifdef HAVE_CONFIG_H
994 #include <config.h>
995 #else
Note: See TracBrowser for help on using the repository browser.