source: issm/branches/trunk-jpl-damage/src/c/objects/Elements/Tria.cpp@ 12280

Last change on this file since 12280 was 12280, checked in by cborstad, 13 years ago

merged trunk-jpl into trunk-jpl-damage through revision 12279
File size: 184.5 KB
Line 
1/*!\file Tria.cpp
2 * \brief: implementation of the Tria object
3 */
4
5/*Headers:*/
6/*{{{1*/
7#ifdef HAVE_CONFIG_H
8 #include <config.h>
9#else
10#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11#endif
12
13#include <stdio.h>
14#include <string.h>
15#include "../objects.h"
16#include "../../shared/shared.h"
17#include "../../Container/Container.h"
18#include "../../include/include.h"
19/*}}}*/
20
21/*Element macros*/
22#define NUMVERTICES 3
23
24/*Constructors/destructor/copy*/
25/*FUNCTION Tria::Tria(){{{1*/
26Tria::Tria(){
27
28 int i;
29
30 this->nodes=NULL;
31 this->matice=NULL;
32 this->matpar=NULL;
33 for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF;
34 this->inputs=NULL;
35 this->parameters=NULL;
36 this->results=NULL;
37
38}
39/*}}}*/
40/*FUNCTION Tria::Tria(int id, int sid,int index, IoModel* iomodel,int nummodels){{{1*/
41Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels)
42 :TriaRef(nummodels)
43 ,TriaHook(nummodels,index+1,iomodel){
44
45 int i;
46 /*id: */
47 this->id=tria_id;
48 this->sid=tria_sid;
49
50 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
51 this->parameters=NULL;
52
53 /*Build horizontalneighborsids list: */
54 _assert_(iomodel->Data(MeshElementconnectivityEnum));
55 //for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;
56
57 /*intialize inputs and results: */
58 this->inputs=new Inputs();
59 this->results=new Results();
60
61 /*initialize pointers:*/
62 this->nodes=NULL;
63 this->matice=NULL;
64 this->matpar=NULL;
65
66}
67/*}}}*/
68/*FUNCTION Tria::~Tria(){{{1*/
69Tria::~Tria(){
70 delete inputs;
71 delete results;
72 this->parameters=NULL;
73}
74/*}}}*/
75/*FUNCTION Tria::copy {{{1*/
76Object* Tria::copy() {
77
78 int i;
79 Tria* tria=NULL;
80
81 tria=new Tria();
82
83 //deal with TriaRef mother class
84 tria->element_type_list=(int*)xmalloc(this->numanalyses*sizeof(int));
85 for(i=0;i<this->numanalyses;i++) tria->element_type_list[i]=this->element_type_list[i];
86
87 //deal with TriaHook mother class
88 tria->numanalyses=this->numanalyses;
89 tria->hnodes=new Hook*[tria->numanalyses];
90 for(i=0;i<tria->numanalyses;i++)tria->hnodes[i]=(Hook*)this->hnodes[i]->copy();
91 tria->hmatice=(Hook*)this->hmatice->copy();
92 tria->hmatpar=(Hook*)this->hmatpar->copy();
93
94 /*deal with Tria fields: */
95 tria->id=this->id;
96 tria->sid=this->sid;
97 if(this->inputs){
98 tria->inputs=(Inputs*)this->inputs->Copy();
99 }
100 else{
101 tria->inputs=new Inputs();
102 }
103 if(this->results){
104 tria->results=(Results*)this->results->Copy();
105 }
106 else{
107 tria->results=new Results();
108 }
109 /*point parameters: */
110 tria->parameters=this->parameters;
111
112 /*recover objects: */
113 tria->nodes=(Node**)xmalloc(3*sizeof(Node*)); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes.
114 for(i=0;i<3;i++)tria->nodes[i]=this->nodes[i];
115 tria->matice=(Matice*)tria->hmatice->delivers();
116 tria->matpar=(Matpar*)tria->hmatpar->delivers();
117
118 /*neighbors: */
119 for(i=0;i<3;i++)tria->horizontalneighborsids[i]=this->horizontalneighborsids[i];
120
121 return tria;
122}
123/*}}}*/
124
125/*Other*/
126/*FUNCTION Tria::AverageOntoPartition {{{1*/
127void Tria::AverageOntoPartition(Vector* partition_contributions,Vector* partition_areas,double* vertex_response,double* qmu_part){
128
129 bool already=false;
130 int i,j;
131 int partition[NUMVERTICES];
132 int offsetsid[NUMVERTICES];
133 int offsetdof[NUMVERTICES];
134 double area;
135 double mean;
136 double values[3];
137
138 /*First, get the area: */
139 area=this->GetArea();
140
141 /*Figure out the average for this element: */
142 this->GetSidList(&offsetsid[0]);
143 this->GetDofList1(&offsetdof[0]);
144 mean=0;
145 for(i=0;i<NUMVERTICES;i++){
146 partition[i]=(int)qmu_part[offsetsid[i]];
147 mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
148 }
149
150 /*Add contribution: */
151 for(i=0;i<NUMVERTICES;i++){
152 already=false;
153 for(j=0;j<i;j++){
154 if (partition[i]==partition[j]){
155 already=true;
156 break;
157 }
158 }
159 if(!already){
160 partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
161 partition_areas->SetValue(partition[i],area,ADD_VAL);
162 };
163 }
164}
165/*}}}*/
166/*FUNCTION Tria::CreateKMatrix {{{1*/
167void Tria::CreateKMatrix(Matrix* Kff, Matrix* Kfs,Vector* df){
168
169 /*retreive parameters: */
170 ElementMatrix* Ke=NULL;
171 int analysis_type;
172 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
173
174 /*Checks in debugging mode{{{2*/
175 _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
176 /*}}}*/
177
178 /*Skip if water element*/
179 if(IsOnWater()) return;
180
181 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
182 switch(analysis_type){
183 #ifdef _HAVE_DIAGNOSTIC_
184 case DiagnosticHorizAnalysisEnum:
185 Ke=CreateKMatrixDiagnosticMacAyeal();
186 break;
187 case AdjointHorizAnalysisEnum:
188 Ke=CreateKMatrixAdjointMacAyeal();
189 break;
190 case DiagnosticHutterAnalysisEnum:
191 Ke=CreateKMatrixDiagnosticHutter();
192 break;
193 #endif
194 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
195 Ke=CreateKMatrixSlope();
196 break;
197 case PrognosticAnalysisEnum:
198 Ke=CreateKMatrixPrognostic();
199 break;
200 #ifdef _HAVE_HYDROLOGY_
201 case HydrologyAnalysisEnum:
202 Ke=CreateKMatrixHydrology();
203 break;
204 #endif
205 #ifdef _HAVE_BALANCED_
206 case BalancethicknessAnalysisEnum:
207 Ke=CreateKMatrixBalancethickness();
208 break;
209 #endif
210 #ifdef _HAVE_CONTROL_
211 case AdjointBalancethicknessAnalysisEnum:
212 Ke=CreateKMatrixAdjointBalancethickness();
213 break;
214 #endif
215 default:
216 _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
217 }
218
219 /*Add to global matrix*/
220 if(Ke){
221 Ke->AddToGlobal(Kff,Kfs);
222 delete Ke;
223 }
224}
225/*}}}*/
226/*FUNCTION Tria::CreateKMatrixMelting {{{1*/
227ElementMatrix* Tria::CreateKMatrixMelting(void){
228
229 /*Constants*/
230 const int numdof=NUMVERTICES*NDOF1;
231
232 /*Intermediaries */
233 int i,j,ig;
234 double heatcapacity,latentheat;
235 double Jdet,D_scalar;
236 double xyz_list[NUMVERTICES][3];
237 double L[3];
238 GaussTria *gauss=NULL;
239
240 /*Initialize Element matrix*/
241 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
242
243 /*Retrieve all inputs and parameters*/
244 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
245 latentheat=matpar->GetLatentHeat();
246 heatcapacity=matpar->GetHeatCapacity();
247
248 /* Start looping on the number of gauss (nodes on the bedrock) */
249 gauss=new GaussTria(2);
250 for (ig=gauss->begin();ig<gauss->end();ig++){
251
252 gauss->GaussPoint(ig);
253
254 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
255 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);
256
257 D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
258
259 TripleMultiply(&L[0],numdof,1,0,
260 &D_scalar,1,1,0,
261 &L[0],1,numdof,0,
262 &Ke->values[0],1);
263 }
264
265 /*Clean up and return*/
266 delete gauss;
267 return Ke;
268}
269/*}}}*/
270/*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/
271ElementMatrix* Tria::CreateKMatrixPrognostic(void){
272
273 switch(GetElementType()){
274 case P1Enum:
275 return CreateKMatrixPrognostic_CG();
276 case P1DGEnum:
277 return CreateKMatrixPrognostic_DG();
278 default:
279 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
280 }
281
282}
283/*}}}*/
284/*FUNCTION Tria::CreateKMatrixPrognostic_CG {{{1*/
285ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){
286
287 /*Constants*/
288 const int numdof=NDOF1*NUMVERTICES;
289
290 /*Intermediaries */
291 int stabilization;
292 int i,j,ig,dim;
293 double Jdettria,DL_scalar,dt,h;
294 double vel,vx,vy,dvxdx,dvydy;
295 double dvx[2],dvy[2];
296 double v_gauss[2]={0.0};
297 double xyz_list[NUMVERTICES][3];
298 double L[NUMVERTICES];
299 double B[2][NUMVERTICES];
300 double Bprime[2][NUMVERTICES];
301 double K[2][2] ={0.0};
302 double KDL[2][2] ={0.0};
303 double DL[2][2] ={0.0};
304 double DLprime[2][2] ={0.0};
305 GaussTria *gauss=NULL;
306
307 /*Initialize Element matrix*/
308 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
309
310 /*Retrieve all inputs and parameters*/
311 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
312 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
313 this->parameters->FindParam(&dim,MeshDimensionEnum);
314 this->parameters->FindParam(&stabilization,PrognosticStabilizationEnum);
315 Input* vxaverage_input=NULL;
316 Input* vyaverage_input=NULL;
317 if(dim==2){
318 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
319 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
320 }
321 else{
322 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
323 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
324 }
325 h=sqrt(2*this->GetArea());
326
327 /* Start looping on the number of gaussian points: */
328 gauss=new GaussTria(2);
329 for (ig=gauss->begin();ig<gauss->end();ig++){
330
331 gauss->GaussPoint(ig);
332
333 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
334 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
335
336 vxaverage_input->GetInputValue(&vx,gauss);
337 vyaverage_input->GetInputValue(&vy,gauss);
338 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
339 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
340
341 DL_scalar=gauss->weight*Jdettria;
342
343 TripleMultiply( &L[0],1,numdof,1,
344 &DL_scalar,1,1,0,
345 &L[0],1,numdof,0,
346 &Ke->values[0],1);
347
348 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
349 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
350
351 dvxdx=dvx[0];
352 dvydy=dvy[1];
353 DL_scalar=dt*gauss->weight*Jdettria;
354
355 DL[0][0]=DL_scalar*dvxdx;
356 DL[1][1]=DL_scalar*dvydy;
357 DLprime[0][0]=DL_scalar*vx;
358 DLprime[1][1]=DL_scalar*vy;
359
360 TripleMultiply( &B[0][0],2,numdof,1,
361 &DL[0][0],2,2,0,
362 &B[0][0],2,numdof,0,
363 &Ke->values[0],1);
364
365 TripleMultiply( &B[0][0],2,numdof,1,
366 &DLprime[0][0],2,2,0,
367 &Bprime[0][0],2,numdof,0,
368 &Ke->values[0],1);
369
370 if(stabilization==2){
371 /*Streamline upwinding*/
372 vel=sqrt(pow(vx,2.)+pow(vy,2.))+1.e-8;
373 K[0][0]=h/(2*vel)*vx*vx;
374 K[1][0]=h/(2*vel)*vy*vx;
375 K[0][1]=h/(2*vel)*vx*vy;
376 K[1][1]=h/(2*vel)*vy*vy;
377 }
378 else if(stabilization==1){
379 /*MacAyeal*/
380 vxaverage_input->GetInputAverage(&vx);
381 vyaverage_input->GetInputAverage(&vy);
382 K[0][0]=h/2.0*fabs(vx);
383 K[0][1]=0.;
384 K[1][0]=0.;
385 K[1][1]=h/2.0*fabs(vy);
386 }
387 if(stabilization==1 || stabilization==2){
388 KDL[0][0]=DL_scalar*K[0][0];
389 KDL[1][0]=DL_scalar*K[1][0];
390 KDL[0][1]=DL_scalar*K[0][1];
391 KDL[1][1]=DL_scalar*K[1][1];
392 TripleMultiply( &Bprime[0][0],2,numdof,1,
393 &KDL[0][0],2,2,0,
394 &Bprime[0][0],2,numdof,0,
395 &Ke->values[0],1);
396 }
397 }
398
399 /*Clean up and return*/
400 delete gauss;
401 return Ke;
402}
403/*}}}*/
404/*FUNCTION Tria::CreateKMatrixPrognostic_DG {{{1*/
405ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){
406
407 /*Constants*/
408 const int numdof=NDOF1*NUMVERTICES;
409
410 /*Intermediaries */
411 int i,j,ig,dim;
412 double xyz_list[NUMVERTICES][3];
413 double Jdettria,dt,vx,vy;
414 double L[NUMVERTICES];
415 double B[2][NUMVERTICES];
416 double Bprime[2][NUMVERTICES];
417 double DL[2][2]={0.0};
418 double DLprime[2][2]={0.0};
419 double DL_scalar;
420 GaussTria *gauss=NULL;
421
422 /*Initialize Element matrix*/
423 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
424
425 /*Retrieve all inputs and parameters*/
426 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
427 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
428 this->parameters->FindParam(&dim,MeshDimensionEnum);
429 Input* vxaverage_input=NULL;
430 Input* vyaverage_input=NULL;
431 if(dim==2){
432 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
433 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
434 }
435 else{
436 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
437 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
438 }
439
440 /* Start looping on the number of gaussian points: */
441 gauss=new GaussTria(2);
442 for (ig=gauss->begin();ig<gauss->end();ig++){
443
444 gauss->GaussPoint(ig);
445
446 vxaverage_input->GetInputValue(&vx,gauss);
447 vyaverage_input->GetInputValue(&vy,gauss);
448
449 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
450 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
451
452 DL_scalar=gauss->weight*Jdettria;
453
454 TripleMultiply( &L[0],1,numdof,1,
455 &DL_scalar,1,1,0,
456 &L[0],1,numdof,0,
457 &Ke->values[0],1);
458
459 /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
460 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
461 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
462
463 DL_scalar=-dt*gauss->weight*Jdettria;
464
465 DLprime[0][0]=DL_scalar*vx;
466 DLprime[1][1]=DL_scalar*vy;
467
468 TripleMultiply( &B[0][0],2,numdof,1,
469 &DLprime[0][0],2,2,0,
470 &Bprime[0][0],2,numdof,0,
471 &Ke->values[0],1);
472 }
473
474 /*Clean up and return*/
475 delete gauss;
476 return Ke;
477}
478/*}}}*/
479/*FUNCTION Tria::CreateKMatrixSlope {{{1*/
480ElementMatrix* Tria::CreateKMatrixSlope(void){
481
482 /*constants: */
483 const int numdof=NDOF1*NUMVERTICES;
484
485 /* Intermediaries */
486 int i,j,ig;
487 double DL_scalar,Jdet;
488 double xyz_list[NUMVERTICES][3];
489 double L[1][3];
490 GaussTria *gauss = NULL;
491
492 /*Initialize Element matrix*/
493 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
494
495 GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
496
497 /* Start looping on the number of gaussian points: */
498 gauss=new GaussTria(2);
499 for (ig=gauss->begin();ig<gauss->end();ig++){
500
501 gauss->GaussPoint(ig);
502
503 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
504 DL_scalar=gauss->weight*Jdet;
505
506 GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF1);
507
508 TripleMultiply(&L[0][0],1,3,1,
509 &DL_scalar,1,1,0,
510 &L[0][0],1,3,0,
511 &Ke->values[0],1);
512 }
513
514 /*Clean up and return*/
515 delete gauss;
516 return Ke;
517}
518/*}}}*/
519/*FUNCTION Tria::CreatePVector {{{1*/
520void Tria::CreatePVector(Vector* pf){
521
522 /*retrive parameters: */
523 ElementVector* pe=NULL;
524 int analysis_type;
525 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
526
527 /*asserts: {{{*/
528 /*if debugging mode, check that all pointers exist*/
529 _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
530 /*}}}*/
531
532 /*Skip if water element*/
533 if(IsOnWater()) return;
534
535 /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
536 switch(analysis_type){
537 #ifdef _HAVE_DIAGNOSTIC_
538 case DiagnosticHorizAnalysisEnum:
539 pe=CreatePVectorDiagnosticMacAyeal();
540 break;
541 case DiagnosticHutterAnalysisEnum:
542 pe=CreatePVectorDiagnosticHutter();
543 break;
544 #endif
545 case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
546 pe=CreatePVectorSlope();
547 break;
548 case PrognosticAnalysisEnum:
549 pe=CreatePVectorPrognostic();
550 break;
551 #ifdef _HAVE_HYDROLOGY_
552 case HydrologyAnalysisEnum:
553 pe=CreatePVectorHydrology();
554 break;
555 #endif
556 #ifdef _HAVE_BALANCED_
557 case BalancethicknessAnalysisEnum:
558 pe=CreatePVectorBalancethickness();
559 break;
560 #endif
561 #ifdef _HAVE_CONTROL_
562 case AdjointBalancethicknessAnalysisEnum:
563 pe=CreatePVectorAdjointBalancethickness();
564 break;
565 case AdjointHorizAnalysisEnum:
566 pe=CreatePVectorAdjointHoriz();
567 break;
568 #endif
569 default:
570 _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
571 }
572
573 /*Add to global Vector*/
574 if(pe){
575 pe->AddToGlobal(pf);
576 delete pe;
577 }
578}
579/*}}}*/
580/*FUNCTION Tria::CreatePVectorPrognostic{{{1*/
581ElementVector* Tria::CreatePVectorPrognostic(void){
582
583 switch(GetElementType()){
584 case P1Enum:
585 return CreatePVectorPrognostic_CG();
586 case P1DGEnum:
587 return CreatePVectorPrognostic_DG();
588 default:
589 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
590 }
591}
592/*}}}*/
593/*FUNCTION Tria::CreatePVectorPrognostic_CG {{{1*/
594ElementVector* Tria::CreatePVectorPrognostic_CG(void){
595
596 /*Constants*/
597 const int numdof=NDOF1*NUMVERTICES;
598
599 /*Intermediaries */
600 int i,j,ig;
601 double Jdettria,dt;
602 double surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
603 double xyz_list[NUMVERTICES][3];
604 double L[NUMVERTICES];
605 GaussTria* gauss=NULL;
606
607 /*Initialize Element vector*/
608 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
609
610 /*Retrieve all inputs and parameters*/
611 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
612 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
613 Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
614 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
615 Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
616 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
617
618 /*Initialize basal_melting_correction_g to 0, do not forget!:*/
619 /* Start looping on the number of gaussian points: */
620 gauss=new GaussTria(2);
621 for(ig=gauss->begin();ig<gauss->end();ig++){
622
623 gauss->GaussPoint(ig);
624
625 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
626 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
627
628 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
629 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
630 thickness_input->GetInputValue(&thickness_g,gauss);
631 if(basal_melting_correction_input)
632 basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss);
633 else
634 basal_melting_correction_g=0.;
635
636 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*L[i];
637 }
638
639 /*Clean up and return*/
640 delete gauss;
641 return pe;
642}
643/*}}}*/
644/*FUNCTION Tria::CreatePVectorPrognostic_DG {{{1*/
645ElementVector* Tria::CreatePVectorPrognostic_DG(void){
646
647 /*Constants*/
648 const int numdof=NDOF1*NUMVERTICES;
649
650 /*Intermediaries */
651 int i,j,ig;
652 double Jdettria,dt;
653 double surface_mass_balance_g,basal_melting_g,thickness_g;
654 double xyz_list[NUMVERTICES][3];
655 double L[NUMVERTICES];
656 GaussTria* gauss=NULL;
657
658 /*Initialize Element vector*/
659 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
660
661 /*Retrieve all inputs and parameters*/
662 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
663 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
664 Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
665 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
666 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
667
668 /* Start looping on the number of gaussian points: */
669 gauss=new GaussTria(2);
670 for(ig=gauss->begin();ig<gauss->end();ig++){
671
672 gauss->GaussPoint(ig);
673
674 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
675 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
676
677 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
678 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
679 thickness_input->GetInputValue(&thickness_g,gauss);
680
681 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*L[i];
682 }
683
684 /*Clean up and return*/
685 delete gauss;
686 return pe;
687}
688/*}}}*/
689/*FUNCTION Tria::CreatePVectorSlope {{{1*/
690ElementVector* Tria::CreatePVectorSlope(void){
691
692 /*Constants*/
693 const int numdof=NDOF1*NUMVERTICES;
694
695 /*Intermediaries */
696 int i,j,ig;
697 int analysis_type;
698 double Jdet;
699 double xyz_list[NUMVERTICES][3];
700 double slope[2];
701 double basis[3];
702 GaussTria* gauss=NULL;
703
704 /*Initialize Element vector*/
705 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
706
707 /*Retrieve all inputs and parameters*/
708 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
709 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
710 Input* slope_input=NULL;
711 if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==SurfaceSlopeYAnalysisEnum)){
712 slope_input=inputs->GetInput(SurfaceEnum); _assert_(slope_input);
713 }
714 if ( (analysis_type==BedSlopeXAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
715 slope_input=inputs->GetInput(BedEnum); _assert_(slope_input);
716 }
717
718 /* Start looping on the number of gaussian points: */
719 gauss=new GaussTria(2);
720 for(ig=gauss->begin();ig<gauss->end();ig++){
721
722 gauss->GaussPoint(ig);
723
724 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
725 GetNodalFunctions(basis, gauss);
726
727 slope_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
728
729 if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
730 for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*basis[i];
731 }
732 if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
733 for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*basis[i];
734 }
735 }
736
737 /*Clean up and return*/
738 delete gauss;
739 return pe;
740}
741/*}}}*/
742/*FUNCTION Tria::CreateJacobianMatrix{{{1*/
743void Tria::CreateJacobianMatrix(Matrix* Jff){
744
745 /*retrieve parameters: */
746 ElementMatrix* Ke=NULL;
747 int analysis_type;
748 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
749
750 /*Checks in debugging {{{2*/
751 _assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
752 /*}}}*/
753
754 /*Skip if water element*/
755 if(IsOnWater()) return;
756
757 /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
758 switch(analysis_type){
759#ifdef _HAVE_DIAGNOSTIC_
760 case DiagnosticHorizAnalysisEnum:
761 Ke=CreateJacobianDiagnosticMacayeal();
762 break;
763#endif
764 default:
765 _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
766 }
767
768 /*Add to global matrix*/
769 if(Ke){
770 Ke->AddToGlobal(Jff);
771 delete Ke;
772 }
773}
774/*}}}*/
775/*FUNCTION Tria::ComputeBasalStress {{{1*/
776void Tria::ComputeBasalStress(Vector* eps){
777 _error_("Not Implemented yet");
778}
779/*}}}*/
780/*FUNCTION Tria::ComputeStrainRate {{{1*/
781void Tria::ComputeStrainRate(Vector* eps){
782 _error_("Not Implemented yet");
783}
784/*}}}*/
785/*FUNCTION Tria::ComputeStressTensor {{{1*/
786void Tria::ComputeStressTensor(){
787
788 int iv;
789 double xyz_list[NUMVERTICES][3];
790 double pressure,viscosity;
791 double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
792 double sigma_xx[NUMVERTICES];
793 double sigma_yy[NUMVERTICES];
794 double sigma_zz[NUMVERTICES]={0,0,0};
795 double sigma_xy[NUMVERTICES];
796 double sigma_xz[NUMVERTICES]={0,0,0};
797 double sigma_yz[NUMVERTICES]={0,0,0};
798 GaussTria* gauss=NULL;
799
800 /* Get node coordinates and dof list: */
801 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
802
803 /*Retrieve all inputs we will be needing: */
804 Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
805 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
806 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
807
808 /* Start looping on the number of vertices: */
809 gauss=new GaussTria();
810 for (int iv=0;iv<NUMVERTICES;iv++){
811 gauss->GaussVertex(iv);
812
813 /*Compute strain rate viscosity and pressure: */
814 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
815 matice->GetViscosity2d(&viscosity,&epsilon[0]);
816 pressure_input->GetInputValue(&pressure,gauss);
817
818 /*Compute Stress*/
819 sigma_xx[iv]=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
820 sigma_yy[iv]=2*viscosity*epsilon[1]-pressure;
821 sigma_xy[iv]=2*viscosity*epsilon[2];
822 }
823
824 /*Add Stress tensor components into inputs*/
825 this->inputs->AddInput(new TriaP1Input(StressTensorxxEnum,&sigma_xx[0]));
826 this->inputs->AddInput(new TriaP1Input(StressTensorxyEnum,&sigma_xy[0]));
827 this->inputs->AddInput(new TriaP1Input(StressTensorxzEnum,&sigma_xz[0]));
828 this->inputs->AddInput(new TriaP1Input(StressTensoryyEnum,&sigma_yy[0]));
829 this->inputs->AddInput(new TriaP1Input(StressTensoryzEnum,&sigma_yz[0]));
830 this->inputs->AddInput(new TriaP1Input(StressTensorzzEnum,&sigma_zz[0]));
831
832 /*Clean up and return*/
833 delete gauss;
834}
835/*}}}*/
836/*FUNCTION Tria::Configure {{{1*/
837void Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
838
839 /*go into parameters and get the analysis_counter: */
840 int analysis_counter;
841 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
842
843 /*Get Element type*/
844 this->element_type=this->element_type_list[analysis_counter];
845
846 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
847 * datasets, using internal ids and offsets hidden in hooks: */
848 if(this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin);
849 this->hmatice->configure(materialsin);
850 this->hmatpar->configure(materialsin);
851
852 /*Now, go pick up the objects inside the hooks: */
853 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
854 else this->nodes=NULL;
855 this->matice=(Matice*)this->hmatice->delivers();
856 this->matpar=(Matpar*)this->hmatpar->delivers();
857
858 /*point parameters to real dataset: */
859 this->parameters=parametersin;
860
861 /*get inputs configured too: */
862 this->inputs->Configure(parameters);
863
864}
865/*}}}*/
866/*FUNCTION Tria::DeepEcho{{{1*/
867void Tria::DeepEcho(void){
868
869 printf("Tria:\n");
870 printf(" id: %i\n",id);
871 if(nodes){
872 nodes[0]->DeepEcho();
873 nodes[1]->DeepEcho();
874 nodes[2]->DeepEcho();
875 }
876 else printf("nodes = NULL\n");
877
878 if (matice) matice->DeepEcho();
879 else printf("matice = NULL\n");
880
881 if (matpar) matpar->DeepEcho();
882 else printf("matpar = NULL\n");
883
884 printf(" parameters\n");
885 if (parameters) parameters->DeepEcho();
886 else printf("parameters = NULL\n");
887
888 printf(" inputs\n");
889 if (inputs) inputs->DeepEcho();
890 else printf("inputs=NULL\n");
891
892 if (results) results->DeepEcho();
893 else printf("results=NULL\n");
894
895 printf("neighboor sids: \n");
896 printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
897
898 return;
899}
900/*}}}*/
901/*FUNCTION Tria::DeleteResults {{{1*/
902void Tria::DeleteResults(void){
903
904 /*Delete and reinitialize results*/
905 delete this->results;
906 this->results=new Results();
907
908}
909/*}}}*/
910/*FUNCTION Tria::Echo{{{1*/
911void Tria::Echo(void){
912 printf("Tria:\n");
913 printf(" id: %i\n",id);
914 if(nodes){
915 nodes[0]->Echo();
916 nodes[1]->Echo();
917 nodes[2]->Echo();
918 }
919 else printf("nodes = NULL\n");
920
921 if (matice) matice->Echo();
922 else printf("matice = NULL\n");
923
924 if (matpar) matpar->Echo();
925 else printf("matpar = NULL\n");
926
927 printf(" parameters\n");
928 if (parameters) parameters->Echo();
929 else printf("parameters = NULL\n");
930
931 printf(" inputs\n");
932 if (inputs) inputs->Echo();
933 else printf("inputs=NULL\n");
934
935 if (results) results->Echo();
936 else printf("results=NULL\n");
937
938 printf("neighboor sids: \n");
939 printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
940}
941/*}}}*/
942/*FUNCTION Tria::ObjectEnum{{{1*/
943int Tria::ObjectEnum(void){
944
945 return TriaEnum;
946
947}
948/*}}}*/
949/*FUNCTION Tria::GetArea {{{1*/
950double Tria::GetArea(void){
951
952 double area=0;
953 double xyz_list[NUMVERTICES][3];
954 double x1,y1,x2,y2,x3,y3;
955
956 /*Get xyz list: */
957 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
958 x1=xyz_list[0][0]; y1=xyz_list[0][1];
959 x2=xyz_list[1][0]; y2=xyz_list[1][1];
960 x3=xyz_list[2][0]; y3=xyz_list[2][1];
961
962 _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
963 return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
964}
965/*}}}*/
966/*FUNCTION Tria::GetDofList {{{1*/
967void Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){
968
969 int i,j;
970 int count=0;
971 int numberofdofs=0;
972 int* doflist=NULL;
973
974 /*First, figure out size of doflist and create it: */
975 for(i=0;i<3;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
976 doflist=(int*)xmalloc(numberofdofs*sizeof(int));
977
978 /*Populate: */
979 count=0;
980 for(i=0;i<3;i++){
981 nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
982 count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
983 }
984
985 /*Assign output pointers:*/
986 *pdoflist=doflist;
987}
988/*}}}*/
989/*FUNCTION Tria::GetDofList1 {{{1*/
990void Tria::GetDofList1(int* doflist){
991
992 int i;
993 for(i=0;i<3;i++) doflist[i]=nodes[i]->GetDofList1();
994
995}
996/*}}}*/
997/*FUNCTION Tria::GetElementType {{{1*/
998int Tria::GetElementType(){
999
1000 /*return TriaRef field*/
1001 return this->element_type;
1002
1003}
1004/*}}}*/
1005/*FUNCTION Tria::GetHorizontalNeighboorSids {{{1*/
1006int* Tria::GetHorizontalNeighboorSids(){
1007
1008 /*return TriaRef field*/
1009 return &this->horizontalneighborsids[0];
1010
1011}
1012/*}}}*/
1013/*FUNCTION Tria::GetNodeIndex {{{1*/
1014int Tria::GetNodeIndex(Node* node){
1015
1016 _assert_(nodes);
1017 for(int i=0;i<NUMVERTICES;i++){
1018 if(node==nodes[i])
1019 return i;
1020 }
1021 _error_("Node provided not found among element nodes");
1022}
1023/*}}}*/
1024/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype) {{{1*/
1025void Tria::GetInputListOnVertices(double* pvalue,int enumtype){
1026
1027 /*Intermediaries*/
1028 double value[NUMVERTICES];
1029 GaussTria *gauss = NULL;
1030
1031 /*Recover input*/
1032 Input* input=inputs->GetInput(enumtype);
1033 if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype));
1034
1035 /*Checks in debugging mode*/
1036 _assert_(pvalue);
1037
1038 /* Start looping on the number of vertices: */
1039 gauss=new GaussTria();
1040 for (int iv=0;iv<NUMVERTICES;iv++){
1041 gauss->GaussVertex(iv);
1042 input->GetInputValue(&pvalue[iv],gauss);
1043 }
1044
1045 /*clean-up*/
1046 delete gauss;
1047}
1048/*}}}*/
1049/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
1050void Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue){
1051
1052 double value[NUMVERTICES];
1053 GaussTria *gauss = NULL;
1054 Input *input = inputs->GetInput(enumtype);
1055
1056 /*Checks in debugging mode*/
1057 _assert_(pvalue);
1058
1059 /* Start looping on the number of vertices: */
1060 if (input){
1061 gauss=new GaussTria();
1062 for (int iv=0;iv<NUMVERTICES;iv++){
1063 gauss->GaussVertex(iv);
1064 input->GetInputValue(&pvalue[iv],gauss);
1065 }
1066 }
1067 else{
1068 for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
1069 }
1070
1071 /*clean-up*/
1072 delete gauss;
1073}
1074/*}}}*/
1075/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index) TO BE REMOVED{{{1*/
1076void Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index){
1077
1078 double value[NUMVERTICES];
1079 GaussTria *gauss = NULL;
1080 Input *input = inputs->GetInput(enumtype);
1081
1082 /*Checks in debugging mode*/
1083 _assert_(pvalue);
1084
1085 /* Start looping on the number of vertices: */
1086 if (input){
1087 gauss=new GaussTria();
1088 for (int iv=0;iv<NUMVERTICES;iv++){
1089 gauss->GaussVertex(iv);
1090 input->GetInputValue(&pvalue[iv],gauss,index);
1091 }
1092 }
1093 else{
1094 for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
1095 }
1096
1097 /*clean-up*/
1098 delete gauss;
1099}
1100/*}}}*/
1101/*FUNCTION Tria::GetInputValue(double* pvalue,Node* node,int enumtype) {{{1*/
1102void Tria::GetInputValue(double* pvalue,Node* node,int enumtype){
1103
1104 Input* input=inputs->GetInput(enumtype);
1105 if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype));
1106
1107 GaussTria* gauss=new GaussTria();
1108 gauss->GaussVertex(this->GetNodeIndex(node));
1109
1110 input->GetInputValue(pvalue,gauss);
1111 delete gauss;
1112}
1113/*}}}*/
1114/*FUNCTION Tria::GetSidList {{{1*/
1115void Tria::GetSidList(int* sidlist){
1116 for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->GetSidList();
1117}
1118/*}}}*/
1119/*FUNCTION Tria::GetConnectivityList {{{1*/
1120void Tria::GetConnectivityList(int* connectivity){
1121 for(int i=0;i<NUMVERTICES;i++) connectivity[i]=nodes[i]->GetConnectivity();
1122}
1123/*}}}*/
1124/*FUNCTION Tria::GetSolutionFromInputs{{{1*/
1125void Tria::GetSolutionFromInputs(Vector* solution){
1126
1127 /*retrive parameters: */
1128 int analysis_type;
1129 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
1130
1131 /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
1132 switch(analysis_type){
1133 #ifdef _HAVE_DIAGNOSTIC_
1134 case DiagnosticHorizAnalysisEnum:
1135 GetSolutionFromInputsDiagnosticHoriz(solution);
1136 break;
1137 case DiagnosticHutterAnalysisEnum:
1138 GetSolutionFromInputsDiagnosticHutter(solution);
1139 break;
1140 #endif
1141 #ifdef _HAVE_HYDROLOGY_
1142 case HydrologyAnalysisEnum:
1143 GetSolutionFromInputsHydrology(solution);
1144 break;
1145 #endif
1146 default:
1147 _error_("analysis: %s not supported yet",EnumToStringx(analysis_type));
1148 }
1149
1150}
1151/*}}}*/
1152/*FUNCTION Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{1*/
1153void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){
1154 /*Compute the 2d Strain Rate (3 components):
1155 * epsilon=[exx eyy exy] */
1156
1157 int i;
1158 double epsilonvx[3];
1159 double epsilonvy[3];
1160
1161 /*Check that both inputs have been found*/
1162 if (!vx_input || !vy_input){
1163 _error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
1164 }
1165
1166 /*Get strain rate assuming that epsilon has been allocated*/
1167 vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
1168 vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
1169
1170 /*Sum all contributions*/
1171 for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
1172}
1173/*}}}*/
1174/*FUNCTION Tria::GetVectorFromInputs{{{1*/
1175void Tria::GetVectorFromInputs(Vector* vector,int input_enum){
1176
1177 int doflist1[NUMVERTICES];
1178
1179 /*Get out if this is not an element input*/
1180 if(!IsInput(input_enum)) return;
1181
1182 /*Prepare index list*/
1183 this->GetDofList1(&doflist1[0]);
1184
1185 /*Get input (either in element or material)*/
1186 Input* input=inputs->GetInput(input_enum);
1187 if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));
1188
1189 /*We found the enum. Use its values to fill into the vector, using the vertices ids: */
1190 input->GetVectorFromInputs(vector,&doflist1[0]);
1191}
1192/*}}}*/
1193/*FUNCTION Tria::GetVectorFromResults{{{1*/
1194void Tria::GetVectorFromResults(Vector* vector,int offset,int interp){
1195
1196 /*Get result*/
1197 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(offset);
1198 if(interp==P1Enum){
1199 int doflist1[NUMVERTICES];
1200 int connectivity[NUMVERTICES];
1201 this->GetSidList(&doflist1[0]);
1202 this->GetConnectivityList(&connectivity[0]);
1203 elementresult->GetVectorFromResults(vector,&doflist1[0],&connectivity[0],NUMVERTICES);
1204 }
1205 else if(interp==P0Enum){
1206 elementresult->GetElementVectorFromResults(vector,sid);
1207 }
1208 else{
1209 printf("Interpolation %s not supported\n",EnumToStringx(interp));
1210 }
1211}
1212/*}}}*/
1213/*FUNCTION Tria::Id {{{1*/
1214int Tria::Id(){
1215
1216 return id;
1217
1218}
1219/*}}}*/
1220/*FUNCTION Tria::Sid {{{1*/
1221int Tria::Sid(){
1222
1223 return sid;
1224
1225}
1226/*}}}*/
1227/*FUNCTION Tria::InputArtificialNoise{{{1*/
1228void Tria::InputArtificialNoise(int enum_type,double min,double max){
1229
1230 Input* input=NULL;
1231
1232 /*Make a copy of the original input: */
1233 input=(Input*)this->inputs->GetInput(enum_type);
1234 if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));
1235
1236 /*ArtificialNoise: */
1237 input->ArtificialNoise(min,max);
1238}
1239/*}}}*/
1240/*FUNCTION Tria::InputConvergence{{{1*/
1241bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){
1242
1243 bool converged=true;
1244 int i;
1245 Input** new_inputs=NULL;
1246 Input** old_inputs=NULL;
1247
1248 new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
1249 old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
1250
1251 for(i=0;i<num_enums/2;i++){
1252 new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
1253 old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
1254 if(!new_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
1255 if(!old_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
1256 }
1257
1258 /*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
1259 for(i=0;i<num_criterionenums;i++){
1260 IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
1261 if(eps[i]>criterionvalues[i]) converged=false;
1262 }
1263
1264 /*clean up and return*/
1265 xfree((void**)&new_inputs);
1266 xfree((void**)&old_inputs);
1267 return converged;
1268}
1269/*}}}*/
1270/*FUNCTION Tria::InputDepthAverageAtBase {{{1*/
1271void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
1272
1273 /*New input*/
1274 Input* oldinput=NULL;
1275 Input* newinput=NULL;
1276
1277 /*copy input of enum_type*/
1278 if (object_enum==MeshElementsEnum)
1279 oldinput=(Input*)this->inputs->GetInput(enum_type);
1280 else if (object_enum==MaterialsEnum)
1281 oldinput=(Input*)this->matice->inputs->GetInput(enum_type);
1282 else
1283 _error_("object %s not supported yet",EnumToStringx(object_enum));
1284 if(!oldinput)_error_("%s%s"," could not find old input with enum: ",EnumToStringx(enum_type));
1285 newinput=(Input*)oldinput->copy();
1286
1287 /*Assign new name (average)*/
1288 newinput->ChangeEnum(average_enum_type);
1289
1290 /*Add new input to current element*/
1291 if (object_enum==MeshElementsEnum)
1292 this->inputs->AddInput((Input*)newinput);
1293 else if (object_enum==MaterialsEnum)
1294 this->matice->inputs->AddInput((Input*)newinput);
1295 else
1296 _error_("object %s not supported yet",EnumToStringx(object_enum));
1297}
1298/*}}}*/
1299/*FUNCTION Tria::InputDuplicate{{{1*/
1300void Tria::InputDuplicate(int original_enum,int new_enum){
1301
1302 /*Call inputs method*/
1303 if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);
1304
1305}
1306/*}}}*/
1307/*FUNCTION Tria::InputScale{{{1*/
1308void Tria::InputScale(int enum_type,double scale_factor){
1309
1310 Input* input=NULL;
1311
1312 /*Make a copy of the original input: */
1313 input=(Input*)this->inputs->GetInput(enum_type);
1314 if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));
1315
1316 /*Scale: */
1317 input->Scale(scale_factor);
1318}
1319/*}}}*/
1320/*FUNCTION Tria::InputToResult{{{1*/
1321void Tria::InputToResult(int enum_type,int step,double time){
1322
1323 int i;
1324 Input *input = NULL;
1325
1326 /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
1327 if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type);
1328 else input=this->inputs->GetInput(enum_type);
1329 //if (!input) _error_("Input %s not found in tria->inputs",EnumToStringx(enum_type));
1330 if(!input)return;
1331
1332 /*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result
1333 * object out of the input, with the additional step and time information: */
1334 this->results->AddObject((Object*)input->SpawnResult(step,time));
1335
1336 #ifdef _HAVE_CONTROL_
1337 if(input->ObjectEnum()==ControlInputEnum){
1338 if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
1339 }
1340 #endif
1341}
1342/*}}}*/
1343/*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{1*/
1344void Tria::InputUpdateFromConstant(int constant, int name){
1345 /*Check that name is an element input*/
1346 if (!IsInput(name)) return;
1347
1348 /*update input*/
1349 this->inputs->AddInput(new IntInput(name,constant));
1350}
1351/*}}}*/
1352/*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{1*/
1353void Tria::InputUpdateFromConstant(double constant, int name){
1354 /*Check that name is an element input*/
1355 if (!IsInput(name)) return;
1356
1357 /*update input*/
1358 this->inputs->AddInput(new DoubleInput(name,constant));
1359}
1360/*}}}*/
1361/*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{1*/
1362void Tria::InputUpdateFromConstant(bool constant, int name){
1363 /*Check that name is an element input*/
1364 if (!IsInput(name)) return;
1365
1366 /*update input*/
1367 this->inputs->AddInput(new BoolInput(name,constant));
1368}
1369/*}}}*/
1370/*FUNCTION Tria::InputUpdateFromIoModel{{{1*/
1371void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index
1372
1373 /*Intermediaries*/
1374 int i,j;
1375 int tria_vertex_ids[3];
1376 double nodeinputs[3];
1377 double cmmininputs[3];
1378 double cmmaxinputs[3];
1379 double cmthreshinputs[3];
1380 bool control_analysis=false;
1381 int num_control_type;
1382 double yts;
1383 int num_cm_responses;
1384
1385 /*Get parameters: */
1386 iomodel->Constant(&yts,ConstantsYtsEnum);
1387 iomodel->Constant(&control_analysis,InversionIscontrolEnum);
1388 if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
1389 if(control_analysis) iomodel->Constant(&num_cm_responses,InversionNumCostFunctionsEnum);
1390
1391 /*Recover vertices ids needed to initialize inputs*/
1392 for(i=0;i<3;i++){
1393 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
1394 }
1395
1396 /*Control Inputs*/
1397 #ifdef _HAVE_CONTROL_
1398 if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
1399 for(i=0;i<num_control_type;i++){
1400 switch((int)iomodel->Data(InversionControlParametersEnum)[i]){
1401 case BalancethicknessThickeningRateEnum:
1402 if (iomodel->Data(BalancethicknessThickeningRateEnum)){
1403 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[tria_vertex_ids[j]-1]/yts;
1404 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1405 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1406 for(j=0;j<3;j++)cmthreshinputs[j]=iomodel->Data(InversionThreshParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1407 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,cmthreshinputs,i+1));
1408 }
1409 break;
1410 case VxEnum:
1411 if (iomodel->Data(VxEnum)){
1412 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VxEnum)[tria_vertex_ids[j]-1]/yts;
1413 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1414 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1415 for(j=0;j<3;j++)cmthreshinputs[j]=iomodel->Data(InversionThreshParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1416 this->inputs->AddInput(new ControlInput(VxEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,cmthreshinputs,i+1));
1417 }
1418 break;
1419 case VyEnum:
1420 if (iomodel->Data(VyEnum)){
1421 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VyEnum)[tria_vertex_ids[j]-1]/yts;
1422 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1423 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1424 for(j=0;j<3;j++)cmthreshinputs[j]=iomodel->Data(InversionThreshParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1425 this->inputs->AddInput(new ControlInput(VyEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,cmthreshinputs,i+1));
1426 }
1427 break;
1428 case FrictionCoefficientEnum:
1429 if (iomodel->Data(FrictionCoefficientEnum)){
1430 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[tria_vertex_ids[j]-1];
1431 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
1432 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
1433 for(j=0;j<3;j++)cmthreshinputs[j]=iomodel->Data(InversionThreshParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1434 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,cmthreshinputs,i+1));
1435 }
1436 break;
1437 case MaterialsRheologyBbarEnum:
1438 case MaterialsRheologyZbarEnum:
1439 /*Matice will take care of it*/ break;
1440 default:
1441 _error_("Control %s not implemented yet",EnumToStringx((int)iomodel->Data(InversionControlParametersEnum)[i]));
1442 }
1443 }
1444 }
1445 #endif
1446
1447 /*DatasetInputs*/
1448 if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)){
1449
1450 /*Create inputs and add to DataSetInput*/
1451 DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
1452 for(i=0;i<num_cm_responses;i++){
1453 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(tria_vertex_ids[j]-1)*num_cm_responses+i];
1454 datasetinput->inputs->AddObject(new TriaP1Input(InversionCostFunctionsCoefficientsEnum,nodeinputs));
1455 }
1456
1457 /*Add datasetinput to element inputs*/
1458 this->inputs->AddInput(datasetinput);
1459 }
1460}
1461/*}}}*/
1462/*FUNCTION Tria::InputUpdateFromSolution {{{1*/
1463void Tria::InputUpdateFromSolution(double* solution){
1464
1465 /*retrive parameters: */
1466 int analysis_type;
1467 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
1468
1469 /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
1470 switch(analysis_type){
1471 #ifdef _HAVE_DIAGNOSTIC_
1472 case DiagnosticHorizAnalysisEnum:
1473 InputUpdateFromSolutionDiagnosticHoriz( solution);
1474 break;
1475 case DiagnosticHutterAnalysisEnum:
1476 InputUpdateFromSolutionDiagnosticHoriz( solution);
1477 break;
1478 #endif
1479 #ifdef _HAVE_CONTROL_
1480 case AdjointHorizAnalysisEnum:
1481 InputUpdateFromSolutionAdjointHoriz( solution);
1482 break;
1483 case AdjointBalancethicknessAnalysisEnum:
1484 InputUpdateFromSolutionAdjointBalancethickness( solution);
1485 break;
1486 #endif
1487 #ifdef _HAVE_HYDROLOGY_
1488 case HydrologyAnalysisEnum:
1489 InputUpdateFromSolutionHydrology(solution);
1490 break ;
1491 #endif
1492 #ifdef _HAVE_BALANCED_
1493 case BalancethicknessAnalysisEnum:
1494 InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
1495 break;
1496 #endif
1497 case BedSlopeXAnalysisEnum:
1498 InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum);
1499 break;
1500 case BedSlopeYAnalysisEnum:
1501 InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum);
1502 break;
1503 case SurfaceSlopeXAnalysisEnum:
1504 InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
1505 break;
1506 case SurfaceSlopeYAnalysisEnum:
1507 InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
1508 break;
1509 case PrognosticAnalysisEnum:
1510 InputUpdateFromSolutionPrognostic(solution);
1511 break;
1512 default:
1513 _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
1514 }
1515}
1516/*}}}*/
1517/*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{1*/
1518void Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
1519
1520 const int numdof = NDOF1*NUMVERTICES;
1521
1522 int* doflist=NULL;
1523 double values[numdof];
1524
1525 /*Get dof list: */
1526 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
1527
1528 /*Use the dof list to index into the solution vector: */
1529 for(int i=0;i<numdof;i++){
1530 values[i]=solution[doflist[i]];
1531 if(isnan(values[i])) _error_("NaN found in solution vector");
1532 }
1533
1534 /*Add input to the element: */
1535 this->inputs->AddInput(new TriaP1Input(enum_type,values));
1536
1537 /*Free ressources:*/
1538 xfree((void**)&doflist);
1539}
1540/*}}}*/
1541/*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{1*/
1542void Tria::InputUpdateFromSolutionPrognostic(double* solution){
1543
1544 /*Intermediaries*/
1545 const int numdof = NDOF1*NUMVERTICES;
1546
1547 int i,hydroadjustment;
1548 int* doflist=NULL;
1549 double rho_ice,rho_water,minthickness;
1550 double newthickness[numdof];
1551 double newbed[numdof];
1552 double newsurface[numdof];
1553 double oldbed[NUMVERTICES];
1554 double oldsurface[NUMVERTICES];
1555 double oldthickness[NUMVERTICES];
1556
1557 /*Get dof list: */
1558 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
1559
1560 /*Use the dof list to index into the solution vector: */
1561 this->parameters->FindParam(&minthickness,PrognosticMinThicknessEnum);
1562 for(i=0;i<numdof;i++){
1563 newthickness[i]=solution[doflist[i]];
1564 if(isnan(newthickness[i])) _error_("NaN found in solution vector");
1565 /*Constrain thickness to be at least 1m*/
1566 if(newthickness[i]<minthickness) newthickness[i]=minthickness;
1567 }
1568
1569 /*Get previous bed, thickness and surface*/
1570 GetInputListOnVertices(&oldbed[0],BedEnum);
1571 GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
1572 GetInputListOnVertices(&oldthickness[0],ThicknessEnum);
1573
1574 /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
1575 this->parameters->FindParam(&hydroadjustment,PrognosticHydrostaticAdjustmentEnum);
1576 rho_ice=matpar->GetRhoIce();
1577 rho_water=matpar->GetRhoWater();
1578
1579 for(i=0;i<numdof;i++) {
1580 /*If shelf: hydrostatic equilibrium*/
1581 if (this->nodes[i]->IsGrounded()){
1582 newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
1583 newbed[i]=oldbed[i]; //same bed: do nothing
1584 }
1585 else{ //this is an ice shelf
1586
1587 if(hydroadjustment==AbsoluteEnum){
1588 newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
1589 newbed[i]=newthickness[i]*(-rho_ice/rho_water);
1590 }
1591 else if(hydroadjustment==IncrementalEnum){
1592 newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
1593 newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
1594 }
1595 else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment));
1596 }
1597 }
1598
1599 /*Add input to the element: */
1600 this->inputs->AddInput(new TriaP1Input(ThicknessEnum,newthickness));
1601 this->inputs->AddInput(new TriaP1Input(SurfaceEnum,newsurface));
1602 this->inputs->AddInput(new TriaP1Input(BedEnum,newbed));
1603
1604 /*Free ressources:*/
1605 xfree((void**)&doflist);
1606}
1607/*}}}*/
1608/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
1609void Tria::InputUpdateFromVector(double* vector, int name, int type){
1610
1611 /*Check that name is an element input*/
1612 if (!IsInput(name)) return;
1613
1614 switch(type){
1615
1616 case VertexEnum:
1617
1618 /*New TriaP1Input*/
1619 double values[3];
1620
1621 /*Get values on the 3 vertices*/
1622 for (int i=0;i<3;i++){
1623 values[i]=vector[this->nodes[i]->GetVertexDof()];
1624 }
1625
1626 /*update input*/
1627 if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
1628 matice->inputs->AddInput(new TriaP1Input(name,values));
1629 }
1630 else{
1631 this->inputs->AddInput(new TriaP1Input(name,values));
1632 }
1633 return;
1634
1635 default:
1636 _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
1637 }
1638}
1639/*}}}*/
1640/*FUNCTION Tria::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
1641void Tria::InputUpdateFromVector(int* vector, int name, int type){
1642 _error_(" not supported yet!");
1643}
1644/*}}}*/
1645/*FUNCTION Tria::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
1646void Tria::InputUpdateFromVector(bool* vector, int name, int type){
1647 _error_(" not supported yet!");
1648}
1649/*}}}*/
1650/*FUNCTION Tria::InputCreate(double scalar,int enum,int code);{{{1*/
1651void Tria::InputCreate(double scalar,int name,int code){
1652
1653 /*Check that name is an element input*/
1654 if (!IsInput(name)) return;
1655
1656 if ((code==5) || (code==1)){ //boolean
1657 this->inputs->AddInput(new BoolInput(name,(bool)scalar));
1658 }
1659 else if ((code==6) || (code==2)){ //integer
1660 this->inputs->AddInput(new IntInput(name,(int)scalar));
1661 }
1662 else if ((code==7) || (code==3)){ //double
1663 this->inputs->AddInput(new DoubleInput(name,(int)scalar));
1664 }
1665 else _error_("%s%i"," could not recognize nature of vector from code ",code);
1666
1667}
1668/*}}}*/
1669/*FUNCTION Tria::InputCreate(double* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{1*/
1670void Tria::InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements
1671
1672 /*Intermediaries*/
1673 int i,j,t;
1674 int tria_vertex_ids[3];
1675 int row;
1676 double nodeinputs[3];
1677 double time;
1678 TransientInput* transientinput=NULL;
1679 int numberofvertices;
1680 int numberofelements;
1681 double yts;
1682
1683
1684 /*Fetch parameters: */
1685 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
1686 iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
1687 iomodel->Constant(&yts,ConstantsYtsEnum);
1688
1689 /*Branch on type of vector: nodal or elementary: */
1690 if(vector_type==1){ //nodal vector
1691
1692 /*Recover vertices ids needed to initialize inputs*/
1693 for(i=0;i<3;i++){
1694 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
1695 }
1696
1697 /*Are we in transient or static? */
1698 if(M==numberofvertices){
1699
1700 /*create input values: */
1701 for(i=0;i<3;i++)nodeinputs[i]=(double)vector[tria_vertex_ids[i]-1];
1702
1703 /*process units: */
1704 UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
1705
1706 /*create static input: */
1707 this->inputs->AddInput(new TriaP1Input(vector_enum,nodeinputs));
1708 }
1709 else if(M==numberofvertices+1){
1710 /*create transient input: */
1711 for(t=0;t<N;t++){ //N is the number of times
1712
1713 /*create input values: */
1714 for(i=0;i<3;i++){
1715 row=tria_vertex_ids[i]-1;
1716 nodeinputs[i]=(double)vector[N*row+t];
1717 }
1718
1719 /*process units: */
1720 UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
1721
1722 /*time? :*/
1723 time=(double)vector[(M-1)*N+t]*yts;
1724
1725 if(t==0) transientinput=new TransientInput(vector_enum);
1726 transientinput->AddTimeInput(new TriaP1Input(vector_enum,nodeinputs),time);
1727 }
1728 this->inputs->AddInput(transientinput);
1729 }
1730 else _error_("nodal vector is either numberofnodes or numberofnodes+1 long. Field provided (%s) is %i long",EnumToStringx(vector_enum),M);
1731 }
1732 else if(vector_type==2){ //element vector
1733 /*Are we in transient or static? */
1734 if(M==numberofelements){
1735
1736 /*static mode: create an input out of the element value: */
1737
1738 if (code==5){ //boolean
1739 this->inputs->AddInput(new BoolInput(vector_enum,(bool)vector[index]));
1740 }
1741 else if (code==6){ //integer
1742 this->inputs->AddInput(new IntInput(vector_enum,(int)vector[index]));
1743 }
1744 else if (code==7){ //double
1745 this->inputs->AddInput(new DoubleInput(vector_enum,(double)vector[index]));
1746 }
1747 else _error_("%s%i"," could not recognize nature of vector from code ",code);
1748 }
1749 else {
1750 _error_("transient elementary inputs not supported yet!");
1751 }
1752 }
1753 else{
1754 _error_("Cannot add input for vector type %i (not supported)",vector_type);
1755 }
1756
1757}
1758/*}}}*/
1759/*FUNCTION Tria::IsInput{{{1*/
1760bool Tria::IsInput(int name){
1761 if (
1762 name==ThicknessEnum ||
1763 name==SurfaceEnum ||
1764 name==BedEnum ||
1765 name==SurfaceSlopeXEnum ||
1766 name==SurfaceSlopeYEnum ||
1767 name==BasalforcingsMeltingRateEnum ||
1768 name==WatercolumnEnum ||
1769 name==SurfaceforcingsMassBalanceEnum ||
1770 name==SurfaceAreaEnum||
1771 name==VxEnum ||
1772 name==VyEnum ||
1773 name==InversionVxObsEnum ||
1774 name==InversionVyObsEnum ||
1775 name==FrictionCoefficientEnum ||
1776 name==MaterialsRheologyBbarEnum ||
1777 name==GradientEnum ||
1778 name==OldGradientEnum ||
1779 name==QmuVxEnum ||
1780 name==QmuVyEnum ||
1781 name==QmuPressureEnum ||
1782 name==QmuBedEnum ||
1783 name==QmuThicknessEnum ||
1784 name==QmuSurfaceEnum ||
1785 name==QmuTemperatureEnum ||
1786 name==QmuMeltingEnum
1787 ){
1788 return true;
1789 }
1790 else return false;
1791}
1792/*}}}*/
1793/*FUNCTION Tria::IsOnBed {{{1*/
1794bool Tria::IsOnBed(){
1795
1796 bool onbed;
1797 inputs->GetInputValue(&onbed,MeshElementonbedEnum);
1798 return onbed;
1799}
1800/*}}}*/
1801/*FUNCTION Tria::IsFloating {{{1*/
1802bool Tria::IsFloating(){
1803
1804 bool shelf;
1805 inputs->GetInputValue(&shelf,MaskElementonfloatingiceEnum);
1806 return shelf;
1807}
1808/*}}}*/
1809/*FUNCTION Tria::IsNodeOnShelf {{{1*/
1810bool Tria::IsNodeOnShelf(){
1811
1812 int i;
1813 bool shelf=false;
1814
1815 for(i=0;i<3;i++){
1816 if (nodes[i]->IsFloating()){
1817 shelf=true;
1818 break;
1819 }
1820 }
1821 return shelf;
1822}
1823/*}}}*/
1824/*FUNCTION Tria::IsNodeOnShelfFromFlags {{{1*/
1825bool Tria::IsNodeOnShelfFromFlags(double* flags){
1826
1827 int i;
1828 bool shelf=false;
1829
1830 for(i=0;i<NUMVERTICES;i++){
1831 if (flags[nodes[i]->Sid()]){
1832 shelf=true;
1833 break;
1834 }
1835 }
1836 return shelf;
1837}
1838/*}}}*/
1839/*FUNCTION Tria::IsOnWater {{{1*/
1840bool Tria::IsOnWater(){
1841
1842 bool water;
1843 inputs->GetInputValue(&water,MaskElementonwaterEnum);
1844 return water;
1845}
1846/*}}}*/
1847/*FUNCTION Tria::ListResultsInfo{{{*/
1848void Tria::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,double** in_resultstimes,int** in_resultssteps,int* in_num_results){
1849
1850 /*Intermediaries*/
1851 int i;
1852 int numberofresults = 0;
1853 int *resultsenums = NULL;
1854 int *resultssizes = NULL;
1855 double *resultstimes = NULL;
1856 int *resultssteps = NULL;
1857
1858 /*Checks*/
1859 _assert_(in_num_results);
1860
1861 /*Count number of results*/
1862 for(i=0;i<this->results->Size();i++){
1863 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
1864 numberofresults++;
1865 }
1866
1867 if(numberofresults){
1868
1869 /*Allocate output*/
1870 resultsenums=(int*)xmalloc(numberofresults*sizeof(int));
1871 resultssizes=(int*)xmalloc(numberofresults*sizeof(int));
1872 resultstimes=(double*)xmalloc(numberofresults*sizeof(double));
1873 resultssteps=(int*)xmalloc(numberofresults*sizeof(int));
1874
1875 /*populate enums*/
1876 for(i=0;i<this->results->Size();i++){
1877 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
1878 resultsenums[i]=elementresult->InstanceEnum();
1879 resultstimes[i]=elementresult->GetTime();
1880 resultssteps[i]=elementresult->GetStep();
1881 if(elementresult->ObjectEnum()==TriaP1ElementResultEnum){
1882 resultssizes[i]=P1Enum;
1883 }
1884 else{
1885 resultssizes[i]=P0Enum;
1886 }
1887 }
1888 }
1889
1890 /*Assign output pointers:*/
1891 *in_num_results=numberofresults;
1892 *in_resultsenums=resultsenums;
1893 *in_resultssizes=resultssizes;
1894 *in_resultstimes=resultstimes;
1895 *in_resultssteps=resultssteps;
1896
1897}/*}}}*/
1898/*FUNCTION Tria::MigrateGroundingLine{{{1*/
1899void Tria::MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding){
1900
1901 int i,migration_style,unground;
1902 bool elementonshelf = false;
1903 double bed_hydro,yts,gl_melting_rate;
1904 double rho_water,rho_ice,density;
1905 double melting[NUMVERTICES];
1906 double h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
1907
1908 /*Recover info at the vertices: */
1909 parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
1910 parameters->FindParam(&yts,ConstantsYtsEnum);
1911 GetInputListOnVertices(&h[0],ThicknessEnum);
1912 GetInputListOnVertices(&s[0],SurfaceEnum);
1913 GetInputListOnVertices(&b[0],BedEnum);
1914 GetInputListOnVertices(&ba[0],BathymetryEnum);
1915 rho_water=matpar->GetRhoWater();
1916 rho_ice=matpar->GetRhoIce();
1917 density=rho_ice/rho_water;
1918
1919 /*go through vertices, and update inputs, considering them to be TriaVertex type: */
1920 for(i=0;i<NUMVERTICES;i++){
1921 /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
1922 if(old_floating_ice[nodes[i]->Sid()]){
1923 if(b[i]<=ba[i]){
1924 b[i]=ba[i];
1925 s[i]=b[i]+h[i];
1926 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
1927 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
1928 }
1929 }
1930 /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
1931 /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
1932 else{
1933 bed_hydro=-density*h[i];
1934 if (bed_hydro>ba[i]){
1935 /*Unground only if the element is connected to the ice shelf*/
1936 if(migration_style==AgressiveMigrationEnum){
1937 s[i]=(1-density)*h[i];
1938 b[i]=-density*h[i];
1939 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
1940 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
1941 }
1942 else if(migration_style==SoftMigrationEnum && sheet_ungrounding[nodes[i]->Sid()]){
1943 s[i]=(1-density)*h[i];
1944 b[i]=-density*h[i];
1945 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
1946 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
1947 }
1948 }
1949 }
1950 }
1951
1952 /*If at least one vertex is now floating, the element is now floating*/
1953 for(i=0;i<NUMVERTICES;i++){
1954 if(nodes[i]->IsFloating()){
1955 elementonshelf=true;
1956 break;
1957 }
1958 }
1959
1960 /*Add basal melting rate if element just ungrounded*/
1961 if(!this->IsFloating() && elementonshelf==true){
1962 for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
1963 this->inputs->AddInput(new TriaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
1964 }
1965
1966 /*Update inputs*/
1967 this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
1968
1969 /*Update inputs*/
1970 this->inputs->AddInput(new TriaP1Input(SurfaceEnum,&s[0]));
1971 this->inputs->AddInput(new TriaP1Input(BedEnum,&b[0]));
1972}
1973/*}}}*/
1974/*FUNCTION Tria::MyRank {{{1*/
1975int Tria::MyRank(void){
1976 extern int my_rank;
1977 return my_rank;
1978}
1979/*}}}*/
1980/*FUNCTION Tria::NodalValue {{{1*/
1981int Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){
1982
1983 int i;
1984 int found=0;
1985 double value;
1986 Input* data=NULL;
1987 GaussTria *gauss = NULL;
1988
1989 /*First, serarch the input: */
1990 data=inputs->GetInput(natureofdataenum);
1991
1992 /*figure out if we have the vertex id: */
1993 found=0;
1994 for(i=0;i<NUMVERTICES;i++){
1995 if(index==nodes[i]->GetVertexId()){
1996 /*Do we have natureofdataenum in our inputs? :*/
1997 if(data){
1998 /*ok, we are good. retrieve value of input at vertex :*/
1999 gauss=new GaussTria(); gauss->GaussVertex(i);
2000 data->GetInputValue(&value,gauss);
2001 found=1;
2002 break;
2003 }
2004 }
2005 }
2006
2007 if(found)*pvalue=value;
2008 return found;
2009}
2010/*}}}*/
2011/*FUNCTION Tria::PatchFill{{{1*/
2012void Tria::PatchFill(int* prow, Patch* patch){
2013
2014 int i,row;
2015 int vertices_ids[3];
2016
2017 /*recover pointer: */
2018 row=*prow;
2019
2020 for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.
2021
2022 for(i=0;i<this->results->Size();i++){
2023 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
2024
2025 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand
2026 *it to the result object, to fill the rest: */
2027 patch->fillelementinfo(row,this->sid+1,vertices_ids,3);
2028 elementresult->PatchFill(row,patch);
2029
2030 /*increment rower: */
2031 row++;
2032 }
2033
2034 /*Assign output pointers:*/
2035 *prow=row;
2036}
2037/*}}}*/
2038/*FUNCTION Tria::PatchSize{{{1*/
2039void Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
2040
2041 int i;
2042 int numrows = 0;
2043 int numnodes = 0;
2044 int temp_numnodes = 0;
2045
2046 /*Go through all the results objects, and update the counters: */
2047 for (i=0;i<this->results->Size();i++){
2048 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
2049 /*first, we have one more result: */
2050 numrows++;
2051 /*now, how many vertices and how many nodal values for this result? :*/
2052 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
2053 if(temp_numnodes>numnodes)numnodes=temp_numnodes;
2054 }
2055
2056 /*Assign output pointers:*/
2057 *pnumrows=numrows;
2058 *pnumvertices=NUMVERTICES;
2059 *pnumnodes=numnodes;
2060}
2061/*}}}*/
2062/*FUNCTION Tria::PotentialSheetUngrounding{{{1*/
2063void Tria::PotentialSheetUngrounding(Vector* potential_sheet_ungrounding){
2064
2065 int i;
2066 double h[NUMVERTICES],ba[NUMVERTICES];
2067 double bed_hydro;
2068 double rho_water,rho_ice,density;
2069 bool elementonshelf = false;
2070
2071 /*material parameters: */
2072 rho_water=matpar->GetRhoWater();
2073 rho_ice=matpar->GetRhoIce();
2074 density=rho_ice/rho_water;
2075 GetInputListOnVertices(&h[0],ThicknessEnum);
2076 GetInputListOnVertices(&ba[0],BathymetryEnum);
2077
2078 /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
2079 for(i=0;i<NUMVERTICES;i++){
2080 /*Find if grounded vertices want to start floating*/
2081 if (!nodes[i]->IsFloating()){
2082 bed_hydro=-density*h[i];
2083 if (bed_hydro>ba[i]){
2084 /*Vertex that could potentially unground, flag it*/
2085 potential_sheet_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
2086 }
2087 }
2088 }
2089}
2090/*}}}*/
2091/*FUNCTION Tria::PositiveDegreeDay{{{1*/
2092void Tria::PositiveDegreeDay(double* pdds,double* pds,double signorm){
2093
2094 int i,iqj,imonth;
2095 double agd[NUMVERTICES]; // surface and basal
2096 double saccu[NUMVERTICES] = {0}; // yearly surface accumulation
2097 double smelt[NUMVERTICES] = {0}; // yearly melt
2098 double precrunoff[NUMVERTICES]; // yearly runoff
2099 double prect; // total precipitation during 1 year taking into account des. ef.
2100 double water; //water=rain + snowmelt
2101 double runoff; //meltwater only, does not include rain
2102 double sconv; //rhow_rain/rhoi / 12 months
2103
2104 double rho_water,rho_ice,density;
2105 double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter.
2106 double desfac = 0.5; //desert elevation factor
2107 double s0p[NUMVERTICES]={0}; //should be set to elevation from precip source
2108 double s0t[NUMVERTICES]={0}; //should be set to elevation from temperature source
2109 double st; // elevation between altitude of the temp record and current altitude
2110 double sp; // elevation between altitude of the prec record and current altitude
2111
2112
2113 // PDD and PD constants and variables
2114 double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
2115 double siglimc=0, siglim0, siglim0c;
2116 double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
2117 double DT = 0.02;
2118 double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
2119
2120 double q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
2121 double qm[NUMVERTICES] = {0}; // snow part of the precipitation
2122 double qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment
2123 double qmp[NUMVERTICES] = {0}; // desertification taken into account
2124 double pdd[NUMVERTICES] = {0};
2125 double frzndd[NUMVERTICES] = {0};
2126
2127 double tstar; // monthly mean surface temp
2128 double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature
2129 double Tsurf[NUMVERTICES] = {0}; // average annual temperature
2130
2131 double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
2132 double t[NUMVERTICES][12],prec[NUMVERTICES][12];
2133 double deltm=1/12;
2134 int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
2135
2136 double snwm; // snow that could have been melted in a year.
2137 double snwmf; // ablation factor for snow per positive degree day.
2138 double smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
2139
2140 double dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
2141 double supice,supcap,diffndd;
2142 double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
2143 double pddtj[NUMVERTICES], hmx2;
2144
2145 /*Recover info at the vertices: */
2146 GetInputListOnVertices(&h[0],ThicknessEnum);
2147 GetInputListOnVertices(&s[0],SurfaceEnum);
2148 GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum);
2149 GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
2150
2151 for(i=0;i<NUMVERTICES;i++) ttmp[i]=ttmp[i]-273.15; // convertion from Kelvin to celcius
2152
2153 for(i=0;i<NUMVERTICES;i++)
2154 for(imonth=0;imonth<12;imonth++){
2155 t[i][imonth]=ttmp[i];
2156 prec[i][imonth]=prectmp[i];
2157 }
2158
2159 /*Get material parameters :*/
2160 rho_ice=matpar->GetRhoIce();
2161 rho_water=matpar->GetRhoWater();
2162 density=rho_ice/rho_water;
2163
2164 sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
2165
2166 /*PDD constant*/
2167 siglim = 2.5*signorm;
2168 siglim0 = siglim/DT + 0.5;
2169 siglim0c = siglimc/DT + 0.5;
2170 PDup = siglimc+PDCUT;
2171
2172 // seasonal loop
2173 for (iqj = 0; iqj < 12; iqj++){
2174 imonth = ismon[iqj];
2175 for (i = 0; i < NUMVERTICES; i++){
2176 st=(s[i]-s0t[i])/1000;
2177 tstar = t[i][imonth] - lapser *max(st,sealev);
2178 Tsurf[i] = tstar*deltm+Tsurf[i];
2179
2180 /*********compute PD ****************/
2181 if (tstar < PDup){
2182 pd = 1;
2183 if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
2184 else {
2185 pd = 0;}
2186
2187 /******exp des/elev precip reduction*******/
2188 sp=(s[i]-s0p[i])/1000; // deselev effect is wrt chng in topo
2189 if (sp>0.0){q = exp(-desfac*sp);}
2190 else {q = 1.0;}
2191
2192 qmt[i]= qmt[i] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent
2193 qmpt= q*prec[i][imonth]*sconv;
2194 qmp[i]= qmp[i] + qmpt;
2195 qm[i]= qm[i] + qmpt*pd;
2196
2197 /*********compute PDD************/
2198 // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
2199 // gaussian=T_m, so ndd=-(Tsurf-pdd)
2200 if (iqj>6 && iqj<10){ Tsum[i]=Tsum[i]+tstar;}
2201 if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
2202 else if (tstar> -siglim){
2203 pddsig=pdds[int(tstar/DT + siglim0)];
2204 pdd[i] = pdd[i] + pddsig*deltm;
2205 frzndd[i] = frzndd[i] - (tstar-pddsig)*deltm;}
2206 else{frzndd[i] = frzndd[i] - tstar*deltm; }
2207 }
2208 } // end of seasonal loop
2209
2210 //******************************************************************
2211 for(i=0;i<NUMVERTICES;i++){
2212 saccu[i] = qm[i];
2213 prect = qmp[i]; // total precipitation during 1 year taking into account des. ef.
2214 Tsum[i]=Tsum[i]/3;
2215
2216 /***** determine PDD factors *****/
2217 if(Tsum[i]< -1.) {
2218 snwmf=2.65*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
2219 smf=17.22*0.001; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
2220 }
2221 else if(Tsum[i]< 10){
2222 snwmf = (0.15*Tsum[i] + 2.8)*0.001;
2223 smf = (0.0067*pow((10.-Tsum[i]),3) + 8.3)*0.001;
2224 }
2225 else{
2226 snwmf=4.3*0.001;
2227 smf=8.3*0.001;
2228 }
2229 snwmf=0.95*snwmf;
2230 smf=0.95*smf;
2231
2232 /***** compute PDD ablation and refreezing *****/
2233 pddt = pdd[i] *365;
2234 snwm = snwmf*pddt; // snow that could have been melted in a year
2235 hmx2 = min(h[i],dfrz); // refreeze active layer max depth: dfrz
2236
2237 if(snwm < saccu[i]) {
2238 water=prect-saccu[i] + snwm; //water=rain + snowmelt
2239 // l 2.2= capillary factor
2240 // Should refreezing be controlled by frzndd or by mean annual Tsurf?
2241 // dCovrLm concept is of warming of active layer (thickness =d@=1-
2242 // >2m)
2243 // problem with water seepage into ice: should be sealed after
2244 // refreezing
2245 // so everything needs to be predicated on 1 year scale, except for
2246 // thermal
2247 // conductivity through ice
2248 // also, need to account that melt season has low accum, so what's
2249 // going to
2250 // hold the meltwater around for refreezing? And melt-time will have
2251 // low seasonal frzndd
2252
2253 // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002
2254
2255 supice= min(hmx2*CovrLm*frzndd[i]+2.2*(saccu[i]-snwm), water); // superimposed ice
2256 supcap=min(2.2*(saccu[i]-snwm),water);
2257 runoff=snwm - supice; //meltwater only, does not include rain
2258 }
2259 else { //all snow melted
2260 supice= min(hmx2*CovrLm*frzndd[i], prect );
2261 runoff= saccu[i] + smf*(pddt-saccu[i]/snwmf) - supice;
2262 supcap=0;
2263 }
2264 // pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
2265 // except pdd melt heat source is atmosphere, while refreeze is
2266 // ground/ice stored interim
2267 // assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
2268 // assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
2269 // <0
2270 // assume ndd>pdd, little melt => little supice
2271 // bottom line: compare for Tsurf<0 : supice and no supice case,
2272 // expect Tsurf difference
2273 // except some of cooling flux comes from atmosphere//
2274 // 1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
2275 // does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
2276 // < 0.1
2277
2278 // make more sense to just use residual pdd-ndd except that pdd
2279 // residual not clear yet
2280 // frzndd should not be used up by refreezing in snow, so stick in
2281 // supcap.
2282 diffndd=0;
2283 if (frzndd[i]>0) {
2284 diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd[i]);
2285 frzndd[i]=frzndd[i]-diffndd;
2286 }
2287 if(runoff<0){
2288 saccu[i]= saccu[i] -runoff;
2289 smelt[i] = 0;
2290 precrunoff[i]=prect-saccu[i];
2291 //here assume pdd residual is 0, =>
2292 Tsurf[i]= max(Tsurf[i],-frzndd[i]);
2293 }
2294 else {
2295 smelt[i] = runoff;
2296 precrunoff[i]=prect-max(0.,supice)-saccu[i];}
2297 //here really need pdd balance, try 0.5 fudge factor?
2298 //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
2299 //yet from site plots, can be ice free with Tsurf=-5.5C
2300 if(Tsurf[i]<0) {
2301 Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
2302
2303 agd[i] = -smelt[i]+saccu[i];
2304 pddtj[i]=pddt;
2305
2306 /*Update inputs*/
2307 this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
2308 // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
2309 this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
2310
2311 } //end of the for loop over the vertices
2312}
2313/*}}}*/
2314/*FUNCTION Tria::ProcessResultsUnits{{{1*/
2315void Tria::ProcessResultsUnits(void){
2316
2317 int i;
2318
2319 for(i=0;i<this->results->Size();i++){
2320 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
2321 elementresult->ProcessUnits(this->parameters);
2322 }
2323}
2324/*}}}*/
2325/*FUNCTION Tria::RequestedOutput{{{1*/
2326void Tria::RequestedOutput(int output_enum,int step,double time){
2327
2328 if(IsInput(output_enum)){
2329 /*just transfer this input to results, and we are done: */
2330 InputToResult(output_enum,step,time);
2331 }
2332 else{
2333 /*this input does not exist, compute it, and then transfer to results: */
2334 switch(output_enum){
2335 case StressTensorEnum:
2336 this->ComputeStressTensor();
2337 InputToResult(StressTensorxxEnum,step,time);
2338 InputToResult(StressTensorxyEnum,step,time);
2339 InputToResult(StressTensorxzEnum,step,time);
2340 InputToResult(StressTensoryyEnum,step,time);
2341 InputToResult(StressTensoryzEnum,step,time);
2342 InputToResult(StressTensorzzEnum,step,time);
2343 break;
2344
2345 default:
2346 /*do nothing, no need to derail the computation because one of the outputs requested cannot be found: */
2347 break;
2348 }
2349 }
2350
2351}
2352/*}}}*/
2353/*FUNCTION Tria::SetClone {{{1*/
2354void Tria::SetClone(int* minranks){
2355
2356 _error_("not implemented yet");
2357}
2358/*}}}1*/
2359/*FUNCTION Tria::SmearFunction {{{1*/
2360void Tria::SmearFunction(Vector* smearedvector,double (*WeightFunction)(double distance,double radius),double radius){
2361 _error_("not implemented yet");
2362
2363}
2364/*}}}1*/
2365/*FUNCTION Tria::SetCurrentConfiguration {{{1*/
2366void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
2367
2368 /*go into parameters and get the analysis_counter: */
2369 int analysis_counter;
2370 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
2371
2372 /*Get Element type*/
2373 this->element_type=this->element_type_list[analysis_counter];
2374
2375 /*Pick up nodes*/
2376 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
2377 else this->nodes=NULL;
2378
2379}
2380/*}}}*/
2381/*FUNCTION Tria::SurfaceArea {{{1*/
2382double Tria::SurfaceArea(void){
2383
2384 int i;
2385 double S;
2386 double normal[3];
2387 double v13[3],v23[3];
2388 double xyz_list[NUMVERTICES][3];
2389
2390 /*If on water, return 0: */
2391 if(IsOnWater())return 0;
2392
2393 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2394
2395 for (i=0;i<3;i++){
2396 v13[i]=xyz_list[0][i]-xyz_list[2][i];
2397 v23[i]=xyz_list[1][i]-xyz_list[2][i];
2398 }
2399
2400 normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
2401 normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
2402 normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
2403
2404 S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
2405
2406 /*Return: */
2407 return S;
2408}
2409/*}}}*/
2410/*FUNCTION Tria::SurfaceNormal{{{1*/
2411void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
2412
2413 int i;
2414 double v13[3],v23[3];
2415 double normal[3];
2416 double normal_norm;
2417
2418 for (i=0;i<3;i++){
2419 v13[i]=xyz_list[0][i]-xyz_list[2][i];
2420 v23[i]=xyz_list[1][i]-xyz_list[2][i];
2421 }
2422
2423 normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
2424 normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
2425 normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
2426
2427 normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
2428
2429 *(surface_normal)=normal[0]/normal_norm;
2430 *(surface_normal+1)=normal[1]/normal_norm;
2431 *(surface_normal+2)=normal[2]/normal_norm;
2432}
2433/*}}}*/
2434/*FUNCTION Tria::TimeAdapt{{{1*/
2435double Tria::TimeAdapt(void){
2436
2437 /*intermediary: */
2438 int i;
2439 double C,dt;
2440 double dx,dy;
2441 double maxx,minx;
2442 double maxy,miny;
2443 double maxabsvx,maxabsvy;
2444 double xyz_list[NUMVERTICES][3];
2445
2446 /*get CFL coefficient:*/
2447 this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum);
2448
2449 /*Get for Vx and Vy, the max of abs value: */
2450 #ifdef _HAVE_RESPONSES_
2451 this->MaxAbsVx(&maxabsvx,false);
2452 this->MaxAbsVy(&maxabsvy,false);
2453 #else
2454 _error_("ISSM was not compiled with responses compiled in, exiting!");
2455 #endif
2456
2457 /* Get node coordinates and dof list: */
2458 GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES);
2459
2460 minx=xyz_list[0][0];
2461 maxx=xyz_list[0][0];
2462 miny=xyz_list[0][1];
2463 maxy=xyz_list[0][1];
2464
2465 for(i=1;i<NUMVERTICES;i++){
2466 if (xyz_list[i][0]<minx)minx=xyz_list[i][0];
2467 if (xyz_list[i][0]>maxx)maxx=xyz_list[i][0];
2468 if (xyz_list[i][1]<miny)miny=xyz_list[i][1];
2469 if (xyz_list[i][1]>maxy)maxy=xyz_list[i][1];
2470 }
2471 dx=maxx-minx;
2472 dy=maxy-miny;
2473
2474 /*CFL criterion: */
2475 dt=C/(maxabsvy/dx+maxabsvy/dy);
2476
2477 return dt;
2478}
2479/*}}}*/
2480/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/
2481void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index
2482
2483 /*Intermediaries*/
2484 int i,j;
2485 int tria_node_ids[3];
2486 int tria_vertex_ids[3];
2487 int tria_type;
2488 double nodeinputs[3];
2489 double yts;
2490 int progstabilization,balancestabilization;
2491 bool dakota_analysis;
2492
2493 /*Checks if debuging*/
2494 /*{{{2*/
2495 _assert_(iomodel->Data(MeshElementsEnum));
2496 /*}}}*/
2497
2498 /*Fetch parameters: */
2499 iomodel->Constant(&yts,ConstantsYtsEnum);
2500 iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
2501 iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
2502 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
2503
2504 /*Recover element type*/
2505 if ((analysis_type==PrognosticAnalysisEnum && progstabilization==3) || (analysis_type==BalancethicknessAnalysisEnum && balancestabilization==3)){
2506 /*P1 Discontinuous Galerkin*/
2507 tria_type=P1DGEnum;
2508 }
2509 else{
2510 /*P1 Continuous Galerkin*/
2511 tria_type=P1Enum;
2512 }
2513 this->SetElementType(tria_type,analysis_counter);
2514
2515 /*Recover vertices ids needed to initialize inputs*/
2516 for(i=0;i<3;i++){
2517 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
2518 }
2519
2520 /*Recover nodes ids needed to initialize the node hook.*/
2521 if (tria_type==P1DGEnum){
2522 /*Discontinuous Galerkin*/
2523 tria_node_ids[0]=iomodel->nodecounter+3*index+1;
2524 tria_node_ids[1]=iomodel->nodecounter+3*index+2;
2525 tria_node_ids[2]=iomodel->nodecounter+3*index+3;
2526 }
2527 else{
2528 /*Continuous Galerkin*/
2529 for(i=0;i<3;i++){
2530 tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->Data(MeshElementsEnum)+3*index+i); //ids for vertices are in the elements array from Matlab
2531 }
2532 }
2533
2534 /*hooks: */
2535 this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
2536
2537 /*Fill with IoModel*/
2538 this->InputUpdateFromIoModel(index,iomodel);
2539
2540 /*Defaults if not provided in iomodel*/
2541 switch(analysis_type){
2542
2543 case DiagnosticHorizAnalysisEnum:
2544
2545 /*default vx,vy and vz: either observation or 0 */
2546 if(!iomodel->Data(VxEnum)){
2547 for(i=0;i<3;i++)nodeinputs[i]=0;
2548 this->inputs->AddInput(new TriaP1Input(VxEnum,nodeinputs));
2549 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVxEnum,nodeinputs));
2550 }
2551 if(!iomodel->Data(VyEnum)){
2552 for(i=0;i<3;i++)nodeinputs[i]=0;
2553 this->inputs->AddInput(new TriaP1Input(VyEnum,nodeinputs));
2554 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVyEnum,nodeinputs));
2555 }
2556 if(!iomodel->Data(VzEnum)){
2557 for(i=0;i<3;i++)nodeinputs[i]=0;
2558 this->inputs->AddInput(new TriaP1Input(VzEnum,nodeinputs));
2559 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVzEnum,nodeinputs));
2560 }
2561 if(!iomodel->Data(PressureEnum)){
2562 for(i=0;i<3;i++)nodeinputs[i]=0;
2563 if(dakota_analysis){
2564 this->inputs->AddInput(new TriaP1Input(PressureEnum,nodeinputs));
2565 this->inputs->AddInput(new TriaP1Input(QmuPressureEnum,nodeinputs));
2566 }
2567 }
2568 break;
2569
2570 default:
2571 /*No update for other solution types*/
2572 break;
2573
2574 }
2575
2576 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
2577 this->parameters=NULL;
2578}
2579/*}}}*/
2580/*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{1*/
2581int Tria::UpdatePotentialSheetUngrounding(double* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,double* nodes_on_iceshelf){
2582
2583 int i;
2584 int nflipped=0;
2585
2586 /*Go through nodes, and whoever is on the potential_sheet_ungrounding, ends up in nodes_on_iceshelf: */
2587 for(i=0;i<3;i++){
2588 if (vertices_potentially_ungrounding[nodes[i]->Sid()]){
2589 vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
2590
2591 /*If node was not on ice shelf, we flipped*/
2592 if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
2593 nflipped++;
2594 }
2595 }
2596 }
2597 return nflipped;
2598}
2599/*}}}*/
2600
2601#ifdef _HAVE_RESPONSES_
2602/*FUNCTION Tria::IceVolume {{{1*/
2603double Tria::IceVolume(void){
2604
2605 /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
2606 double base,surface,bed;
2607 double xyz_list[NUMVERTICES][3];
2608
2609 if(IsOnWater())return 0;
2610
2611 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2612
2613 /*First calculate the area of the base (cross section triangle)
2614 * http://en.wikipedia.org/wiki/Triangle
2615 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
2616 base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
2617
2618 /*Now get the average height*/
2619 Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
2620 Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
2621 surface_input->GetInputAverage(&surface);
2622 bed_input->GetInputAverage(&bed);
2623
2624 /*Return: */
2625 return base*(surface-bed);
2626}
2627/*}}}*/
2628/*FUNCTION Tria::MassFlux {{{1*/
2629double Tria::MassFlux( double* segment,bool process_units){
2630
2631 const int numdofs=2;
2632
2633 int i,dim;
2634 double mass_flux=0;
2635 double xyz_list[NUMVERTICES][3];
2636 double normal[2];
2637 double length,rho_ice;
2638 double x1,y1,x2,y2,h1,h2;
2639 double vx1,vx2,vy1,vy2;
2640 GaussTria* gauss_1=NULL;
2641 GaussTria* gauss_2=NULL;
2642
2643 /*Get material parameters :*/
2644 rho_ice=matpar->GetRhoIce();
2645
2646 /*First off, check that this segment belongs to this element: */
2647 if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id);
2648
2649 /*Recover segment node locations: */
2650 x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
2651
2652 /*Get xyz list: */
2653 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2654
2655 /*get area coordinates of 0 and 1 locations: */
2656 gauss_1=new GaussTria();
2657 gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
2658 gauss_2=new GaussTria();
2659 gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
2660
2661 normal[0]=cos(atan2(x1-x2,y2-y1));
2662 normal[1]=sin(atan2(x1-x2,y2-y1));
2663
2664 length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
2665
2666 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
2667 this->parameters->FindParam(&dim,MeshDimensionEnum);
2668 Input* vx_input=NULL;
2669 Input* vy_input=NULL;
2670 if(dim==2){
2671 vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2672 vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2673 }
2674 else{
2675 vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
2676 vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
2677 }
2678
2679 thickness_input->GetInputValue(&h1, gauss_1);
2680 thickness_input->GetInputValue(&h2, gauss_2);
2681 vx_input->GetInputValue(&vx1,gauss_1);
2682 vx_input->GetInputValue(&vx2,gauss_2);
2683 vy_input->GetInputValue(&vy1,gauss_1);
2684 vy_input->GetInputValue(&vy2,gauss_2);
2685
2686 mass_flux= rho_ice*length*(
2687 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
2688 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
2689 );
2690
2691 /*Process units: */
2692 mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum);
2693
2694 /*clean up and return:*/
2695 delete gauss_1;
2696 delete gauss_2;
2697 return mass_flux;
2698}
2699/*}}}*/
2700/*FUNCTION Tria::MaxAbsVx{{{1*/
2701void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){
2702
2703 /*Get maximum:*/
2704 double maxabsvx=this->inputs->MaxAbs(VxEnum);
2705
2706 /*process units if requested: */
2707 if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum);
2708
2709 /*Assign output pointers:*/
2710 *pmaxabsvx=maxabsvx;
2711}
2712/*}}}*/
2713/*FUNCTION Tria::MaxAbsVy{{{1*/
2714void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){
2715
2716 /*Get maximum:*/
2717 double maxabsvy=this->inputs->MaxAbs(VyEnum);
2718
2719 /*process units if requested: */
2720 if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum);
2721
2722 /*Assign output pointers:*/
2723 *pmaxabsvy=maxabsvy;
2724}
2725/*}}}*/
2726/*FUNCTION Tria::MaxAbsVz{{{1*/
2727void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){
2728
2729 /*Get maximum:*/
2730 double maxabsvz=this->inputs->MaxAbs(VzEnum);
2731
2732 /*process units if requested: */
2733 if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum);
2734
2735 /*Assign output pointers:*/
2736 *pmaxabsvz=maxabsvz;
2737}
2738/*}}}*/
2739/*FUNCTION Tria::MaxVel{{{1*/
2740void Tria::MaxVel(double* pmaxvel, bool process_units){
2741
2742 /*Get maximum:*/
2743 double maxvel=this->inputs->Max(VelEnum);
2744
2745 /*process units if requested: */
2746 if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum);
2747
2748 /*Assign output pointers:*/
2749 *pmaxvel=maxvel;
2750}
2751/*}}}*/
2752/*FUNCTION Tria::MaxVx{{{1*/
2753void Tria::MaxVx(double* pmaxvx, bool process_units){
2754
2755 /*Get maximum:*/
2756 double maxvx=this->inputs->Max(VxEnum);
2757
2758 /*process units if requested: */
2759 if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum);
2760
2761 /*Assign output pointers:*/
2762 *pmaxvx=maxvx;
2763}
2764/*}}}*/
2765/*FUNCTION Tria::MaxVy{{{1*/
2766void Tria::MaxVy(double* pmaxvy, bool process_units){
2767
2768 /*Get maximum:*/
2769 double maxvy=this->inputs->Max(VyEnum);
2770
2771 /*process units if requested: */
2772 if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum);
2773
2774 /*Assign output pointers:*/
2775 *pmaxvy=maxvy;
2776
2777}
2778/*}}}*/
2779/*FUNCTION Tria::MaxVz{{{1*/
2780void Tria::MaxVz(double* pmaxvz, bool process_units){
2781
2782 /*Get maximum:*/
2783 double maxvz=this->inputs->Max(VzEnum);
2784
2785 /*process units if requested: */
2786 if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum);
2787
2788 /*Assign output pointers:*/
2789 *pmaxvz=maxvz;
2790}
2791/*}}}*/
2792/*FUNCTION Tria::MinVel{{{1*/
2793void Tria::MinVel(double* pminvel, bool process_units){
2794
2795 /*Get minimum:*/
2796 double minvel=this->inputs->Min(VelEnum);
2797
2798 /*process units if requested: */
2799 if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum);
2800
2801 /*Assign output pointers:*/
2802 *pminvel=minvel;
2803}
2804/*}}}*/
2805/*FUNCTION Tria::MinVx{{{1*/
2806void Tria::MinVx(double* pminvx, bool process_units){
2807
2808 /*Get minimum:*/
2809 double minvx=this->inputs->Min(VxEnum);
2810
2811 /*process units if requested: */
2812 if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum);
2813
2814 /*Assign output pointers:*/
2815 *pminvx=minvx;
2816}
2817/*}}}*/
2818/*FUNCTION Tria::MinVy{{{1*/
2819void Tria::MinVy(double* pminvy, bool process_units){
2820
2821 /*Get minimum:*/
2822 double minvy=this->inputs->Min(VyEnum);
2823
2824 /*process units if requested: */
2825 if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum);
2826
2827 /*Assign output pointers:*/
2828 *pminvy=minvy;
2829}
2830/*}}}*/
2831/*FUNCTION Tria::MinVz{{{1*/
2832void Tria::MinVz(double* pminvz, bool process_units){
2833
2834 /*Get minimum:*/
2835 double minvz=this->inputs->Min(VzEnum);
2836
2837 /*process units if requested: */
2838 if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum);
2839
2840 /*Assign output pointers:*/
2841 *pminvz=minvz;
2842}
2843/*}}}*/
2844/*FUNCTION Tria::ElementResponse{{{1*/
2845void Tria::ElementResponse(double* presponse,int response_enum,bool process_units){
2846
2847 switch(response_enum){
2848 case MaterialsRheologyBbarEnum:
2849 *presponse=this->matice->GetBbar();
2850 break;
2851 case MaterialsRheologyZbarEnum:
2852 *presponse=this->matice->GetZbar();
2853 break;
2854 case VelEnum:
2855
2856 /*Get input:*/
2857 double vel;
2858 Input* vel_input;
2859
2860 vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
2861 vel_input->GetInputAverage(&vel);
2862
2863 /*process units if requested: */
2864 if(process_units) vel=UnitConversion(vel,IuToExtEnum,VelEnum);
2865
2866 /*Assign output pointers:*/
2867 *presponse=vel;
2868 default:
2869 _error_("Response type %s not supported yet!",EnumToStringx(response_enum));
2870 }
2871
2872}
2873/*}}}*/
2874#endif
2875
2876#ifdef _HAVE_DIAGNOSTIC_
2877/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/
2878ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyeal(void){
2879
2880 /*compute all stiffness matrices for this element*/
2881 ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyealViscous();
2882 ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction();
2883 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
2884
2885 /*clean-up and return*/
2886 delete Ke1;
2887 delete Ke2;
2888 return Ke;
2889}
2890/*}}}*/
2891/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
2892ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
2893
2894 /*Constants*/
2895 const int numdof=NDOF2*NUMVERTICES;
2896
2897 /*Intermediaries*/
2898 int i,j,ig;
2899 double xyz_list[NUMVERTICES][3];
2900 double viscosity,newviscosity,oldviscosity;
2901 double viscosity_overshoot,thickness,Jdet;
2902 double epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */
2903 double B[3][numdof];
2904 double Bprime[3][numdof];
2905 double D[3][3] = {0.0};
2906 double D_scalar;
2907 GaussTria *gauss = NULL;
2908
2909 /*Initialize Element matrix*/
2910 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
2911
2912 /*Retrieve all inputs and parameters*/
2913 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2914 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
2915 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2916 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2917 Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input);
2918 Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input);
2919 this->parameters->FindParam(&viscosity_overshoot,DiagnosticViscosityOvershootEnum);
2920
2921 /* Start looping on the number of gaussian points: */
2922 gauss=new GaussTria(2);
2923 for (ig=gauss->begin();ig<gauss->end();ig++){
2924
2925 gauss->GaussPoint(ig);
2926
2927 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
2928 GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
2929 GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
2930
2931 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
2932 this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
2933 matice->GetViscosity2d(&viscosity, &epsilon[0]);
2934 matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
2935 thickness_input->GetInputValue(&thickness, gauss);
2936
2937 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
2938 D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
2939 for (i=0;i<3;i++) D[i][i]=D_scalar;
2940
2941 TripleMultiply(&B[0][0],3,numdof,1,
2942 &D[0][0],3,3,0,
2943 &Bprime[0][0],3,numdof,0,
2944 &Ke->values[0],1);
2945 }
2946
2947 /*Transform Coordinate System*/
2948 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
2949
2950 /*Clean up and return*/
2951 delete gauss;
2952 return Ke;
2953}
2954/*}}}*/
2955/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
2956ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
2957
2958 /*Constants*/
2959 const int numdof=NDOF2*NUMVERTICES;
2960
2961 /*Intermediaries*/
2962 int i,j,ig;
2963 int analysis_type;
2964 double MAXSLOPE = .06; // 6 %
2965 double MOUNTAINKEXPONENT = 10;
2966 double slope_magnitude,alpha2;
2967 double Jdet;
2968 double L[2][numdof];
2969 double DL[2][2] = {{ 0,0 },{0,0}};
2970 double DL_scalar;
2971 double slope[2] = {0.0,0.0};
2972 double xyz_list[NUMVERTICES][3];
2973 Friction *friction = NULL;
2974 GaussTria *gauss = NULL;
2975
2976 /*Initialize Element matrix and return if necessary*/
2977 if(IsFloating()) return NULL;
2978 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
2979
2980 /*Retrieve all inputs and parameters*/
2981 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2982 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
2983 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2984 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2985 Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
2986 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
2987
2988 /*build friction object, used later on: */
2989 friction=new Friction("2d",inputs,matpar,analysis_type);
2990
2991 /* Start looping on the number of gaussian points: */
2992 gauss=new GaussTria(2);
2993 for (ig=gauss->begin();ig<gauss->end();ig++){
2994
2995 gauss->GaussPoint(ig);
2996
2997 // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case,
2998 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
2999 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
3000 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
3001 if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT);
3002 else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
3003
3004 GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
3005 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3006 DL_scalar=alpha2*gauss->weight*Jdet;
3007 for (i=0;i<2;i++) DL[i][i]=DL_scalar;
3008
3009 TripleMultiply( &L[0][0],2,numdof,1,
3010 &DL[0][0],2,2,0,
3011 &L[0][0],2,numdof,0,
3012 &Ke->values[0],1);
3013 }
3014
3015 /*Transform Coordinate System*/
3016 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
3017
3018 /*Clean up and return*/
3019 delete gauss;
3020 delete friction;
3021 return Ke;
3022}
3023/*}}}*/
3024/*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/
3025ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){
3026
3027 /*Intermediaries*/
3028 const int numdof=NUMVERTICES*NDOF2;
3029 int i,connectivity;
3030
3031 /*Initialize Element matrix*/
3032 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
3033
3034 /*Create Element matrix*/
3035 for(i=0;i<NUMVERTICES;i++){
3036 connectivity=nodes[i]->GetConnectivity();
3037 Ke->values[(2*i)*numdof +(2*i) ]=1/(double)connectivity;
3038 Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity;
3039 }
3040
3041 /*Clean up and return*/
3042 return Ke;
3043}
3044/*}}}*/
3045/*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/
3046ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
3047
3048 /*Constants*/
3049 const int numdof=NDOF2*NUMVERTICES;
3050
3051 /*Intermediaries */
3052 int i,j,ig;
3053 double driving_stress_baseline,thickness;
3054 double Jdet;
3055 double xyz_list[NUMVERTICES][3];
3056 double slope[2];
3057 double basis[3];
3058 double pe_g_gaussian[numdof];
3059 GaussTria* gauss=NULL;
3060
3061 /*Initialize Element vector*/
3062 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
3063
3064 /*Retrieve all inputs and parameters*/
3065 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3066 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3067 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
3068 Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
3069
3070 /* Start looping on the number of gaussian points: */
3071 gauss=new GaussTria(2);
3072 for(ig=gauss->begin();ig<gauss->end();ig++){
3073
3074 gauss->GaussPoint(ig);
3075
3076 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3077 GetNodalFunctions(basis, gauss);
3078
3079 thickness_input->GetInputValue(&thickness,gauss);
3080 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
3081 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
3082
3083 /*Build pe_g_gaussian vector: */
3084 for (i=0;i<NUMVERTICES;i++){
3085 for (j=0;j<NDOF2;j++){
3086 pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
3087 }
3088 }
3089 }
3090
3091 /*Transform coordinate system*/
3092 TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);
3093
3094 /*Clean up and return*/
3095 delete gauss;
3096 return pe;
3097}
3098/*}}}*/
3099/*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
3100ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
3101
3102 /*Intermediaries */
3103 int i,connectivity;
3104 double constant_part,ub,vb;
3105 double rho_ice,gravity,n,B;
3106 double slope2,thickness;
3107 double slope[2];
3108 GaussTria* gauss=NULL;
3109
3110 /*Initialize Element vector*/
3111 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
3112
3113 /*Retrieve all inputs and parameters*/
3114 rho_ice=matpar->GetRhoIce();
3115 gravity=matpar->GetG();
3116 n=matice->GetN();
3117 B=matice->GetBbar();
3118 Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
3119 Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
3120 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3121
3122 /*Spawn 3 sing elements: */
3123 gauss=new GaussTria();
3124 for(i=0;i<NUMVERTICES;i++){
3125
3126 gauss->GaussVertex(i);
3127
3128 connectivity=nodes[i]->GetConnectivity();
3129
3130 thickness_input->GetInputValue(&thickness,gauss);
3131 slopex_input->GetInputValue(&slope[0],gauss);
3132 slopey_input->GetInputValue(&slope[1],gauss);
3133 slope2=pow(slope[0],2)+pow(slope[1],2);
3134
3135 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
3136
3137 ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
3138 vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1];
3139
3140 pe->values[2*i] =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity;
3141 pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity;
3142 }
3143
3144 /*Clean up and return*/
3145 delete gauss;
3146 return pe;
3147}
3148/*}}}*/
3149/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{1*/
3150ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
3151
3152 /*Constants*/
3153 const int numdof=NDOF2*NUMVERTICES;
3154
3155 /*Intermediaries */
3156 int i,j,ig;
3157 double xyz_list[NUMVERTICES][3];
3158 double Jdet,thickness;
3159 double eps1dotdphii,eps1dotdphij;
3160 double eps2dotdphii,eps2dotdphij;
3161 double mu_prime;
3162 double epsilon[3];/* epsilon=[exx,eyy,exy];*/
3163 double eps1[2],eps2[2];
3164 double phi[NUMVERTICES];
3165 double dphi[2][NUMVERTICES];
3166 GaussTria *gauss=NULL;
3167
3168 /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
3169 ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
3170
3171 /*Retrieve all inputs and parameters*/
3172 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3173 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3174 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3175 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3176
3177 /* Start looping on the number of gaussian points: */
3178 gauss=new GaussTria(2);
3179 for (ig=gauss->begin();ig<gauss->end();ig++){
3180
3181 gauss->GaussPoint(ig);
3182
3183 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3184 GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
3185
3186 thickness_input->GetInputValue(&thickness, gauss);
3187 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
3188 matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
3189 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
3190 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
3191
3192 for(i=0;i<3;i++){
3193 for(j=0;j<3;j++){
3194 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
3195 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
3196 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
3197 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
3198
3199 Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
3200 Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
3201 Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
3202 Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
3203 }
3204 }
3205 }
3206
3207 /*Transform Coordinate System*/
3208 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
3209
3210 /*Clean up and return*/
3211 delete gauss;
3212 return Ke;
3213}
3214/*}}}*/
3215/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
3216void Tria::GetSolutionFromInputsDiagnosticHoriz(Vector* solution){
3217
3218 const int numdof=NDOF2*NUMVERTICES;
3219
3220 int i;
3221 int* doflist=NULL;
3222 double vx,vy;
3223 double values[numdof];
3224 GaussTria* gauss=NULL;
3225
3226 /*Get dof list: */
3227 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3228
3229 /*Get inputs*/
3230 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3231 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3232
3233 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3234 /*P1 element only for now*/
3235 gauss=new GaussTria();
3236 for(i=0;i<NUMVERTICES;i++){
3237
3238 gauss->GaussVertex(i);
3239
3240 /*Recover vx and vy*/
3241 vx_input->GetInputValue(&vx,gauss);
3242 vy_input->GetInputValue(&vy,gauss);
3243 values[i*NDOF2+0]=vx;
3244 values[i*NDOF2+1]=vy;
3245 }
3246
3247 solution->SetValues(numdof,doflist,values,INS_VAL);
3248
3249 /*Free ressources:*/
3250 delete gauss;
3251 xfree((void**)&doflist);
3252}
3253/*}}}*/
3254/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
3255void Tria::GetSolutionFromInputsDiagnosticHutter(Vector* solution){
3256
3257 const int numdof=NDOF2*NUMVERTICES;
3258
3259 int i;
3260 double vx,vy;
3261 double values[numdof];
3262 int *doflist = NULL;
3263 GaussTria *gauss = NULL;
3264
3265 /*Get dof list: */
3266 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3267
3268 /*Get inputs*/
3269 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3270 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3271
3272 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3273 /*P1 element only for now*/
3274 gauss=new GaussTria();
3275 for(i=0;i<NUMVERTICES;i++){
3276
3277 gauss->GaussVertex(i);
3278
3279 /*Recover vx and vy*/
3280 vx_input->GetInputValue(&vx,gauss);
3281 vy_input->GetInputValue(&vy,gauss);
3282 values[i*NDOF2+0]=vx;
3283 values[i*NDOF2+1]=vy;
3284 }
3285
3286 solution->SetValues(numdof,doflist,values,INS_VAL);
3287
3288 /*Free ressources:*/
3289 delete gauss;
3290 xfree((void**)&doflist);
3291}
3292/*}}}*/
3293/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
3294void Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
3295
3296 const int numdof=NDOF2*NUMVERTICES;
3297
3298 int i;
3299 int* doflist=NULL;
3300 double rho_ice,g;
3301 double values[numdof];
3302 double vx[NUMVERTICES];
3303 double vy[NUMVERTICES];
3304 double vz[NUMVERTICES];
3305 double vel[NUMVERTICES];
3306 double pressure[NUMVERTICES];
3307 double thickness[NUMVERTICES];
3308
3309 /*Get dof list: */
3310 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3311
3312 /*Use the dof list to index into the solution vector: */
3313 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
3314
3315 /*Transform solution in Cartesian Space*/
3316 TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
3317
3318 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3319 for(i=0;i<NUMVERTICES;i++){
3320 vx[i]=values[i*NDOF2+0];
3321 vy[i]=values[i*NDOF2+1];
3322
3323 /*Check solution*/
3324 if(isnan(vx[i])) _error_("NaN found in solution vector");
3325 if(isnan(vy[i])) _error_("NaN found in solution vector");
3326 }
3327
3328 /*Get Vz and compute vel*/
3329 GetInputListOnVertices(&vz[0],VzEnum,0);
3330 for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
3331
3332 /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
3333 *so the pressure is just the pressure at the bedrock: */
3334 rho_ice=matpar->GetRhoIce();
3335 g=matpar->GetG();
3336 GetInputListOnVertices(&thickness[0],ThicknessEnum);
3337 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
3338
3339 /*Now, we have to move the previous Vx and Vy inputs to old
3340 * status, otherwise, we'll wipe them off: */
3341 this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
3342 this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
3343 this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
3344
3345 /*Add vx and vy as inputs to the tria element: */
3346 this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
3347 this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
3348 this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
3349 this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));
3350
3351 /*Free ressources:*/
3352 xfree((void**)&doflist);
3353
3354}
3355/*}}}*/
3356/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{1*/
3357void Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){
3358
3359 const int numdof=NDOF2*NUMVERTICES;
3360
3361 int i;
3362 int* doflist=NULL;
3363 double rho_ice,g;
3364 double values[numdof];
3365 double vx[NUMVERTICES];
3366 double vy[NUMVERTICES];
3367 double vz[NUMVERTICES];
3368 double vel[NUMVERTICES];
3369 double pressure[NUMVERTICES];
3370 double thickness[NUMVERTICES];
3371
3372 /*Get dof list: */
3373 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3374
3375 /*Use the dof list to index into the solution vector: */
3376 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
3377
3378 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3379 for(i=0;i<NUMVERTICES;i++){
3380 vx[i]=values[i*NDOF2+0];
3381 vy[i]=values[i*NDOF2+1];
3382
3383 /*Check solution*/
3384 if(isnan(vx[i])) _error_("NaN found in solution vector");
3385 if(isnan(vy[i])) _error_("NaN found in solution vector");
3386 }
3387
3388 /*Now Compute vel*/
3389 GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
3390 for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
3391
3392 /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
3393 *so the pressure is just the pressure at the bedrock: */
3394 rho_ice=matpar->GetRhoIce();
3395 g=matpar->GetG();
3396 GetInputListOnVertices(&thickness[0],ThicknessEnum);
3397 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
3398
3399 /*Now, we have to move the previous Vx and Vy inputs to old
3400 * status, otherwise, we'll wipe them off: */
3401 this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
3402 this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
3403 this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
3404
3405 /*Add vx and vy as inputs to the tria element: */
3406 this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
3407 this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
3408 this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
3409 this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));
3410
3411 /*Free ressources:*/
3412 xfree((void**)&doflist);
3413}
3414/*}}}*/
3415#endif
3416
3417#ifdef _HAVE_CONTROL_
3418/*FUNCTION Tria::InputControlUpdate{{{1*/
3419void Tria::InputControlUpdate(double scalar,bool save_parameter){
3420
3421 /*Intermediary*/
3422 int num_controls;
3423 int* control_type=NULL;
3424 Input* input=NULL;
3425
3426 /*retrieve some parameters: */
3427 this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
3428 this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
3429
3430 for(int i=0;i<num_controls;i++){
3431
3432 if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){
3433 input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
3434 }
3435 else{
3436 input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
3437 }
3438
3439 if (input->ObjectEnum()!=ControlInputEnum){
3440 _error_("input %s is not a ControlInput",EnumToStringx(control_type[i]));
3441 }
3442
3443 ((ControlInput*)input)->UpdateValue(scalar);
3444 ((ControlInput*)input)->Constrain();
3445 if (save_parameter) ((ControlInput*)input)->SaveValue();
3446
3447 }
3448
3449 /*Clean up and return*/
3450 xfree((void**)&control_type);
3451}
3452/*}}}*/
3453/*FUNCTION Tria::ControlInputGetGradient{{{1*/
3454void Tria::ControlInputGetGradient(Vector* gradient,int enum_type,int control_index){
3455
3456 int doflist1[NUMVERTICES];
3457 Input* input=NULL;
3458
3459 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
3460 input=(Input*)matice->inputs->GetInput(enum_type);
3461 }
3462 else{
3463 input=inputs->GetInput(enum_type);
3464 }
3465 if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
3466 if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
3467
3468 GradientIndexing(&doflist1[0],control_index);
3469 ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]);
3470
3471}/*}}}*/
3472/*FUNCTION Tria::ControlInputScaleGradient{{{1*/
3473void Tria::ControlInputScaleGradient(int enum_type,double scale){
3474
3475 Input* input=NULL;
3476
3477 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
3478 input=(Input*)matice->inputs->GetInput(enum_type);
3479 }
3480 else{
3481 input=inputs->GetInput(enum_type);
3482 }
3483 if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
3484 if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
3485
3486 ((ControlInput*)input)->ScaleGradient(scale);
3487}/*}}}*/
3488/*FUNCTION Tria::ControlInputSetGradient{{{1*/
3489void Tria::ControlInputSetGradient(double* gradient,int enum_type,int control_index){
3490
3491 int doflist1[NUMVERTICES];
3492 double grad_list[NUMVERTICES];
3493 Input* grad_input=NULL;
3494 Input* input=NULL;
3495
3496 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
3497 input=(Input*)matice->inputs->GetInput(enum_type);
3498 }
3499 else{
3500 input=inputs->GetInput(enum_type);
3501 }
3502 if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
3503 if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
3504
3505 GradientIndexing(&doflist1[0],control_index);
3506 for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[doflist1[i]];
3507 grad_input=new TriaP1Input(GradientEnum,grad_list);
3508
3509 ((ControlInput*)input)->SetGradient(grad_input);
3510
3511}/*}}}*/
3512/*FUNCTION Tria::Gradj {{{1*/
3513void Tria::Gradj(Vector* gradient,int control_type,int control_index){
3514 /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
3515
3516 /*If on water, grad = 0: */
3517 if(IsOnWater()) return;
3518
3519 /*First deal with ∂/∂alpha(KU-F)*/
3520 switch(control_type){
3521 case FrictionCoefficientEnum:
3522 GradjDragMacAyeal(gradient,control_index);
3523 break;
3524 case MaterialsRheologyBbarEnum:
3525 GradjBMacAyeal(gradient,control_index);
3526 break;
3527 case MaterialsRheologyZbarEnum:
3528 GradjZMacAyeal(gradient,control_index);
3529 break;
3530 case BalancethicknessThickeningRateEnum:
3531 GradjDhDtBalancedthickness(gradient,control_index);
3532 break;
3533 case VxEnum:
3534 GradjVxBalancedthickness(gradient,control_index);
3535 break;
3536 case VyEnum:
3537 GradjVyBalancedthickness(gradient,control_index);
3538 break;
3539 default:
3540 _error_("%s%i","control type not supported yet: ",control_type);
3541 }
3542
3543 /*Now deal with ∂J/∂alpha*/
3544 int *responses = NULL;
3545 int num_responses,resp;
3546 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
3547 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
3548
3549 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
3550 //FIXME: the control type should be checked somewhere (with respect to what variable are we taking the gradient!)
3551
3552 case ThicknessAbsMisfitEnum:
3553 case ThicknessAbsGradientEnum:
3554 case SurfaceAbsVelMisfitEnum:
3555 case SurfaceRelVelMisfitEnum:
3556 case SurfaceLogVelMisfitEnum:
3557 case SurfaceLogVxVyMisfitEnum:
3558 case SurfaceAverageVelMisfitEnum:
3559 /*Nothing, J does not depends on the parameter being inverted for*/
3560 break;
3561 case DragCoefficientAbsGradientEnum:
3562 GradjDragGradient(gradient,resp,control_index);
3563 break;
3564 case RheologyBbarAbsGradientEnum:
3565 GradjBGradient(gradient,resp,control_index);
3566 break;
3567 default:
3568 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
3569 }
3570
3571 xfree((void**)&responses);
3572}
3573/*}}}*/
3574/*FUNCTION Tria::GradjBGradient{{{1*/
3575void Tria::GradjBGradient(Vector* gradient,int weight_index,int control_index){
3576
3577 int i,ig;
3578 int doflist1[NUMVERTICES];
3579 double Jdet,weight;
3580 double xyz_list[NUMVERTICES][3];
3581 double dbasis[NDOF2][NUMVERTICES];
3582 double dk[NDOF2];
3583 double grade_g[NUMVERTICES]={0.0};
3584 GaussTria *gauss=NULL;
3585
3586 /*Retrieve all inputs we will be needing: */
3587 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3588 GradientIndexing(&doflist1[0],control_index);
3589 Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
3590 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
3591
3592 /* Start looping on the number of gaussian points: */
3593 gauss=new GaussTria(2);
3594 for (ig=gauss->begin();ig<gauss->end();ig++){
3595
3596 gauss->GaussPoint(ig);
3597
3598 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3599 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
3600 weights_input->GetInputValue(&weight,gauss,weight_index);
3601
3602 /*Build alpha_complement_list: */
3603 rheologyb_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3604
3605 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3606 for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
3607 }
3608 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3609
3610 /*Clean up and return*/
3611 delete gauss;
3612}
3613/*}}}*/
3614/*FUNCTION Tria::GradjZGradient{{{1*/
3615void Tria::GradjZGradient(Vector* gradient,int weight_index,int control_index){
3616
3617 int i,ig;
3618 int doflist1[NUMVERTICES];
3619 double Jdet,weight;
3620 double xyz_list[NUMVERTICES][3];
3621 double dbasis[NDOF2][NUMVERTICES];
3622 double dk[NDOF2];
3623 double grade_g[NUMVERTICES]={0.0};
3624 GaussTria *gauss=NULL;
3625
3626 /*Retrieve all inputs we will be needing: */
3627 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3628 GradientIndexing(&doflist1[0],control_index);
3629 Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
3630 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
3631
3632 /* Start looping on the number of gaussian points: */
3633 gauss=new GaussTria(2);
3634 for (ig=gauss->begin();ig<gauss->end();ig++){
3635
3636 gauss->GaussPoint(ig);
3637
3638 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3639 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
3640 weights_input->GetInputValue(&weight,gauss,weight_index);
3641
3642 /*Build alpha_complement_list: */
3643 rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3644
3645 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3646 for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
3647 }
3648 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3649
3650 /*Clean up and return*/
3651 delete gauss;
3652}
3653/*}}}*/
3654/*FUNCTION Tria::GradjBMacAyeal{{{1*/
3655void Tria::GradjBMacAyeal(Vector* gradient,int control_index){
3656
3657 /*Intermediaries*/
3658 int i,ig;
3659 int doflist[NUMVERTICES];
3660 double vx,vy,lambda,mu,thickness,Jdet;
3661 double viscosity_complement;
3662 double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2];
3663 double xyz_list[NUMVERTICES][3];
3664 double basis[3],epsilon[3];
3665 double grad[NUMVERTICES]={0.0};
3666 GaussTria *gauss = NULL;
3667
3668 /* Get node coordinates and dof list: */
3669 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3670 GradientIndexing(&doflist[0],control_index);
3671
3672 /*Retrieve all inputs*/
3673 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3674 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3675 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3676 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input);
3677 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input);
3678 Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
3679
3680 /* Start looping on the number of gaussian points: */
3681 gauss=new GaussTria(4);
3682 for (ig=gauss->begin();ig<gauss->end();ig++){
3683
3684 gauss->GaussPoint(ig);
3685
3686 thickness_input->GetInputValue(&thickness,gauss);
3687 rheologyb_input->GetInputDerivativeValue(&dB[0],&xyz_list[0][0],gauss);
3688 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
3689 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
3690 adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
3691 adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
3692
3693 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
3694 matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);
3695
3696 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3697 GetNodalFunctions(basis,gauss);
3698
3699 /*standard gradient dJ/dki*/
3700 for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
3701 (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
3702 )*Jdet*gauss->weight*basis[i];
3703 }
3704
3705 gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
3706
3707 /*clean-up*/
3708 delete gauss;
3709}
3710/*}}}*/
3711/*FUNCTION Tria::GradjZMacAyeal{{{1*/
3712void Tria::GradjZMacAyeal(Vector* gradient,int control_index){
3713
3714 /*Intermediaries*/
3715 int i,ig;
3716 int doflist[NUMVERTICES];
3717 double vx,vy,lambda,mu,thickness,Jdet;
3718 double viscosity_complement;
3719 double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2];
3720 double xyz_list[NUMVERTICES][3];
3721 double basis[3],epsilon[3];
3722 double grad[NUMVERTICES]={0.0};
3723 GaussTria *gauss = NULL;
3724
3725 /* Get node coordinates and dof list: */
3726 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3727 GradientIndexing(&doflist[0],control_index);
3728
3729 /*Retrieve all inputs*/
3730 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3731 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3732 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3733 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input);
3734 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input);
3735 Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
3736
3737 /* Start looping on the number of gaussian points: */
3738 gauss=new GaussTria(4);
3739 for (ig=gauss->begin();ig<gauss->end();ig++){
3740
3741 gauss->GaussPoint(ig);
3742
3743 thickness_input->GetInputValue(&thickness,gauss);
3744 rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss);
3745 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
3746 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
3747 adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
3748 adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
3749
3750 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
3751 matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]);
3752
3753 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3754 GetNodalFunctions(basis,gauss);
3755
3756 /*standard gradient dJ/dki*/
3757 for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
3758 (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
3759 )*Jdet*gauss->weight*basis[i];
3760 }
3761
3762 gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
3763
3764 /*clean-up*/
3765 delete gauss;
3766}
3767/*}}}*/
3768/*FUNCTION Tria::GradjDragMacAyeal {{{1*/
3769void Tria::GradjDragMacAyeal(Vector* gradient,int control_index){
3770
3771 int i,ig;
3772 int analysis_type;
3773 int doflist1[NUMVERTICES];
3774 int connectivity[NUMVERTICES];
3775 double vx,vy,lambda,mu,alpha_complement,Jdet;
3776 double bed,thickness,Neff,drag;
3777 double xyz_list[NUMVERTICES][3];
3778 double dk[NDOF2];
3779 double grade_g[NUMVERTICES]={0.0};
3780 double grade_g_gaussian[NUMVERTICES];
3781 double basis[3];
3782 double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
3783 Friction* friction=NULL;
3784 GaussTria *gauss=NULL;
3785
3786 if(IsFloating())return;
3787
3788 /*retrive parameters: */
3789 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
3790 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3791 GradientIndexing(&doflist1[0],control_index);
3792 this->GetConnectivityList(&connectivity[0]);
3793
3794 /*Build friction element, needed later: */
3795 friction=new Friction("2d",inputs,matpar,analysis_type);
3796
3797 /*Retrieve all inputs we will be needing: */
3798 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input);
3799 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input);
3800 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3801 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3802 Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
3803
3804 /* Start looping on the number of gaussian points: */
3805 gauss=new GaussTria(4);
3806 for (ig=gauss->begin();ig<gauss->end();ig++){
3807
3808 gauss->GaussPoint(ig);
3809
3810 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3811 GetNodalFunctions(basis, gauss);
3812
3813 /*Build alpha_complement_list: */
3814 friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
3815
3816 dragcoefficient_input->GetInputValue(&drag, gauss);
3817 adjointx_input->GetInputValue(&lambda, gauss);
3818 adjointy_input->GetInputValue(&mu, gauss);
3819 vx_input->GetInputValue(&vx,gauss);
3820 vy_input->GetInputValue(&vy,gauss);
3821 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3822
3823 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3824 for (i=0;i<NUMVERTICES;i++){
3825 grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
3826 }
3827
3828 /*Add gradje_g_gaussian vector to gradje_g: */
3829 for(i=0;i<NUMVERTICES;i++){
3830 _assert_(!isnan(grade_g[i]));
3831 grade_g[i]+=grade_g_gaussian[i];
3832 }
3833 }
3834 /*Analytical gradient*/
3835 //delete gauss;
3836 //gauss=new GaussTria();
3837 //for (int iv=0;iv<NUMVERTICES;iv++){
3838 // gauss->GaussVertex(iv);
3839 // friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
3840 // dragcoefficient_input->GetInputValue(&drag, gauss);
3841 // adjointx_input->GetInputValue(&lambda, gauss);
3842 // adjointy_input->GetInputValue(&mu, gauss);
3843 // vx_input->GetInputValue(&vx,gauss);
3844 // vy_input->GetInputValue(&vy,gauss);
3845 // grade_g[iv] = -2*1.e+7*drag*alpha_complement*(lambda*vx+mu*vy)/((double)connectivity[iv]);
3846 //}
3847 /*End Analytical gradient*/
3848
3849 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3850
3851 /*Clean up and return*/
3852 delete gauss;
3853 delete friction;
3854}
3855/*}}}*/
3856/*FUNCTION Tria::GradjDragGradient{{{1*/
3857void Tria::GradjDragGradient(Vector* gradient, int weight_index,int control_index){
3858
3859 int i,ig;
3860 int doflist1[NUMVERTICES];
3861 double Jdet,weight;
3862 double xyz_list[NUMVERTICES][3];
3863 double dbasis[NDOF2][NUMVERTICES];
3864 double dk[NDOF2];
3865 double grade_g[NUMVERTICES]={0.0};
3866 GaussTria *gauss=NULL;
3867
3868 /*Retrieve all inputs we will be needing: */
3869 if(IsFloating())return;
3870 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3871 GradientIndexing(&doflist1[0],control_index);
3872 Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
3873 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
3874
3875 /* Start looping on the number of gaussian points: */
3876 gauss=new GaussTria(2);
3877 for (ig=gauss->begin();ig<gauss->end();ig++){
3878
3879 gauss->GaussPoint(ig);
3880
3881 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3882 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
3883 weights_input->GetInputValue(&weight,gauss,weight_index);
3884
3885 /*Build alpha_complement_list: */
3886 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3887
3888 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3889 for (i=0;i<NUMVERTICES;i++){
3890 grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
3891 _assert_(!isnan(grade_g[i]));
3892 }
3893 }
3894 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3895
3896 /*Clean up and return*/
3897 delete gauss;
3898}
3899/*}}}*/
3900/*FUNCTION Tria::GradjDhDtBalancedthickness{{{1*/
3901void Tria::GradjDhDtBalancedthickness(Vector* gradient,int control_index){
3902
3903 /*Intermediaries*/
3904 int doflist1[NUMVERTICES];
3905 double lambda[NUMVERTICES];
3906 double gradient_g[NUMVERTICES];
3907
3908 /*Compute Gradient*/
3909 GradientIndexing(&doflist1[0],control_index);
3910 GetInputListOnVertices(&lambda[0],AdjointEnum);
3911 for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i];
3912
3913 gradient->SetValues(NUMVERTICES,doflist1,gradient_g,INS_VAL);
3914}
3915/*}}}*/
3916/*FUNCTION Tria::GradjVxBalancedthickness{{{1*/
3917void Tria::GradjVxBalancedthickness(Vector* gradient,int control_index){
3918
3919 /*Intermediaries*/
3920 int i,ig;
3921 int doflist1[NUMVERTICES];
3922 double thickness,Jdet;
3923 double basis[3];
3924 double Dlambda[2],dp[2];
3925 double xyz_list[NUMVERTICES][3];
3926 double grade_g[NUMVERTICES] = {0.0};
3927 GaussTria *gauss = NULL;
3928
3929 /* Get node coordinates and dof list: */
3930 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3931 GradientIndexing(&doflist1[0],control_index);
3932
3933 /*Retrieve all inputs we will be needing: */
3934 Input* adjoint_input=inputs->GetInput(AdjointEnum); _assert_(adjoint_input);
3935 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3936
3937 /* Start looping on the number of gaussian points: */
3938 gauss=new GaussTria(2);
3939 for (ig=gauss->begin();ig<gauss->end();ig++){
3940
3941 gauss->GaussPoint(ig);
3942
3943 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3944 GetNodalFunctions(basis, gauss);
3945
3946 adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
3947 thickness_input->GetInputValue(&thickness, gauss);
3948 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
3949
3950 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
3951 }
3952
3953 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3954
3955 /*Clean up and return*/
3956 delete gauss;
3957}
3958/*}}}*/
3959/*FUNCTION Tria::GradjVyBalancedthickness{{{1*/
3960void Tria::GradjVyBalancedthickness(Vector* gradient,int control_index){
3961
3962 /*Intermediaries*/
3963 int i,ig;
3964 int doflist1[NUMVERTICES];
3965 double thickness,Jdet;
3966 double basis[3];
3967 double Dlambda[2],dp[2];
3968 double xyz_list[NUMVERTICES][3];
3969 double grade_g[NUMVERTICES] = {0.0};
3970 GaussTria *gauss = NULL;
3971
3972 /* Get node coordinates and dof list: */
3973 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3974 GradientIndexing(&doflist1[0],control_index);
3975
3976 /*Retrieve all inputs we will be needing: */
3977 Input* adjoint_input=inputs->GetInput(AdjointEnum); _assert_(adjoint_input);
3978 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3979
3980 /* Start looping on the number of gaussian points: */
3981 gauss=new GaussTria(2);
3982 for (ig=gauss->begin();ig<gauss->end();ig++){
3983
3984 gauss->GaussPoint(ig);
3985
3986 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3987 GetNodalFunctions(basis, gauss);
3988
3989 adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
3990 thickness_input->GetInputValue(&thickness, gauss);
3991 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
3992
3993 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
3994 }
3995 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3996
3997 /*Clean up and return*/
3998 delete gauss;
3999}
4000/*}}}*/
4001/*FUNCTION Tria::GradientIndexing{{{1*/
4002void Tria::GradientIndexing(int* indexing,int control_index){
4003
4004 /*Get some parameters*/
4005 int num_controls;
4006 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
4007
4008 /*get gradient indices*/
4009 for(int i=0;i<NUMVERTICES;i++){
4010 indexing[i]=num_controls*this->nodes[i]->GetVertexDof() + control_index;
4011 }
4012
4013}
4014/*}}}*/
4015/*FUNCTION Tria::RheologyBbarAbsGradient{{{1*/
4016double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){
4017
4018 /* Intermediaries */
4019 int ig;
4020 double Jelem = 0;
4021 double weight;
4022 double Jdet;
4023 double xyz_list[NUMVERTICES][3];
4024 double dp[NDOF2];
4025 GaussTria *gauss = NULL;
4026
4027 /*retrieve parameters and inputs*/
4028
4029 /*If on water, return 0: */
4030 if(IsOnWater()) return 0;
4031
4032 /*Retrieve all inputs we will be needing: */
4033 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4034 Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4035 Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
4036
4037 /* Start looping on the number of gaussian points: */
4038 gauss=new GaussTria(2);
4039 for (ig=gauss->begin();ig<gauss->end();ig++){
4040
4041 gauss->GaussPoint(ig);
4042
4043 /* Get Jacobian determinant: */
4044 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4045
4046 /*Get all parameters at gaussian point*/
4047 weights_input->GetInputValue(&weight,gauss,weight_index);
4048 rheologyb_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
4049
4050 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
4051 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
4052 }
4053
4054 /*Clean up and return*/
4055 delete gauss;
4056 return Jelem;
4057}
4058/*}}}*/
4059/*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/
4060double Tria::SurfaceAverageVelMisfit(bool process_units,int weight_index){
4061
4062 const int numdof=2*NUMVERTICES;
4063
4064 int i,ig;
4065 double Jelem=0,S,Jdet;
4066 double misfit;
4067 double vx,vy,vxobs,vyobs,weight;
4068 double xyz_list[NUMVERTICES][3];
4069 GaussTria *gauss=NULL;
4070
4071 /*If on water, return 0: */
4072 if(IsOnWater())return 0;
4073
4074 /* Get node coordinates and dof list: */
4075 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4076
4077 /*Retrieve all inputs we will be needing: */
4078 inputs->GetInputValue(&S,SurfaceAreaEnum);
4079 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4080 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4081 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4082 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4083 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4084
4085 /* Start looping on the number of gaussian points: */
4086 gauss=new GaussTria(3);
4087 for (ig=gauss->begin();ig<gauss->end();ig++){
4088
4089 gauss->GaussPoint(ig);
4090
4091 /* Get Jacobian determinant: */
4092 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4093
4094 /*Get all parameters at gaussian point*/
4095 weights_input->GetInputValue(&weight,gauss,weight_index);
4096 vx_input->GetInputValue(&vx,gauss);
4097 vy_input->GetInputValue(&vy,gauss);
4098 vxobs_input->GetInputValue(&vxobs,gauss);
4099 vyobs_input->GetInputValue(&vyobs,gauss);
4100
4101 /*Compute SurfaceAverageVelMisfitEnum:
4102 *
4103 * 1 2 2
4104 * J = --- sqrt( (u - u ) + (v - v ) )
4105 * S obs obs
4106 */
4107 misfit=1/S*pow( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) ,0.5);
4108
4109 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
4110
4111 /*Add to cost function*/
4112 Jelem+=misfit*weight*Jdet*gauss->weight;
4113 }
4114
4115 /*clean-up and Return: */
4116 delete gauss;
4117 return Jelem;
4118}
4119/*}}}*/
4120/*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/
4121double Tria::SurfaceLogVelMisfit(bool process_units,int weight_index){
4122
4123 const int numdof=NDOF2*NUMVERTICES;
4124
4125 int i,ig;
4126 double Jelem=0;
4127 double misfit,Jdet;
4128 double epsvel=2.220446049250313e-16;
4129 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4130 double velocity_mag,obs_velocity_mag;
4131 double xyz_list[NUMVERTICES][3];
4132 double vx,vy,vxobs,vyobs,weight;
4133 GaussTria *gauss=NULL;
4134
4135 /*If on water, return 0: */
4136 if(IsOnWater())return 0;
4137
4138 /* Get node coordinates and dof list: */
4139 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4140
4141 /*Retrieve all inputs we will be needing: */
4142 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4143 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4144 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4145 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4146 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4147
4148 /* Start looping on the number of gaussian points: */
4149 gauss=new GaussTria(4);
4150 for (ig=gauss->begin();ig<gauss->end();ig++){
4151
4152 gauss->GaussPoint(ig);
4153
4154 /* Get Jacobian determinant: */
4155 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4156
4157 /*Get all parameters at gaussian point*/
4158 weights_input->GetInputValue(&weight,gauss,weight_index);
4159 vx_input->GetInputValue(&vx,gauss);
4160 vy_input->GetInputValue(&vy,gauss);
4161 vxobs_input->GetInputValue(&vxobs,gauss);
4162 vyobs_input->GetInputValue(&vyobs,gauss);
4163
4164 /*Compute SurfaceLogVelMisfit:
4165 * [ vel + eps ] 2
4166 * J = 4 \bar{v}^2 | log ( ----------- ) |
4167 * [ vel + eps ]
4168 * obs
4169 */
4170 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel;
4171 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
4172 misfit=4*pow(meanvel,2.)*pow(log(velocity_mag/obs_velocity_mag),2.);
4173
4174 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVelMisfitEnum);
4175
4176 /*Add to cost function*/
4177 Jelem+=misfit*weight*Jdet*gauss->weight;
4178 }
4179
4180 /*clean-up and Return: */
4181 delete gauss;
4182 return Jelem;
4183}
4184/*}}}*/
4185/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
4186double Tria::SurfaceLogVxVyMisfit(bool process_units,int weight_index){
4187
4188 const int numdof=NDOF2*NUMVERTICES;
4189
4190 int i,ig;
4191 int fit=-1;
4192 double Jelem=0, S=0;
4193 double epsvel=2.220446049250313e-16;
4194 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4195 double misfit, Jdet;
4196 double vx,vy,vxobs,vyobs,weight;
4197 double xyz_list[NUMVERTICES][3];
4198 GaussTria *gauss=NULL;
4199
4200 /*If on water, return 0: */
4201 if(IsOnWater())return 0;
4202
4203 /* Get node coordinates and dof list: */
4204 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4205
4206 /*Retrieve all inputs we will be needing: */
4207 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4208 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4209 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4210 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4211 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4212
4213 /* Start looping on the number of gaussian points: */
4214 gauss=new GaussTria(4);
4215 for (ig=gauss->begin();ig<gauss->end();ig++){
4216
4217 gauss->GaussPoint(ig);
4218
4219 /* Get Jacobian determinant: */
4220 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4221
4222 /*Get all parameters at gaussian point*/
4223 weights_input->GetInputValue(&weight,gauss,weight_index);
4224 vx_input->GetInputValue(&vx,gauss);
4225 vy_input->GetInputValue(&vy,gauss);
4226 vxobs_input->GetInputValue(&vxobs,gauss);
4227 vyobs_input->GetInputValue(&vyobs,gauss);
4228
4229 /*Compute SurfaceRelVelMisfit:
4230 *
4231 * 1 [ |u| + eps 2 |v| + eps 2 ]
4232 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
4233 * 2 [ |u |+ eps |v |+ eps ]
4234 * obs obs
4235 */
4236 misfit=0.5*pow(meanvel,2.)*(
4237 pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2.) +
4238 pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2.) );
4239
4240 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVxVyMisfitEnum);
4241
4242 /*Add to cost function*/
4243 Jelem+=misfit*weight*Jdet*gauss->weight;
4244 }
4245
4246 /*clean-up and Return: */
4247 delete gauss;
4248 return Jelem;
4249}
4250/*}}}*/
4251/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/
4252double Tria::SurfaceAbsVelMisfit(bool process_units,int weight_index){
4253
4254 const int numdof=NDOF2*NUMVERTICES;
4255
4256 int i,ig;
4257 double Jelem=0;
4258 double misfit,Jdet;
4259 double vx,vy,vxobs,vyobs,weight;
4260 double xyz_list[NUMVERTICES][3];
4261 GaussTria *gauss=NULL;
4262
4263 /*If on water, return 0: */
4264 if(IsOnWater())return 0;
4265
4266 /* Get node coordinates and dof list: */
4267 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4268
4269 /*Retrieve all inputs we will be needing: */
4270 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4271 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4272 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4273 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4274 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4275
4276 /* Start looping on the number of gaussian points: */
4277 gauss=new GaussTria(2);
4278 for (ig=gauss->begin();ig<gauss->end();ig++){
4279
4280 gauss->GaussPoint(ig);
4281
4282 /* Get Jacobian determinant: */
4283 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4284
4285 /*Get all parameters at gaussian point*/
4286 weights_input->GetInputValue(&weight,gauss,weight_index);
4287 vx_input->GetInputValue(&vx,gauss);
4288 vy_input->GetInputValue(&vy,gauss);
4289 vxobs_input->GetInputValue(&vxobs,gauss);
4290 vyobs_input->GetInputValue(&vyobs,gauss);
4291
4292 /*Compute SurfaceAbsVelMisfitEnum:
4293 *
4294 * 1 [ 2 2 ]
4295 * J = --- | (u - u ) + (v - v ) |
4296 * 2 [ obs obs ]
4297 *
4298 */
4299 misfit=0.5*( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) );
4300
4301 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
4302
4303 /*Add to cost function*/
4304 Jelem+=misfit*weight*Jdet*gauss->weight;
4305 }
4306
4307 /*clean up and Return: */
4308 delete gauss;
4309 return Jelem;
4310}
4311/*}}}*/
4312/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
4313double Tria::SurfaceRelVelMisfit(bool process_units,int weight_index){
4314 const int numdof=2*NUMVERTICES;
4315
4316 int i,ig;
4317 double Jelem=0;
4318 double scalex=1,scaley=1;
4319 double misfit,Jdet;
4320 double epsvel=2.220446049250313e-16;
4321 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4322 double vx,vy,vxobs,vyobs,weight;
4323 double xyz_list[NUMVERTICES][3];
4324 GaussTria *gauss=NULL;
4325
4326 /*If on water, return 0: */
4327 if(IsOnWater())return 0;
4328
4329 /* Get node coordinates and dof list: */
4330 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4331
4332 /*Retrieve all inputs we will be needing: */
4333 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4334 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4335 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4336 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4337 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4338
4339 /* Start looping on the number of gaussian points: */
4340 gauss=new GaussTria(4);
4341 for (ig=gauss->begin();ig<gauss->end();ig++){
4342
4343 gauss->GaussPoint(ig);
4344
4345 /* Get Jacobian determinant: */
4346 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4347
4348 /*Get all parameters at gaussian point*/
4349 weights_input->GetInputValue(&weight,gauss,weight_index);
4350 vx_input->GetInputValue(&vx,gauss);
4351 vy_input->GetInputValue(&vy,gauss);
4352 vxobs_input->GetInputValue(&vxobs,gauss);
4353 vyobs_input->GetInputValue(&vyobs,gauss);
4354
4355 /*Compute SurfaceRelVelMisfit:
4356 *
4357 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
4358 * J = --- | ------------- (u - u ) + ------------- (v - v ) |
4359 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
4360 * obs obs
4361 */
4362 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
4363 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
4364 misfit=0.5*(scalex*pow((vx-vxobs),2.)+scaley*pow((vy-vyobs),2.));
4365 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceRelVelMisfitEnum);
4366
4367 /*Add to cost function*/
4368 Jelem+=misfit*weight*Jdet*gauss->weight;
4369 }
4370
4371 /*clean up and Return: */
4372 delete gauss;
4373 return Jelem;
4374}
4375/*}}}*/
4376/*FUNCTION Tria::ThicknessAbsGradient{{{1*/
4377double Tria::ThicknessAbsGradient(bool process_units,int weight_index){
4378
4379 /* Intermediaries */
4380 int ig;
4381 double Jelem = 0;
4382 double weight;
4383 double Jdet;
4384 double xyz_list[NUMVERTICES][3];
4385 double dp[NDOF2];
4386 GaussTria *gauss = NULL;
4387
4388 /*retrieve parameters and inputs*/
4389
4390 /*If on water, return 0: */
4391 if(IsOnWater()) return 0;
4392
4393 /*Retrieve all inputs we will be needing: */
4394 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4395 Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4396 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4397
4398 /* Start looping on the number of gaussian points: */
4399 gauss=new GaussTria(2);
4400 for (ig=gauss->begin();ig<gauss->end();ig++){
4401
4402 gauss->GaussPoint(ig);
4403
4404 /* Get Jacobian determinant: */
4405 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4406
4407 /*Get all parameters at gaussian point*/
4408 weights_input->GetInputValue(&weight,gauss,weight_index);
4409 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
4410
4411 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
4412 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
4413 }
4414
4415 /*Clean up and return*/
4416 delete gauss;
4417 return Jelem;
4418}
4419/*}}}*/
4420/*FUNCTION Tria::ThicknessAbsMisfit {{{1*/
4421double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){
4422
4423 /*Intermediaries*/
4424 int i,ig;
4425 double thickness,thicknessobs,weight;
4426 double Jdet;
4427 double Jelem = 0;
4428 double xyz_list[NUMVERTICES][3];
4429 GaussTria *gauss = NULL;
4430 double dH[2];
4431
4432 /*If on water, return 0: */
4433 if(IsOnWater())return 0;
4434
4435 /*Retrieve all inputs we will be needing: */
4436 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4437 Input* thickness_input =inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4438 Input* thicknessobs_input=inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
4439 Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4440
4441 /* Start looping on the number of gaussian points: */
4442 gauss=new GaussTria(2);
4443 for (ig=gauss->begin();ig<gauss->end();ig++){
4444
4445 gauss->GaussPoint(ig);
4446
4447 /* Get Jacobian determinant: */
4448 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4449
4450 /*Get parameters at gauss point*/
4451 thickness_input->GetInputValue(&thickness,gauss);
4452 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
4453 thicknessobs_input->GetInputValue(&thicknessobs,gauss);
4454 weights_input->GetInputValue(&weight,gauss,weight_index);
4455
4456 /*compute ThicknessAbsMisfit*/
4457 Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight;
4458 }
4459
4460 /* clean up and Return: */
4461 delete gauss;
4462 return Jelem;
4463}
4464/*}}}*/
4465/*FUNCTION Tria::CreatePVectorAdjointBalancethickness{{{1*/
4466ElementVector* Tria::CreatePVectorAdjointBalancethickness(void){
4467
4468 /*Constants*/
4469 const int numdof=1*NUMVERTICES;
4470
4471 /*Intermediaries */
4472 int i,ig,resp;
4473 double Jdet;
4474 double thickness,thicknessobs,weight;
4475 int *responses = NULL;
4476 int num_responses;
4477 double xyz_list[NUMVERTICES][3];
4478 double basis[3];
4479 double dbasis[NDOF2][NUMVERTICES];
4480 double dH[2];
4481 GaussTria* gauss=NULL;
4482
4483 /*Initialize Element vector*/
4484 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
4485
4486 /*Retrieve all inputs and parameters*/
4487 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4488 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
4489 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
4490 Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4491 Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
4492 Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4493
4494 /* Start looping on the number of gaussian points: */
4495 gauss=new GaussTria(2);
4496 for(ig=gauss->begin();ig<gauss->end();ig++){
4497
4498 gauss->GaussPoint(ig);
4499
4500 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4501 GetNodalFunctions(basis, gauss);
4502 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
4503
4504 thickness_input->GetInputValue(&thickness, gauss);
4505 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
4506 thicknessobs_input->GetInputValue(&thicknessobs, gauss);
4507
4508 /*Loop over all requested responses*/
4509 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
4510
4511 case ThicknessAbsMisfitEnum:
4512 weights_input->GetInputValue(&weight, gauss,resp);
4513 for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
4514 break;
4515 case ThicknessAbsGradientEnum:
4516 weights_input->GetInputValue(&weight, gauss,resp);
4517 for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
4518 for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
4519 break;
4520 default:
4521 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
4522 }
4523 }
4524
4525 /*Clean up and return*/
4526 delete gauss;
4527 xfree((void**)&responses);
4528 return pe;
4529}
4530/*}}}*/
4531/*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/
4532ElementVector* Tria::CreatePVectorAdjointHoriz(void){
4533
4534 /*Constants*/
4535 const int numdof=NDOF2*NUMVERTICES;
4536
4537 /*Intermediaries */
4538 int i,resp,ig;
4539 int *responses=NULL;
4540 int num_responses;
4541 double Jdet;
4542 double obs_velocity_mag,velocity_mag;
4543 double dux,duy;
4544 double epsvel=2.220446049250313e-16;
4545 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4546 double scalex=0,scaley=0,scale=0,S=0;
4547 double vx,vy,vxobs,vyobs,weight;
4548 double xyz_list[NUMVERTICES][3];
4549 double basis[3];
4550 GaussTria* gauss=NULL;
4551
4552 /*Initialize Element vector*/
4553 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
4554
4555 /*Retrieve all inputs and parameters*/
4556 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4557 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
4558 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
4559 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4560 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4561 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4562 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4563 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4564
4565 /*Get Surface if required by one response*/
4566 for(resp=0;resp<num_responses;resp++){
4567 if(responses[resp]==SurfaceAverageVelMisfitEnum){
4568 inputs->GetInputValue(&S,SurfaceAreaEnum); break;
4569 }
4570 }
4571
4572 /* Start looping on the number of gaussian points: */
4573 gauss=new GaussTria(4);
4574 for (ig=gauss->begin();ig<gauss->end();ig++){
4575
4576 gauss->GaussPoint(ig);
4577
4578 /* Get Jacobian determinant: */
4579 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4580
4581 /*Get all parameters at gaussian point*/
4582 vx_input->GetInputValue(&vx,gauss);
4583 vy_input->GetInputValue(&vy,gauss);
4584 vxobs_input->GetInputValue(&vxobs,gauss);
4585 vyobs_input->GetInputValue(&vyobs,gauss);
4586 GetNodalFunctions(basis, gauss);
4587
4588 /*Loop over all requested responses*/
4589 for(resp=0;resp<num_responses;resp++){
4590
4591 weights_input->GetInputValue(&weight,gauss,resp);
4592
4593 switch(responses[resp]){
4594 case SurfaceAbsVelMisfitEnum:
4595 /*
4596 * 1 [ 2 2 ]
4597 * J = --- | (u - u ) + (v - v ) |
4598 * 2 [ obs obs ]
4599 *
4600 * dJ
4601 * DU = - -- = (u - u )
4602 * du obs
4603 */
4604 for (i=0;i<NUMVERTICES;i++){
4605 dux=vxobs-vx;
4606 duy=vyobs-vy;
4607 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4608 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4609 }
4610 break;
4611 case SurfaceRelVelMisfitEnum:
4612 /*
4613 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
4614 * J = --- | ------------- (u - u ) + ------------- (v - v ) |
4615 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
4616 * obs obs
4617 *
4618 * dJ \bar{v}^2
4619 * DU = - -- = ------------- (u - u )
4620 * du (u + eps)^2 obs
4621 * obs
4622 */
4623 for (i=0;i<NUMVERTICES;i++){
4624 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
4625 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
4626 dux=scalex*(vxobs-vx);
4627 duy=scaley*(vyobs-vy);
4628 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4629 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4630 }
4631 break;
4632 case SurfaceLogVelMisfitEnum:
4633 /*
4634 * [ vel + eps ] 2
4635 * J = 4 \bar{v}^2 | log ( ----------- ) |
4636 * [ vel + eps ]
4637 * obs
4638 *
4639 * dJ 2 * log(...)
4640 * DU = - -- = - 4 \bar{v}^2 ------------- u
4641 * du vel^2 + eps
4642 *
4643 */
4644 for (i=0;i<NUMVERTICES;i++){
4645 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel;
4646 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
4647 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
4648 dux=scale*vx;
4649 duy=scale*vy;
4650 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4651 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4652 }
4653 break;
4654 case SurfaceAverageVelMisfitEnum:
4655 /*
4656 * 1 2 2
4657 * J = --- sqrt( (u - u ) + (v - v ) )
4658 * S obs obs
4659 *
4660 * dJ 1 1
4661 * DU = - -- = - --- ----------- * 2 (u - u )
4662 * du S 2 sqrt(...) obs
4663 */
4664 for (i=0;i<NUMVERTICES;i++){
4665 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
4666 dux=scale*(vxobs-vx);
4667 duy=scale*(vyobs-vy);
4668 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4669 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4670 }
4671 break;
4672 case SurfaceLogVxVyMisfitEnum:
4673 /*
4674 * 1 [ |u| + eps 2 |v| + eps 2 ]
4675 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
4676 * 2 [ |u |+ eps |v |+ eps ]
4677 * obs obs
4678 * dJ 1 u 1
4679 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------
4680 * du |u| + eps |u| u + eps
4681 */
4682 for (i=0;i<NUMVERTICES;i++){
4683 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
4684 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
4685 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4686 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4687 }
4688 break;
4689 case DragCoefficientAbsGradientEnum:
4690 /*Nothing in P vector*/
4691 break;
4692 case ThicknessAbsGradientEnum:
4693 /*Nothing in P vector*/
4694 break;
4695 case RheologyBbarAbsGradientEnum:
4696 /*Nothing in P vector*/
4697 break;
4698 default:
4699 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
4700 }
4701 }
4702 }
4703
4704 /*Clean up and return*/
4705 delete gauss;
4706 xfree((void**)&responses);
4707 return pe;
4708}
4709/*}}}*/
4710/*FUNCTION Tria::CreatePVectorAdjointStokes{{{1*/
4711ElementVector* Tria::CreatePVectorAdjointStokes(void){
4712
4713 /*Intermediaries */
4714 int i,resp,ig;
4715 int *responses=NULL;
4716 int num_responses;
4717 double Jdet;
4718 double obs_velocity_mag,velocity_mag;
4719 double dux,duy;
4720 double epsvel=2.220446049250313e-16;
4721 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4722 double scalex=0,scaley=0,scale=0,S=0;
4723 double vx,vy,vxobs,vyobs,weight;
4724 double xyz_list[NUMVERTICES][3];
4725 double basis[3];
4726 GaussTria* gauss=NULL;
4727
4728 /*Initialize Element vector*/
4729 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
4730
4731 /*Retrieve all inputs and parameters*/
4732 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4733 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
4734 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
4735 Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4736 Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
4737 Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
4738 Input* vxobs_input = inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4739 Input* vyobs_input = inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4740
4741 /*Get Surface if required by one response*/
4742 for(resp=0;resp<num_responses;resp++){
4743 if(responses[resp]==SurfaceAverageVelMisfitEnum){
4744 inputs->GetInputValue(&S,SurfaceAreaEnum); break;
4745 }
4746 }
4747
4748 /* Start looping on the number of gaussian points: */
4749 gauss=new GaussTria(4);
4750 for (ig=gauss->begin();ig<gauss->end();ig++){
4751
4752 gauss->GaussPoint(ig);
4753
4754 /* Get Jacobian determinant: */
4755 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4756
4757 /*Get all parameters at gaussian point*/
4758 vx_input->GetInputValue(&vx,gauss);
4759 vy_input->GetInputValue(&vy,gauss);
4760 vxobs_input->GetInputValue(&vxobs,gauss);
4761 vyobs_input->GetInputValue(&vyobs,gauss);
4762 GetNodalFunctions(basis, gauss);
4763
4764 /*Loop over all requested responses*/
4765 for(resp=0;resp<num_responses;resp++){
4766
4767 weights_input->GetInputValue(&weight,gauss,resp);
4768
4769 switch(responses[resp]){
4770
4771 case SurfaceAbsVelMisfitEnum:
4772 /*
4773 * 1 [ 2 2 ]
4774 * J = --- | (u - u ) + (v - v ) |
4775 * 2 [ obs obs ]
4776 *
4777 * dJ
4778 * DU = - -- = (u - u )
4779 * du obs
4780 */
4781 for (i=0;i<NUMVERTICES;i++){
4782 dux=vxobs-vx;
4783 duy=vyobs-vy;
4784 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4785 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4786 }
4787 break;
4788 case SurfaceRelVelMisfitEnum:
4789 /*
4790 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
4791 * J = --- | ------------- (u - u ) + ------------- (v - v ) |
4792 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
4793 * obs obs
4794 *
4795 * dJ \bar{v}^2
4796 * DU = - -- = ------------- (u - u )
4797 * du (u + eps)^2 obs
4798 * obs
4799 */
4800 for (i=0;i<NUMVERTICES;i++){
4801 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
4802 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
4803 dux=scalex*(vxobs-vx);
4804 duy=scaley*(vyobs-vy);
4805 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4806 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4807 }
4808 break;
4809 case SurfaceLogVelMisfitEnum:
4810 /*
4811 * [ vel + eps ] 2
4812 * J = 4 \bar{v}^2 | log ( ----------- ) |
4813 * [ vel + eps ]
4814 * obs
4815 *
4816 * dJ 2 * log(...)
4817 * DU = - -- = - 4 \bar{v}^2 ------------- u
4818 * du vel^2 + eps
4819 *
4820 */
4821 for (i=0;i<NUMVERTICES;i++){
4822 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel;
4823 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
4824 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
4825 dux=scale*vx;
4826 duy=scale*vy;
4827 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4828 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4829 }
4830 break;
4831 case SurfaceAverageVelMisfitEnum:
4832 /*
4833 * 1 2 2
4834 * J = --- sqrt( (u - u ) + (v - v ) )
4835 * S obs obs
4836 *
4837 * dJ 1 1
4838 * DU = - -- = - --- ----------- * 2 (u - u )
4839 * du S 2 sqrt(...) obs
4840 */
4841 for (i=0;i<NUMVERTICES;i++){
4842 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
4843 dux=scale*(vxobs-vx);
4844 duy=scale*(vyobs-vy);
4845 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4846 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4847 }
4848 break;
4849 case SurfaceLogVxVyMisfitEnum:
4850 /*
4851 * 1 [ |u| + eps 2 |v| + eps 2 ]
4852 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
4853 * 2 [ |u |+ eps |v |+ eps ]
4854 * obs obs
4855 * dJ 1 u 1
4856 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------
4857 * du |u| + eps |u| u + eps
4858 */
4859 for (i=0;i<NUMVERTICES;i++){
4860 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
4861 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
4862 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4863 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4864 }
4865 break;
4866 case DragCoefficientAbsGradientEnum:
4867 /*Nothing in P vector*/
4868 break;
4869 case ThicknessAbsGradientEnum:
4870 /*Nothing in P vector*/
4871 break;
4872 case RheologyBbarAbsGradientEnum:
4873 /*Nothing in P vector*/
4874 break;
4875 default:
4876 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
4877 }
4878 }
4879 }
4880
4881 /*Clean up and return*/
4882 delete gauss;
4883 xfree((void**)&responses);
4884 return pe;
4885}
4886/*}}}*/
4887/*FUNCTION Tria::DragCoefficientAbsGradient{{{1*/
4888double Tria::DragCoefficientAbsGradient(bool process_units,int weight_index){
4889
4890 /* Intermediaries */
4891 int ig;
4892 double Jelem = 0;
4893 double weight;
4894 double Jdet;
4895 double xyz_list[NUMVERTICES][3];
4896 double dp[NDOF2];
4897 GaussTria *gauss = NULL;
4898
4899 /*retrieve parameters and inputs*/
4900
4901 /*If on water, return 0: */
4902 if(IsOnWater()) return 0;
4903
4904 /*Retrieve all inputs we will be needing: */
4905 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4906 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4907 Input* drag_input =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input);
4908
4909 /* Start looping on the number of gaussian points: */
4910 gauss=new GaussTria(2);
4911 for (ig=gauss->begin();ig<gauss->end();ig++){
4912
4913 gauss->GaussPoint(ig);
4914
4915 /* Get Jacobian determinant: */
4916 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4917
4918 /*Get all parameters at gaussian point*/
4919 weights_input->GetInputValue(&weight,gauss,weight_index);
4920 drag_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
4921
4922 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
4923 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
4924 }
4925
4926 /*Clean up and return*/
4927 delete gauss;
4928 return Jelem;
4929}
4930/*}}}*/
4931/*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{1*/
4932ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){
4933
4934 ElementMatrix* Ke=NULL;
4935
4936 /*Get Element Matrix of the forward model*/
4937 switch(GetElementType()){
4938 case P1Enum:
4939 Ke=CreateKMatrixBalancethickness_CG();
4940 break;
4941 case P1DGEnum:
4942 Ke=CreateKMatrixBalancethickness_DG();
4943 break;
4944 default:
4945 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
4946 }
4947
4948 /*Transpose and return Ke*/
4949 Ke->Transpose();
4950 return Ke;
4951}
4952/*}}}*/
4953/*FUNCTION Tria::CreateKMatrixAdjointMacAyeal{{{1*/
4954ElementMatrix* Tria::CreateKMatrixAdjointMacAyeal(void){
4955
4956 /*Constants*/
4957 const int numdof=NDOF2*NUMVERTICES;
4958
4959 /*Intermediaries */
4960 int i,j,ig;
4961 bool incomplete_adjoint;
4962 double xyz_list[NUMVERTICES][3];
4963 double Jdet,thickness;
4964 double eps1dotdphii,eps1dotdphij;
4965 double eps2dotdphii,eps2dotdphij;
4966 double mu_prime;
4967 double epsilon[3];/* epsilon=[exx,eyy,exy];*/
4968 double eps1[2],eps2[2];
4969 double phi[NUMVERTICES];
4970 double dphi[2][NUMVERTICES];
4971 GaussTria *gauss=NULL;
4972
4973 /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
4974 parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
4975 ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
4976 if(incomplete_adjoint) return Ke;
4977
4978 /*Retrieve all inputs and parameters*/
4979 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4980 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
4981 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
4982 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4983
4984 /* Start looping on the number of gaussian points: */
4985 gauss=new GaussTria(2);
4986 for (ig=gauss->begin();ig<gauss->end();ig++){
4987
4988 gauss->GaussPoint(ig);
4989
4990 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4991 GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
4992
4993 thickness_input->GetInputValue(&thickness, gauss);
4994 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
4995 matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
4996 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
4997 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
4998
4999 for(i=0;i<3;i++){
5000 for(j=0;j<3;j++){
5001 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
5002 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
5003 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
5004 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
5005
5006 Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
5007 Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
5008 Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
5009 Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
5010 }
5011 }
5012 }
5013
5014 /*Transform Coordinate System*/
5015 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
5016
5017 /*Clean up and return*/
5018 delete gauss;
5019 //Ke->Transpose();
5020 return Ke;
5021}
5022/*}}}*/
5023/*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/
5024void Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){
5025
5026 const int numdof=NDOF2*NUMVERTICES;
5027
5028 int i;
5029 int* doflist=NULL;
5030 double values[numdof];
5031 double lambdax[NUMVERTICES];
5032 double lambday[NUMVERTICES];
5033
5034 /*Get dof list: */
5035 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5036
5037 /*Use the dof list to index into the solution vector: */
5038 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
5039
5040 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
5041 for(i=0;i<NUMVERTICES;i++){
5042 lambdax[i]=values[i*NDOF2+0];
5043 lambday[i]=values[i*NDOF2+1];
5044
5045 /*Check solution*/
5046 if(isnan(lambdax[i])) _error_("NaN found in solution vector");
5047 if(isnan(lambday[i])) _error_("NaN found in solution vector");
5048 }
5049
5050 /*Add vx and vy as inputs to the tria element: */
5051 this->inputs->AddInput(new TriaP1Input(AdjointxEnum,lambdax));
5052 this->inputs->AddInput(new TriaP1Input(AdjointyEnum,lambday));
5053
5054 /*Free ressources:*/
5055 xfree((void**)&doflist);
5056}
5057/*}}}*/
5058
5059/*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{1*/
5060void Tria::InputUpdateFromSolutionAdjointBalancethickness(double* solution){
5061
5062 const int numdof=NDOF1*NUMVERTICES;
5063
5064 int i;
5065 int* doflist=NULL;
5066 double values[numdof];
5067 double lambda[NUMVERTICES];
5068
5069 /*Get dof list: */
5070 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5071
5072 /*Use the dof list to index into the solution vector: */
5073 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
5074
5075 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
5076 for(i=0;i<numdof;i++){
5077 lambda[i]=values[i];
5078 if(isnan(lambda[i])) _error_("NaN found in solution vector");
5079 }
5080
5081 /*Add vx and vy as inputs to the tria element: */
5082 this->inputs->AddInput(new TriaP1Input(AdjointEnum,lambda));
5083
5084 /*Free ressources:*/
5085 xfree((void**)&doflist);
5086}
5087/*}}}*/
5088/*FUNCTION Tria::GetVectorFromControlInputs{{{1*/
5089void Tria::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data){
5090
5091 int doflist1[NUMVERTICES];
5092 Input *input=NULL;
5093
5094 /*Get out if this is not an element input*/
5095 if(!IsInput(control_enum)) return;
5096
5097 /*Prepare index list*/
5098 GradientIndexing(&doflist1[0],control_index);
5099
5100 /*Get input (either in element or material)*/
5101 if(control_enum==MaterialsRheologyBbarEnum){
5102 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
5103 }
5104 else{
5105 input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
5106 }
5107
5108 /*Check that it is a ControlInput*/
5109 if (input->ObjectEnum()!=ControlInputEnum){
5110 _error_("input %s is not a ControlInput",EnumToStringx(control_enum));
5111 }
5112
5113 ((ControlInput*)input)->GetVectorFromInputs(vector,&doflist1[0],data);
5114}
5115/*}}}*/
5116/*FUNCTION Tria::SetControlInputsFromVector{{{1*/
5117void Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){
5118
5119 double values[NUMVERTICES];
5120 int doflist1[NUMVERTICES];
5121 Input *input = NULL;
5122 Input *new_input = NULL;
5123
5124 /*Get out if this is not an element input*/
5125 if(!IsInput(control_enum)) return;
5126
5127 /*Prepare index list*/
5128 GradientIndexing(&doflist1[0],control_index);
5129
5130 /*Get values on vertices*/
5131 for (int i=0;i<NUMVERTICES;i++){
5132 values[i]=vector[doflist1[i]];
5133 }
5134 new_input = new TriaP1Input(control_enum,values);
5135
5136 if(control_enum==MaterialsRheologyBbarEnum){
5137 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
5138 }
5139 else{
5140 input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
5141 }
5142
5143 if (input->ObjectEnum()!=ControlInputEnum){
5144 _error_("input %s is not a ControlInput",EnumToStringx(control_enum));
5145 }
5146
5147 ((ControlInput*)input)->SetInput(new_input);
5148}
5149/*}}}*/
5150#endif
5151
5152#ifdef _HAVE_HYDROLOGY_
5153/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{1*/
5154void Tria::CreateHydrologyWaterVelocityInput(void){
5155
5156 /*material parameters: */
5157 double mu_water;
5158 double VelocityFactor; // This factor represents the number 12 in laminar flow velocity which can vary by differnt hydrology.CR
5159 double n_man,CR;
5160 double w;
5161 double rho_ice, rho_water, g;
5162 double dsdx,dsdy,dbdx,dbdy;
5163 double vx[NUMVERTICES];
5164 double vy[NUMVERTICES];
5165 GaussTria *gauss = NULL;
5166
5167 /*Retrieve all inputs and parameters*/
5168 rho_ice=matpar->GetRhoIce();
5169 rho_water=matpar->GetRhoWater();
5170 g=matpar->GetG();
5171 CR=matpar->GetHydrologyCR(); // To have Lebrocq equavalent equation: CR=0.01,n_man=0.02
5172 n_man=matpar->GetHydrologyN();
5173 mu_water=matpar->GetMuWater();
5174 Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
5175 Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
5176 Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(bedslopex_input);
5177 Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(bedslopey_input);
5178 Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
5179
5180 /* compute VelocityFactor */
5181 VelocityFactor= n_man*pow(CR,2)*rho_water*g/mu_water;
5182
5183 gauss=new GaussTria();
5184 for (int iv=0;iv<NUMVERTICES;iv++){
5185 gauss->GaussVertex(iv);
5186 surfaceslopex_input->GetInputValue(&dsdx,gauss);
5187 surfaceslopey_input->GetInputValue(&dsdy,gauss);
5188 bedslopex_input->GetInputValue(&dbdx,gauss);
5189 bedslopey_input->GetInputValue(&dbdy,gauss);
5190 watercolumn_input->GetInputValue(&w,gauss);
5191
5192 /* Water velocity x and y components */
5193 // vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
5194 // vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
5195
5196 vx[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
5197 vy[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
5198 }
5199
5200 /*clean-up*/
5201 delete gauss;
5202
5203 /*Add to inputs*/
5204 this->inputs->AddInput(new TriaP1Input(HydrologyWaterVxEnum,vx));
5205 this->inputs->AddInput(new TriaP1Input(HydrologyWaterVyEnum,vy));
5206}
5207/*}}}*/
5208/*FUNCTION Tria::CreateKMatrixHydrology{{{1*/
5209ElementMatrix* Tria::CreateKMatrixHydrology(void){
5210
5211 /*Constants*/
5212 const int numdof=NDOF1*NUMVERTICES;
5213
5214 /*Intermediaries */
5215 double diffusivity;
5216 int i,j,ig;
5217 double Jdettria,DL_scalar,dt,h;
5218 double vx,vy,vel,dvxdx,dvydy;
5219 double dvx[2],dvy[2];
5220 double v_gauss[2]={0.0};
5221 double xyz_list[NUMVERTICES][3];
5222 double L[NUMVERTICES];
5223 double B[2][NUMVERTICES];
5224 double Bprime[2][NUMVERTICES];
5225 double K[2][2] ={0.0};
5226 double KDL[2][2] ={0.0};
5227 double DL[2][2] ={0.0};
5228 double DLprime[2][2] ={0.0};
5229 GaussTria *gauss=NULL;
5230
5231 /*Skip if water or ice shelf element*/
5232 if(IsOnWater() | IsFloating()) return NULL;
5233
5234 /*Initialize Element matrix*/
5235 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
5236
5237 /*Create water velocity vx and vy from current inputs*/
5238 CreateHydrologyWaterVelocityInput();
5239
5240 /*Retrieve all inputs and parameters*/
5241 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5242 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
5243 this->parameters->FindParam(&diffusivity,HydrologyStabilizationEnum);
5244 Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
5245 Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
5246 h=sqrt(2*this->GetArea());
5247
5248 /* Start looping on the number of gaussian points: */
5249 gauss=new GaussTria(2);
5250 for (ig=gauss->begin();ig<gauss->end();ig++){
5251
5252 gauss->GaussPoint(ig);
5253
5254 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5255 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
5256
5257 vx_input->GetInputValue(&vx,gauss);
5258 vy_input->GetInputValue(&vy,gauss);
5259 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
5260 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
5261
5262 DL_scalar=gauss->weight*Jdettria;
5263
5264 TripleMultiply( &L[0],1,numdof,1,
5265 &DL_scalar,1,1,0,
5266 &L[0],1,numdof,0,
5267 &Ke->values[0],1);
5268
5269 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
5270 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
5271
5272 dvxdx=dvx[0];
5273 dvydy=dvy[1];
5274 DL_scalar=dt*gauss->weight*Jdettria;
5275
5276 DL[0][0]=DL_scalar*dvxdx;
5277 DL[1][1]=DL_scalar*dvydy;
5278 DLprime[0][0]=DL_scalar*vx;
5279 DLprime[1][1]=DL_scalar*vy;
5280
5281 TripleMultiply( &B[0][0],2,numdof,1,
5282 &DL[0][0],2,2,0,
5283 &B[0][0],2,numdof,0,
5284 &Ke->values[0],1);
5285
5286 TripleMultiply( &B[0][0],2,numdof,1,
5287 &DLprime[0][0],2,2,0,
5288 &Bprime[0][0],2,numdof,0,
5289 &Ke->values[0],1);
5290
5291 /*Artificial diffusivity*/
5292 vel=sqrt(pow(vx,2.)+pow(vy,2.));
5293 K[0][0]=diffusivity*h/(2*vel)*vx*vx;
5294 K[1][0]=diffusivity*h/(2*vel)*vy*vx;
5295 K[0][1]=diffusivity*h/(2*vel)*vx*vy;
5296 K[1][1]=diffusivity*h/(2*vel)*vy*vy;
5297 KDL[0][0]=DL_scalar*K[0][0];
5298 KDL[1][0]=DL_scalar*K[1][0];
5299 KDL[0][1]=DL_scalar*K[0][1];
5300 KDL[1][1]=DL_scalar*K[1][1];
5301
5302 TripleMultiply( &Bprime[0][0],2,numdof,1,
5303 &KDL[0][0],2,2,0,
5304 &Bprime[0][0],2,numdof,0,
5305 &Ke->values[0],1);
5306 }
5307
5308 /*Clean up and return*/
5309 delete gauss;
5310 return Ke;
5311}
5312/*}}}*/
5313/*FUNCTION Tria::CreatePVectorHydrology {{{1*/
5314ElementVector* Tria::CreatePVectorHydrology(void){
5315
5316 /*Constants*/
5317 const int numdof=NDOF1*NUMVERTICES;
5318
5319 /*Intermediaries */
5320 int i,j,ig;
5321 double Jdettria,dt;
5322 double basal_melting_g;
5323 double old_watercolumn_g;
5324 double xyz_list[NUMVERTICES][3];
5325 double basis[numdof];
5326 GaussTria* gauss=NULL;
5327
5328 /*Skip if water or ice shelf element*/
5329 if(IsOnWater() | IsFloating()) return NULL;
5330
5331 /*Initialize Element vector*/
5332 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
5333
5334 /*Retrieve all inputs and parameters*/
5335 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5336 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
5337 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
5338 Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum); _assert_(old_watercolumn_input);
5339
5340 /*Initialize basal_melting_correction_g to 0, do not forget!:*/
5341 /* Start looping on the number of gaussian points: */
5342 gauss=new GaussTria(2);
5343 for(ig=gauss->begin();ig<gauss->end();ig++){
5344
5345 gauss->GaussPoint(ig);
5346
5347 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5348 GetNodalFunctions(basis, gauss);
5349
5350 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
5351 old_watercolumn_input->GetInputValue(&old_watercolumn_g,gauss);
5352
5353 if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
5354 else for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
5355 }
5356
5357 /*Clean up and return*/
5358 delete gauss;
5359 return pe;
5360}
5361/*}}}*/
5362/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
5363void Tria::GetSolutionFromInputsHydrology(Vector* solution){
5364
5365 const int numdof=NDOF1*NUMVERTICES;
5366
5367 int i;
5368 int* doflist=NULL;
5369 double watercolumn;
5370 double values[numdof];
5371 GaussTria* gauss=NULL;
5372
5373 /*Get dof list: */
5374 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5375
5376 /*Get inputs*/
5377 Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
5378
5379 /*Ok, we have watercolumn values, fill in watercolumn array: */
5380 /*P1 element only for now*/
5381 gauss=new GaussTria();
5382 for(i=0;i<NUMVERTICES;i++){
5383
5384 gauss->GaussVertex(i);
5385
5386 /*Recover watercolumn*/
5387 watercolumn_input->GetInputValue(&watercolumn,gauss);
5388 values[i]=watercolumn;
5389 }
5390
5391 solution->SetValues(numdof,doflist,values,INS_VAL);
5392
5393 /*Free ressources:*/
5394 delete gauss;
5395 xfree((void**)&doflist);
5396}
5397/*}}}*/
5398/*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{1*/
5399void Tria::InputUpdateFromSolutionHydrology(double* solution){
5400
5401 /*Intermediaries*/
5402 const int numdof = NDOF1*NUMVERTICES;
5403
5404 int i;
5405 int* doflist=NULL;
5406 double values[numdof];
5407
5408 /*Get dof list: */
5409 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5410
5411 /*Use the dof list to index into the solution vector: */
5412 for(i=0;i<numdof;i++){
5413 values[i]=solution[doflist[i]];
5414 if(isnan(values[i])) _error_("NaN found in solution vector");
5415 if (values[i]<pow((double)10,(double)-10))values[i]=pow((double)10,(double)-10); //correcting the water column to positive values
5416
5417 }
5418
5419 /*Add input to the element: */
5420 this->inputs->AddInput(new TriaP1Input(WatercolumnEnum,values));
5421
5422 /*Free ressources:*/
5423 xfree((void**)&doflist);
5424}
5425/*}}}*/
5426#endif
5427
5428#ifdef _HAVE_DAKOTA_
5429/*FUNCTION Tria::InputUpdateFromVectorDakota(double* vector, int name, int type);{{{1*/
5430void Tria::InputUpdateFromVectorDakota(double* vector, int name, int type){
5431
5432 int i,j;
5433
5434 /*Check that name is an element input*/
5435 if (!IsInput(name)) return;
5436
5437 switch(type){
5438
5439 case VertexEnum:
5440
5441 /*New TriaP1Input*/
5442 double values[3];
5443
5444 /*Get values on the 3 vertices*/
5445 for (i=0;i<3;i++){
5446 values[i]=vector[this->nodes[i]->GetSidList()]; //careful, vector of values here is not parallel distributed, but serial distributed (from a serial Dakota core!)
5447 }
5448
5449 /*Branch on the specified type of update: */
5450 switch(name){
5451 case ThicknessEnum:
5452 /*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium {{{2*/
5453 double thickness[3];
5454 double thickness_init[3];
5455 double hydrostatic_ratio[3];
5456 double surface[3];
5457 double bed[3];
5458
5459 /*retrieve inputs: */
5460 GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
5461 GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
5462 GetInputListOnVertices(&bed[0],BedEnum);
5463 GetInputListOnVertices(&surface[0],SurfaceEnum);
5464
5465 /*build new thickness: */
5466// for(j=0;j<3;j++)thickness[j]=values[j];
5467
5468 /*build new bed and surface: */
5469 if (this->IsFloating()){
5470 /*hydrostatic equilibrium: */
5471 double rho_ice,rho_water,di;
5472 rho_ice=this->matpar->GetRhoIce();
5473 rho_water=this->matpar->GetRhoWater();
5474
5475 di=rho_ice/rho_water;
5476
5477 /*build new thickness: */
5478 for (j=0; j<3; j++) {
5479 /* for observed/interpolated/hydrostatic thickness, remove scaling from any hydrostatic thickness */
5480 if (hydrostatic_ratio[j] >= 0.)
5481 thickness[j]=values[j]-(values[j]/thickness_init[j]-1.)*hydrostatic_ratio[j]*surface[j]/(1.-di);
5482 /* for minimum thickness, don't scale */
5483 else
5484 thickness[j]=thickness_init[j];
5485
5486 /* check the computed thickness and update bed */
5487 if (thickness[j] < 0.)
5488 thickness[j]=1./(1.-di);
5489 bed[j]=surface[j]-thickness[j];
5490 }
5491
5492// for(j=0;j<3;j++){
5493// surface[j]=(1-di)*thickness[j];
5494// bed[j]=-di*thickness[j];
5495// }
5496 }
5497 else{
5498 /*build new thickness: */
5499 for (j=0; j<3; j++) {
5500 /* for observed thickness, use scaled value */
5501 if (hydrostatic_ratio[j] >= 0.)
5502 thickness[j]=values[j];
5503 /* for minimum thickness, don't scale */
5504 else
5505 thickness[j]=thickness_init[j];
5506 }
5507
5508 /*update bed on grounded ice: */
5509// for(j=0;j<3;j++)surface[j]=bed[j]+thickness[j];
5510 for(j=0;j<3;j++)bed[j]=surface[j]-thickness[j];
5511 }
5512
5513 /*Add new inputs: */
5514 this->inputs->AddInput(new TriaP1Input(ThicknessEnum,thickness));
5515 this->inputs->AddInput(new TriaP1Input(BedEnum,bed));
5516 this->inputs->AddInput(new TriaP1Input(SurfaceEnum,surface));
5517
5518 /*}}}*/
5519 break;
5520 default:
5521 this->inputs->AddInput(new TriaP1Input(name,values));
5522 }
5523 break;
5524
5525 default:
5526 _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
5527 }
5528
5529}
5530/*}}}*/
5531/*FUNCTION Tria::InputUpdateFromVectorDakota(int* vector, int name, int type);{{{1*/
5532void Tria::InputUpdateFromVectorDakota(int* vector, int name, int type){
5533 _error_(" not supported yet!");
5534}
5535/*}}}*/
5536/*FUNCTION Tria::InputUpdateFromVectorDakota(bool* vector, int name, int type);{{{1*/
5537void Tria::InputUpdateFromVectorDakota(bool* vector, int name, int type){
5538 _error_(" not supported yet!");
5539}
5540/*}}}*/
5541/*FUNCTION Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type);{{{1*/
5542void Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){
5543
5544 int i,j,t;
5545 TransientInput* transientinput=NULL;
5546 double values[3];
5547 double time;
5548 int row;
5549 double yts;
5550
5551 /*Check that name is an element input*/
5552 if (!IsInput(name)) return;
5553
5554 switch(type){
5555
5556 case VertexEnum:
5557
5558 /*Create transient input: */
5559
5560 parameters->FindParam(&yts,ConstantsYtsEnum);
5561 for(t=0;t<ncols;t++){ //ncols is the number of times
5562
5563 /*create input values: */
5564 for(i=0;i<3;i++){
5565 row=this->nodes[i]->GetSidList();
5566 values[i]=(double)matrix[ncols*row+t];
5567 }
5568
5569 /*time? :*/
5570 time=(double)matrix[(nrows-1)*ncols+t]*yts;
5571
5572 if(t==0) transientinput=new TransientInput(name);
5573 transientinput->AddTimeInput(new TriaP1Input(name,values),time);
5574 transientinput->Configure(parameters);
5575 }
5576 this->inputs->AddInput(transientinput);
5577 break;
5578
5579 default:
5580 _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
5581 }
5582
5583}
5584/*}}}*/
5585#endif
5586
5587#ifdef _HAVE_BALANCED_
5588/*FUNCTION Tria::CreateKMatrixBalancethickness {{{1*/
5589ElementMatrix* Tria::CreateKMatrixBalancethickness(void){
5590
5591 switch(GetElementType()){
5592 case P1Enum:
5593 return CreateKMatrixBalancethickness_CG();
5594 case P1DGEnum:
5595 return CreateKMatrixBalancethickness_DG();
5596 default:
5597 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
5598 }
5599
5600}
5601/*}}}*/
5602/*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{1*/
5603ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){
5604
5605 /*Constants*/
5606 const int numdof=NDOF1*NUMVERTICES;
5607
5608 /*Intermediaries */
5609 int stabilization;
5610 int i,j,ig,dim;
5611 double Jdettria,vx,vy,dvxdx,dvydy,vel,h;
5612 double dvx[2],dvy[2];
5613 double xyz_list[NUMVERTICES][3];
5614 double L[NUMVERTICES];
5615 double B[2][NUMVERTICES];
5616 double Bprime[2][NUMVERTICES];
5617 double K[2][2] = {0.0};
5618 double KDL[2][2] = {0.0};
5619 double DL[2][2] = {0.0};
5620 double DLprime[2][2] = {0.0};
5621 double DL_scalar;
5622 GaussTria *gauss = NULL;
5623
5624 /*Initialize Element matrix*/
5625 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
5626
5627 /*Retrieve all Inputs and parameters: */
5628 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5629 this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum);
5630 this->parameters->FindParam(&dim,MeshDimensionEnum);
5631 Input* vxaverage_input=NULL;
5632 Input* vyaverage_input=NULL;
5633 if(dim==2){
5634 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
5635 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
5636 }
5637 else{
5638 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
5639 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
5640 }
5641 h=sqrt(2*this->GetArea());
5642
5643 /*Start looping on the number of gaussian points:*/
5644 gauss=new GaussTria(2);
5645 for (ig=gauss->begin();ig<gauss->end();ig++){
5646
5647 gauss->GaussPoint(ig);
5648
5649 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5650 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
5651 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
5652
5653 vxaverage_input->GetInputValue(&vx,gauss);
5654 vyaverage_input->GetInputValue(&vy,gauss);
5655 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
5656 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
5657
5658 dvxdx=dvx[0];
5659 dvydy=dvy[1];
5660 DL_scalar=gauss->weight*Jdettria;
5661
5662 DL[0][0]=DL_scalar*dvxdx;
5663 DL[1][1]=DL_scalar*dvydy;
5664
5665 DLprime[0][0]=DL_scalar*vx;
5666 DLprime[1][1]=DL_scalar*vy;
5667
5668 TripleMultiply( &B[0][0],2,numdof,1,
5669 &DL[0][0],2,2,0,
5670 &B[0][0],2,numdof,0,
5671 &Ke->values[0],1);
5672
5673 TripleMultiply( &B[0][0],2,numdof,1,
5674 &DLprime[0][0],2,2,0,
5675 &Bprime[0][0],2,numdof,0,
5676 &Ke->values[0],1);
5677
5678 if(stabilization==1){
5679 /*Streamline upwinding*/
5680 vel=sqrt(pow(vx,2.)+pow(vy,2.));
5681 K[0][0]=h/(2*vel)*vx*vx;
5682 K[1][0]=h/(2*vel)*vy*vx;
5683 K[0][1]=h/(2*vel)*vx*vy;
5684 K[1][1]=h/(2*vel)*vy*vy;
5685 }
5686 else if(stabilization==2){
5687 /*MacAyeal*/
5688 vxaverage_input->GetInputAverage(&vx);
5689 vyaverage_input->GetInputAverage(&vy);
5690 K[0][0]=h/2.0*fabs(vx);
5691 K[0][1]=0.;
5692 K[1][0]=0.;
5693 K[1][1]=h/2.0*fabs(vy);
5694 }
5695 if(stabilization==1 || stabilization==2){
5696 KDL[0][0]=DL_scalar*K[0][0];
5697 KDL[1][0]=DL_scalar*K[1][0];
5698 KDL[0][1]=DL_scalar*K[0][1];
5699 KDL[1][1]=DL_scalar*K[1][1];
5700 TripleMultiply( &Bprime[0][0],2,numdof,1,
5701 &KDL[0][0],2,2,0,
5702 &Bprime[0][0],2,numdof,0,
5703 &Ke->values[0],1);
5704 }
5705 }
5706
5707 /*Clean up and return*/
5708 delete gauss;
5709 return Ke;
5710}
5711/*}}}*/
5712/*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{1*/
5713ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){
5714
5715 /*Constants*/
5716 const int numdof=NDOF1*NUMVERTICES;
5717
5718 /*Intermediaries*/
5719 int i,j,ig,dim;
5720 double vx,vy,Jdettria;
5721 double xyz_list[NUMVERTICES][3];
5722 double B[2][NUMVERTICES];
5723 double Bprime[2][NUMVERTICES];
5724 double DL[2][2]={0.0};
5725 double DL_scalar;
5726 GaussTria *gauss=NULL;
5727
5728 /*Initialize Element matrix*/
5729 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
5730
5731 /*Retrieve all inputs and parameters*/
5732 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5733 this->parameters->FindParam(&dim,MeshDimensionEnum);
5734 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
5735 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
5736
5737 /*Start looping on the number of gaussian points:*/
5738 gauss=new GaussTria(2);
5739 for (ig=gauss->begin();ig<gauss->end();ig++){
5740
5741 gauss->GaussPoint(ig);
5742
5743 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5744 /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
5745 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
5746 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
5747
5748 vx_input->GetInputValue(&vx,gauss);
5749 vy_input->GetInputValue(&vy,gauss);
5750
5751 DL_scalar=-gauss->weight*Jdettria;
5752 DL[0][0]=DL_scalar*vx;
5753 DL[1][1]=DL_scalar*vy;
5754
5755 TripleMultiply( &B[0][0],2,numdof,1,
5756 &DL[0][0],2,2,0,
5757 &Bprime[0][0],2,numdof,0,
5758 &Ke->values[0],1);
5759 }
5760
5761 /*Clean up and return*/
5762 delete gauss;
5763 return Ke;
5764}
5765/*}}}*/
5766/*FUNCTION Tria::CreatePVectorBalancethickness{{{1*/
5767ElementVector* Tria::CreatePVectorBalancethickness(void){
5768
5769 switch(GetElementType()){
5770 case P1Enum:
5771 return CreatePVectorBalancethickness_CG();
5772 break;
5773 case P1DGEnum:
5774 return CreatePVectorBalancethickness_DG();
5775 default:
5776 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
5777 }
5778}
5779/*}}}*/
5780/*FUNCTION Tria::CreatePVectorBalancethickness_CG{{{1*/
5781ElementVector* Tria::CreatePVectorBalancethickness_CG(void){
5782
5783 /*Constants*/
5784 const int numdof=NDOF1*NUMVERTICES;
5785
5786 /*Intermediaries */
5787 int i,j,ig;
5788 double xyz_list[NUMVERTICES][3];
5789 double dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
5790 double L[NUMVERTICES];
5791 GaussTria* gauss=NULL;
5792
5793 /*Initialize Element vector*/
5794 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
5795
5796 /*Retrieve all inputs and parameters*/
5797 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5798 Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
5799 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
5800 Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
5801
5802 /* Start looping on the number of gaussian points: */
5803 gauss=new GaussTria(2);
5804 for(ig=gauss->begin();ig<gauss->end();ig++){
5805
5806 gauss->GaussPoint(ig);
5807
5808 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
5809 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
5810 dhdt_input->GetInputValue(&dhdt_g,gauss);
5811
5812 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5813 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
5814
5815 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
5816 }
5817
5818 /*Clean up and return*/
5819 delete gauss;
5820 return pe;
5821}
5822/*}}}*/
5823/*FUNCTION Tria::CreatePVectorBalancethickness_DG {{{1*/
5824ElementVector* Tria::CreatePVectorBalancethickness_DG(void){
5825
5826 /*Constants*/
5827 const int numdof=NDOF1*NUMVERTICES;
5828
5829 /*Intermediaries */
5830 int i,j,ig;
5831 double xyz_list[NUMVERTICES][3];
5832 double basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
5833 double L[NUMVERTICES];
5834 GaussTria* gauss=NULL;
5835
5836 /*Initialize Element vector*/
5837 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
5838
5839 /*Retrieve all inputs and parameters*/
5840 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5841 Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
5842 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
5843 Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
5844
5845 /* Start looping on the number of gaussian points: */
5846 gauss=new GaussTria(2);
5847 for(ig=gauss->begin();ig<gauss->end();ig++){
5848
5849 gauss->GaussPoint(ig);
5850
5851 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
5852 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
5853 dhdt_input->GetInputValue(&dhdt_g,gauss);
5854
5855 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5856 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
5857
5858 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
5859 }
5860
5861 /*Clean up and return*/
5862 delete gauss;
5863 return pe;
5864}
5865/*}}}*/
5866#endif
Note: See TracBrowser for help on using the repository browser.