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

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

stripped inversion threshold out of code, not used
File size: 183.9 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 bool control_analysis=false;
1380 int num_control_type;
1381 double yts;
1382 int num_cm_responses;
1383
1384 /*Get parameters: */
1385 iomodel->Constant(&yts,ConstantsYtsEnum);
1386 iomodel->Constant(&control_analysis,InversionIscontrolEnum);
1387 if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
1388 if(control_analysis) iomodel->Constant(&num_cm_responses,InversionNumCostFunctionsEnum);
1389
1390 /*Recover vertices ids needed to initialize inputs*/
1391 for(i=0;i<3;i++){
1392 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
1393 }
1394
1395 /*Control Inputs*/
1396 #ifdef _HAVE_CONTROL_
1397 if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
1398 for(i=0;i<num_control_type;i++){
1399 switch((int)iomodel->Data(InversionControlParametersEnum)[i]){
1400 case BalancethicknessThickeningRateEnum:
1401 if (iomodel->Data(BalancethicknessThickeningRateEnum)){
1402 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[tria_vertex_ids[j]-1]/yts;
1403 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1404 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1405 this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
1406 }
1407 break;
1408 case VxEnum:
1409 if (iomodel->Data(VxEnum)){
1410 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VxEnum)[tria_vertex_ids[j]-1]/yts;
1411 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1412 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1413 this->inputs->AddInput(new ControlInput(VxEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
1414 }
1415 break;
1416 case VyEnum:
1417 if (iomodel->Data(VyEnum)){
1418 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VyEnum)[tria_vertex_ids[j]-1]/yts;
1419 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1420 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
1421 this->inputs->AddInput(new ControlInput(VyEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
1422 }
1423 break;
1424 case FrictionCoefficientEnum:
1425 if (iomodel->Data(FrictionCoefficientEnum)){
1426 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[tria_vertex_ids[j]-1];
1427 for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
1428 for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
1429 this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
1430 }
1431 break;
1432 case MaterialsRheologyBbarEnum:
1433 case MaterialsRheologyZbarEnum:
1434 /*Matice will take care of it*/ break;
1435 default:
1436 _error_("Control %s not implemented yet",EnumToStringx((int)iomodel->Data(InversionControlParametersEnum)[i]));
1437 }
1438 }
1439 }
1440 #endif
1441
1442 /*DatasetInputs*/
1443 if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)){
1444
1445 /*Create inputs and add to DataSetInput*/
1446 DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
1447 for(i=0;i<num_cm_responses;i++){
1448 for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(tria_vertex_ids[j]-1)*num_cm_responses+i];
1449 datasetinput->inputs->AddObject(new TriaP1Input(InversionCostFunctionsCoefficientsEnum,nodeinputs));
1450 }
1451
1452 /*Add datasetinput to element inputs*/
1453 this->inputs->AddInput(datasetinput);
1454 }
1455}
1456/*}}}*/
1457/*FUNCTION Tria::InputUpdateFromSolution {{{1*/
1458void Tria::InputUpdateFromSolution(double* solution){
1459
1460 /*retrive parameters: */
1461 int analysis_type;
1462 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
1463
1464 /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
1465 switch(analysis_type){
1466 #ifdef _HAVE_DIAGNOSTIC_
1467 case DiagnosticHorizAnalysisEnum:
1468 InputUpdateFromSolutionDiagnosticHoriz( solution);
1469 break;
1470 case DiagnosticHutterAnalysisEnum:
1471 InputUpdateFromSolutionDiagnosticHoriz( solution);
1472 break;
1473 #endif
1474 #ifdef _HAVE_CONTROL_
1475 case AdjointHorizAnalysisEnum:
1476 InputUpdateFromSolutionAdjointHoriz( solution);
1477 break;
1478 case AdjointBalancethicknessAnalysisEnum:
1479 InputUpdateFromSolutionAdjointBalancethickness( solution);
1480 break;
1481 #endif
1482 #ifdef _HAVE_HYDROLOGY_
1483 case HydrologyAnalysisEnum:
1484 InputUpdateFromSolutionHydrology(solution);
1485 break ;
1486 #endif
1487 #ifdef _HAVE_BALANCED_
1488 case BalancethicknessAnalysisEnum:
1489 InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
1490 break;
1491 #endif
1492 case BedSlopeXAnalysisEnum:
1493 InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum);
1494 break;
1495 case BedSlopeYAnalysisEnum:
1496 InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum);
1497 break;
1498 case SurfaceSlopeXAnalysisEnum:
1499 InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
1500 break;
1501 case SurfaceSlopeYAnalysisEnum:
1502 InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
1503 break;
1504 case PrognosticAnalysisEnum:
1505 InputUpdateFromSolutionPrognostic(solution);
1506 break;
1507 default:
1508 _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
1509 }
1510}
1511/*}}}*/
1512/*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{1*/
1513void Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){
1514
1515 const int numdof = NDOF1*NUMVERTICES;
1516
1517 int* doflist=NULL;
1518 double values[numdof];
1519
1520 /*Get dof list: */
1521 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
1522
1523 /*Use the dof list to index into the solution vector: */
1524 for(int i=0;i<numdof;i++){
1525 values[i]=solution[doflist[i]];
1526 if(isnan(values[i])) _error_("NaN found in solution vector");
1527 }
1528
1529 /*Add input to the element: */
1530 this->inputs->AddInput(new TriaP1Input(enum_type,values));
1531
1532 /*Free ressources:*/
1533 xfree((void**)&doflist);
1534}
1535/*}}}*/
1536/*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{1*/
1537void Tria::InputUpdateFromSolutionPrognostic(double* solution){
1538
1539 /*Intermediaries*/
1540 const int numdof = NDOF1*NUMVERTICES;
1541
1542 int i,hydroadjustment;
1543 int* doflist=NULL;
1544 double rho_ice,rho_water,minthickness;
1545 double newthickness[numdof];
1546 double newbed[numdof];
1547 double newsurface[numdof];
1548 double oldbed[NUMVERTICES];
1549 double oldsurface[NUMVERTICES];
1550 double oldthickness[NUMVERTICES];
1551
1552 /*Get dof list: */
1553 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
1554
1555 /*Use the dof list to index into the solution vector: */
1556 this->parameters->FindParam(&minthickness,PrognosticMinThicknessEnum);
1557 for(i=0;i<numdof;i++){
1558 newthickness[i]=solution[doflist[i]];
1559 if(isnan(newthickness[i])) _error_("NaN found in solution vector");
1560 /*Constrain thickness to be at least 1m*/
1561 if(newthickness[i]<minthickness) newthickness[i]=minthickness;
1562 }
1563
1564 /*Get previous bed, thickness and surface*/
1565 GetInputListOnVertices(&oldbed[0],BedEnum);
1566 GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
1567 GetInputListOnVertices(&oldthickness[0],ThicknessEnum);
1568
1569 /*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
1570 this->parameters->FindParam(&hydroadjustment,PrognosticHydrostaticAdjustmentEnum);
1571 rho_ice=matpar->GetRhoIce();
1572 rho_water=matpar->GetRhoWater();
1573
1574 for(i=0;i<numdof;i++) {
1575 /*If shelf: hydrostatic equilibrium*/
1576 if (this->nodes[i]->IsGrounded()){
1577 newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
1578 newbed[i]=oldbed[i]; //same bed: do nothing
1579 }
1580 else{ //this is an ice shelf
1581
1582 if(hydroadjustment==AbsoluteEnum){
1583 newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
1584 newbed[i]=newthickness[i]*(-rho_ice/rho_water);
1585 }
1586 else if(hydroadjustment==IncrementalEnum){
1587 newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH
1588 newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
1589 }
1590 else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment));
1591 }
1592 }
1593
1594 /*Add input to the element: */
1595 this->inputs->AddInput(new TriaP1Input(ThicknessEnum,newthickness));
1596 this->inputs->AddInput(new TriaP1Input(SurfaceEnum,newsurface));
1597 this->inputs->AddInput(new TriaP1Input(BedEnum,newbed));
1598
1599 /*Free ressources:*/
1600 xfree((void**)&doflist);
1601}
1602/*}}}*/
1603/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
1604void Tria::InputUpdateFromVector(double* vector, int name, int type){
1605
1606 /*Check that name is an element input*/
1607 if (!IsInput(name)) return;
1608
1609 switch(type){
1610
1611 case VertexEnum:
1612
1613 /*New TriaP1Input*/
1614 double values[3];
1615
1616 /*Get values on the 3 vertices*/
1617 for (int i=0;i<3;i++){
1618 values[i]=vector[this->nodes[i]->GetVertexDof()];
1619 }
1620
1621 /*update input*/
1622 if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
1623 matice->inputs->AddInput(new TriaP1Input(name,values));
1624 }
1625 else{
1626 this->inputs->AddInput(new TriaP1Input(name,values));
1627 }
1628 return;
1629
1630 default:
1631 _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
1632 }
1633}
1634/*}}}*/
1635/*FUNCTION Tria::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
1636void Tria::InputUpdateFromVector(int* vector, int name, int type){
1637 _error_(" not supported yet!");
1638}
1639/*}}}*/
1640/*FUNCTION Tria::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
1641void Tria::InputUpdateFromVector(bool* vector, int name, int type){
1642 _error_(" not supported yet!");
1643}
1644/*}}}*/
1645/*FUNCTION Tria::InputCreate(double scalar,int enum,int code);{{{1*/
1646void Tria::InputCreate(double scalar,int name,int code){
1647
1648 /*Check that name is an element input*/
1649 if (!IsInput(name)) return;
1650
1651 if ((code==5) || (code==1)){ //boolean
1652 this->inputs->AddInput(new BoolInput(name,(bool)scalar));
1653 }
1654 else if ((code==6) || (code==2)){ //integer
1655 this->inputs->AddInput(new IntInput(name,(int)scalar));
1656 }
1657 else if ((code==7) || (code==3)){ //double
1658 this->inputs->AddInput(new DoubleInput(name,(int)scalar));
1659 }
1660 else _error_("%s%i"," could not recognize nature of vector from code ",code);
1661
1662}
1663/*}}}*/
1664/*FUNCTION Tria::InputCreate(double* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{1*/
1665void Tria::InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements
1666
1667 /*Intermediaries*/
1668 int i,j,t;
1669 int tria_vertex_ids[3];
1670 int row;
1671 double nodeinputs[3];
1672 double time;
1673 TransientInput* transientinput=NULL;
1674 int numberofvertices;
1675 int numberofelements;
1676 double yts;
1677
1678
1679 /*Fetch parameters: */
1680 iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
1681 iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
1682 iomodel->Constant(&yts,ConstantsYtsEnum);
1683
1684 /*Branch on type of vector: nodal or elementary: */
1685 if(vector_type==1){ //nodal vector
1686
1687 /*Recover vertices ids needed to initialize inputs*/
1688 for(i=0;i<3;i++){
1689 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
1690 }
1691
1692 /*Are we in transient or static? */
1693 if(M==numberofvertices){
1694
1695 /*create input values: */
1696 for(i=0;i<3;i++)nodeinputs[i]=(double)vector[tria_vertex_ids[i]-1];
1697
1698 /*process units: */
1699 UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
1700
1701 /*create static input: */
1702 this->inputs->AddInput(new TriaP1Input(vector_enum,nodeinputs));
1703 }
1704 else if(M==numberofvertices+1){
1705 /*create transient input: */
1706 for(t=0;t<N;t++){ //N is the number of times
1707
1708 /*create input values: */
1709 for(i=0;i<3;i++){
1710 row=tria_vertex_ids[i]-1;
1711 nodeinputs[i]=(double)vector[N*row+t];
1712 }
1713
1714 /*process units: */
1715 UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);
1716
1717 /*time? :*/
1718 time=(double)vector[(M-1)*N+t]*yts;
1719
1720 if(t==0) transientinput=new TransientInput(vector_enum);
1721 transientinput->AddTimeInput(new TriaP1Input(vector_enum,nodeinputs),time);
1722 }
1723 this->inputs->AddInput(transientinput);
1724 }
1725 else _error_("nodal vector is either numberofnodes or numberofnodes+1 long. Field provided (%s) is %i long",EnumToStringx(vector_enum),M);
1726 }
1727 else if(vector_type==2){ //element vector
1728 /*Are we in transient or static? */
1729 if(M==numberofelements){
1730
1731 /*static mode: create an input out of the element value: */
1732
1733 if (code==5){ //boolean
1734 this->inputs->AddInput(new BoolInput(vector_enum,(bool)vector[index]));
1735 }
1736 else if (code==6){ //integer
1737 this->inputs->AddInput(new IntInput(vector_enum,(int)vector[index]));
1738 }
1739 else if (code==7){ //double
1740 this->inputs->AddInput(new DoubleInput(vector_enum,(double)vector[index]));
1741 }
1742 else _error_("%s%i"," could not recognize nature of vector from code ",code);
1743 }
1744 else {
1745 _error_("transient elementary inputs not supported yet!");
1746 }
1747 }
1748 else{
1749 _error_("Cannot add input for vector type %i (not supported)",vector_type);
1750 }
1751
1752}
1753/*}}}*/
1754/*FUNCTION Tria::IsInput{{{1*/
1755bool Tria::IsInput(int name){
1756 if (
1757 name==ThicknessEnum ||
1758 name==SurfaceEnum ||
1759 name==BedEnum ||
1760 name==SurfaceSlopeXEnum ||
1761 name==SurfaceSlopeYEnum ||
1762 name==BasalforcingsMeltingRateEnum ||
1763 name==WatercolumnEnum ||
1764 name==SurfaceforcingsMassBalanceEnum ||
1765 name==SurfaceAreaEnum||
1766 name==VxEnum ||
1767 name==VyEnum ||
1768 name==InversionVxObsEnum ||
1769 name==InversionVyObsEnum ||
1770 name==FrictionCoefficientEnum ||
1771 name==MaterialsRheologyBbarEnum ||
1772 name==GradientEnum ||
1773 name==OldGradientEnum ||
1774 name==QmuVxEnum ||
1775 name==QmuVyEnum ||
1776 name==QmuPressureEnum ||
1777 name==QmuBedEnum ||
1778 name==QmuThicknessEnum ||
1779 name==QmuSurfaceEnum ||
1780 name==QmuTemperatureEnum ||
1781 name==QmuMeltingEnum
1782 ){
1783 return true;
1784 }
1785 else return false;
1786}
1787/*}}}*/
1788/*FUNCTION Tria::IsOnBed {{{1*/
1789bool Tria::IsOnBed(){
1790
1791 bool onbed;
1792 inputs->GetInputValue(&onbed,MeshElementonbedEnum);
1793 return onbed;
1794}
1795/*}}}*/
1796/*FUNCTION Tria::IsFloating {{{1*/
1797bool Tria::IsFloating(){
1798
1799 bool shelf;
1800 inputs->GetInputValue(&shelf,MaskElementonfloatingiceEnum);
1801 return shelf;
1802}
1803/*}}}*/
1804/*FUNCTION Tria::IsNodeOnShelf {{{1*/
1805bool Tria::IsNodeOnShelf(){
1806
1807 int i;
1808 bool shelf=false;
1809
1810 for(i=0;i<3;i++){
1811 if (nodes[i]->IsFloating()){
1812 shelf=true;
1813 break;
1814 }
1815 }
1816 return shelf;
1817}
1818/*}}}*/
1819/*FUNCTION Tria::IsNodeOnShelfFromFlags {{{1*/
1820bool Tria::IsNodeOnShelfFromFlags(double* flags){
1821
1822 int i;
1823 bool shelf=false;
1824
1825 for(i=0;i<NUMVERTICES;i++){
1826 if (flags[nodes[i]->Sid()]){
1827 shelf=true;
1828 break;
1829 }
1830 }
1831 return shelf;
1832}
1833/*}}}*/
1834/*FUNCTION Tria::IsOnWater {{{1*/
1835bool Tria::IsOnWater(){
1836
1837 bool water;
1838 inputs->GetInputValue(&water,MaskElementonwaterEnum);
1839 return water;
1840}
1841/*}}}*/
1842/*FUNCTION Tria::ListResultsInfo{{{*/
1843void Tria::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,double** in_resultstimes,int** in_resultssteps,int* in_num_results){
1844
1845 /*Intermediaries*/
1846 int i;
1847 int numberofresults = 0;
1848 int *resultsenums = NULL;
1849 int *resultssizes = NULL;
1850 double *resultstimes = NULL;
1851 int *resultssteps = NULL;
1852
1853 /*Checks*/
1854 _assert_(in_num_results);
1855
1856 /*Count number of results*/
1857 for(i=0;i<this->results->Size();i++){
1858 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
1859 numberofresults++;
1860 }
1861
1862 if(numberofresults){
1863
1864 /*Allocate output*/
1865 resultsenums=(int*)xmalloc(numberofresults*sizeof(int));
1866 resultssizes=(int*)xmalloc(numberofresults*sizeof(int));
1867 resultstimes=(double*)xmalloc(numberofresults*sizeof(double));
1868 resultssteps=(int*)xmalloc(numberofresults*sizeof(int));
1869
1870 /*populate enums*/
1871 for(i=0;i<this->results->Size();i++){
1872 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
1873 resultsenums[i]=elementresult->InstanceEnum();
1874 resultstimes[i]=elementresult->GetTime();
1875 resultssteps[i]=elementresult->GetStep();
1876 if(elementresult->ObjectEnum()==TriaP1ElementResultEnum){
1877 resultssizes[i]=P1Enum;
1878 }
1879 else{
1880 resultssizes[i]=P0Enum;
1881 }
1882 }
1883 }
1884
1885 /*Assign output pointers:*/
1886 *in_num_results=numberofresults;
1887 *in_resultsenums=resultsenums;
1888 *in_resultssizes=resultssizes;
1889 *in_resultstimes=resultstimes;
1890 *in_resultssteps=resultssteps;
1891
1892}/*}}}*/
1893/*FUNCTION Tria::MigrateGroundingLine{{{1*/
1894void Tria::MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding){
1895
1896 int i,migration_style,unground;
1897 bool elementonshelf = false;
1898 double bed_hydro,yts,gl_melting_rate;
1899 double rho_water,rho_ice,density;
1900 double melting[NUMVERTICES];
1901 double h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];
1902
1903 /*Recover info at the vertices: */
1904 parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
1905 parameters->FindParam(&yts,ConstantsYtsEnum);
1906 GetInputListOnVertices(&h[0],ThicknessEnum);
1907 GetInputListOnVertices(&s[0],SurfaceEnum);
1908 GetInputListOnVertices(&b[0],BedEnum);
1909 GetInputListOnVertices(&ba[0],BathymetryEnum);
1910 rho_water=matpar->GetRhoWater();
1911 rho_ice=matpar->GetRhoIce();
1912 density=rho_ice/rho_water;
1913
1914 /*go through vertices, and update inputs, considering them to be TriaVertex type: */
1915 for(i=0;i<NUMVERTICES;i++){
1916 /*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
1917 if(old_floating_ice[nodes[i]->Sid()]){
1918 if(b[i]<=ba[i]){
1919 b[i]=ba[i];
1920 s[i]=b[i]+h[i];
1921 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
1922 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
1923 }
1924 }
1925 /*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
1926 /*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
1927 else{
1928 bed_hydro=-density*h[i];
1929 if (bed_hydro>ba[i]){
1930 /*Unground only if the element is connected to the ice shelf*/
1931 if(migration_style==AgressiveMigrationEnum){
1932 s[i]=(1-density)*h[i];
1933 b[i]=-density*h[i];
1934 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
1935 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
1936 }
1937 else if(migration_style==SoftMigrationEnum && sheet_ungrounding[nodes[i]->Sid()]){
1938 s[i]=(1-density)*h[i];
1939 b[i]=-density*h[i];
1940 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
1941 nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
1942 }
1943 }
1944 }
1945 }
1946
1947 /*If at least one vertex is now floating, the element is now floating*/
1948 for(i=0;i<NUMVERTICES;i++){
1949 if(nodes[i]->IsFloating()){
1950 elementonshelf=true;
1951 break;
1952 }
1953 }
1954
1955 /*Add basal melting rate if element just ungrounded*/
1956 if(!this->IsFloating() && elementonshelf==true){
1957 for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
1958 this->inputs->AddInput(new TriaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
1959 }
1960
1961 /*Update inputs*/
1962 this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));
1963
1964 /*Update inputs*/
1965 this->inputs->AddInput(new TriaP1Input(SurfaceEnum,&s[0]));
1966 this->inputs->AddInput(new TriaP1Input(BedEnum,&b[0]));
1967}
1968/*}}}*/
1969/*FUNCTION Tria::MyRank {{{1*/
1970int Tria::MyRank(void){
1971 extern int my_rank;
1972 return my_rank;
1973}
1974/*}}}*/
1975/*FUNCTION Tria::NodalValue {{{1*/
1976int Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){
1977
1978 int i;
1979 int found=0;
1980 double value;
1981 Input* data=NULL;
1982 GaussTria *gauss = NULL;
1983
1984 /*First, serarch the input: */
1985 data=inputs->GetInput(natureofdataenum);
1986
1987 /*figure out if we have the vertex id: */
1988 found=0;
1989 for(i=0;i<NUMVERTICES;i++){
1990 if(index==nodes[i]->GetVertexId()){
1991 /*Do we have natureofdataenum in our inputs? :*/
1992 if(data){
1993 /*ok, we are good. retrieve value of input at vertex :*/
1994 gauss=new GaussTria(); gauss->GaussVertex(i);
1995 data->GetInputValue(&value,gauss);
1996 found=1;
1997 break;
1998 }
1999 }
2000 }
2001
2002 if(found)*pvalue=value;
2003 return found;
2004}
2005/*}}}*/
2006/*FUNCTION Tria::PatchFill{{{1*/
2007void Tria::PatchFill(int* prow, Patch* patch){
2008
2009 int i,row;
2010 int vertices_ids[3];
2011
2012 /*recover pointer: */
2013 row=*prow;
2014
2015 for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.
2016
2017 for(i=0;i<this->results->Size();i++){
2018 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
2019
2020 /*For this result,fill the information in the Patch object (element id + vertices ids), and then hand
2021 *it to the result object, to fill the rest: */
2022 patch->fillelementinfo(row,this->sid+1,vertices_ids,3);
2023 elementresult->PatchFill(row,patch);
2024
2025 /*increment rower: */
2026 row++;
2027 }
2028
2029 /*Assign output pointers:*/
2030 *prow=row;
2031}
2032/*}}}*/
2033/*FUNCTION Tria::PatchSize{{{1*/
2034void Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){
2035
2036 int i;
2037 int numrows = 0;
2038 int numnodes = 0;
2039 int temp_numnodes = 0;
2040
2041 /*Go through all the results objects, and update the counters: */
2042 for (i=0;i<this->results->Size();i++){
2043 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
2044 /*first, we have one more result: */
2045 numrows++;
2046 /*now, how many vertices and how many nodal values for this result? :*/
2047 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
2048 if(temp_numnodes>numnodes)numnodes=temp_numnodes;
2049 }
2050
2051 /*Assign output pointers:*/
2052 *pnumrows=numrows;
2053 *pnumvertices=NUMVERTICES;
2054 *pnumnodes=numnodes;
2055}
2056/*}}}*/
2057/*FUNCTION Tria::PotentialSheetUngrounding{{{1*/
2058void Tria::PotentialSheetUngrounding(Vector* potential_sheet_ungrounding){
2059
2060 int i;
2061 double h[NUMVERTICES],ba[NUMVERTICES];
2062 double bed_hydro;
2063 double rho_water,rho_ice,density;
2064 bool elementonshelf = false;
2065
2066 /*material parameters: */
2067 rho_water=matpar->GetRhoWater();
2068 rho_ice=matpar->GetRhoIce();
2069 density=rho_ice/rho_water;
2070 GetInputListOnVertices(&h[0],ThicknessEnum);
2071 GetInputListOnVertices(&ba[0],BathymetryEnum);
2072
2073 /*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
2074 for(i=0;i<NUMVERTICES;i++){
2075 /*Find if grounded vertices want to start floating*/
2076 if (!nodes[i]->IsFloating()){
2077 bed_hydro=-density*h[i];
2078 if (bed_hydro>ba[i]){
2079 /*Vertex that could potentially unground, flag it*/
2080 potential_sheet_ungrounding->SetValue(nodes[i]->Sid(),1,INS_VAL);
2081 }
2082 }
2083 }
2084}
2085/*}}}*/
2086/*FUNCTION Tria::PositiveDegreeDay{{{1*/
2087void Tria::PositiveDegreeDay(double* pdds,double* pds,double signorm){
2088
2089 int i,iqj,imonth;
2090 double agd[NUMVERTICES]; // surface and basal
2091 double saccu[NUMVERTICES] = {0}; // yearly surface accumulation
2092 double smelt[NUMVERTICES] = {0}; // yearly melt
2093 double precrunoff[NUMVERTICES]; // yearly runoff
2094 double prect; // total precipitation during 1 year taking into account des. ef.
2095 double water; //water=rain + snowmelt
2096 double runoff; //meltwater only, does not include rain
2097 double sconv; //rhow_rain/rhoi / 12 months
2098
2099 double rho_water,rho_ice,density;
2100 double lapser=6.5/1000, sealev=0; // lapse rate. degrees per meter.
2101 double desfac = 0.5; //desert elevation factor
2102 double s0p[NUMVERTICES]={0}; //should be set to elevation from precip source
2103 double s0t[NUMVERTICES]={0}; //should be set to elevation from temperature source
2104 double st; // elevation between altitude of the temp record and current altitude
2105 double sp; // elevation between altitude of the prec record and current altitude
2106
2107
2108 // PDD and PD constants and variables
2109 double siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
2110 double siglimc=0, siglim0, siglim0c;
2111 double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
2112 double DT = 0.02;
2113 double pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
2114
2115 double q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
2116 double qm[NUMVERTICES] = {0}; // snow part of the precipitation
2117 double qmt[NUMVERTICES] = {0}; // precipitation without desertification effect adjustment
2118 double qmp[NUMVERTICES] = {0}; // desertification taken into account
2119 double pdd[NUMVERTICES] = {0};
2120 double frzndd[NUMVERTICES] = {0};
2121
2122 double tstar; // monthly mean surface temp
2123 double Tsum[NUMVERTICES]= {0}; // average summer (JJA) temperature
2124 double Tsurf[NUMVERTICES] = {0}; // average annual temperature
2125
2126 double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
2127 double t[NUMVERTICES][12],prec[NUMVERTICES][12];
2128 double deltm=1/12;
2129 int ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
2130
2131 double snwm; // snow that could have been melted in a year.
2132 double snwmf; // ablation factor for snow per positive degree day.
2133 double smf; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
2134
2135 double dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
2136 double supice,supcap,diffndd;
2137 double fsupT=0.5, fsupndd=0.5; // Tsurf mode factors for supice
2138 double pddtj[NUMVERTICES], hmx2;
2139
2140 /*Recover info at the vertices: */
2141 GetInputListOnVertices(&h[0],ThicknessEnum);
2142 GetInputListOnVertices(&s[0],SurfaceEnum);
2143 GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum);
2144 GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
2145
2146 for(i=0;i<NUMVERTICES;i++) ttmp[i]=ttmp[i]-273.15; // convertion from Kelvin to celcius
2147
2148 for(i=0;i<NUMVERTICES;i++)
2149 for(imonth=0;imonth<12;imonth++){
2150 t[i][imonth]=ttmp[i];
2151 prec[i][imonth]=prectmp[i];
2152 }
2153
2154 /*Get material parameters :*/
2155 rho_ice=matpar->GetRhoIce();
2156 rho_water=matpar->GetRhoWater();
2157 density=rho_ice/rho_water;
2158
2159 sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
2160
2161 /*PDD constant*/
2162 siglim = 2.5*signorm;
2163 siglim0 = siglim/DT + 0.5;
2164 siglim0c = siglimc/DT + 0.5;
2165 PDup = siglimc+PDCUT;
2166
2167 // seasonal loop
2168 for (iqj = 0; iqj < 12; iqj++){
2169 imonth = ismon[iqj];
2170 for (i = 0; i < NUMVERTICES; i++){
2171 st=(s[i]-s0t[i])/1000;
2172 tstar = t[i][imonth] - lapser *max(st,sealev);
2173 Tsurf[i] = tstar*deltm+Tsurf[i];
2174
2175 /*********compute PD ****************/
2176 if (tstar < PDup){
2177 pd = 1;
2178 if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
2179 else {
2180 pd = 0;}
2181
2182 /******exp des/elev precip reduction*******/
2183 sp=(s[i]-s0p[i])/1000; // deselev effect is wrt chng in topo
2184 if (sp>0.0){q = exp(-desfac*sp);}
2185 else {q = 1.0;}
2186
2187 qmt[i]= qmt[i] + prec[i][imonth]*sconv; //*sconv to convert in m of ice equivalent
2188 qmpt= q*prec[i][imonth]*sconv;
2189 qmp[i]= qmp[i] + qmpt;
2190 qm[i]= qm[i] + qmpt*pd;
2191
2192 /*********compute PDD************/
2193 // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
2194 // gaussian=T_m, so ndd=-(Tsurf-pdd)
2195 if (iqj>6 && iqj<10){ Tsum[i]=Tsum[i]+tstar;}
2196 if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
2197 else if (tstar> -siglim){
2198 pddsig=pdds[int(tstar/DT + siglim0)];
2199 pdd[i] = pdd[i] + pddsig*deltm;
2200 frzndd[i] = frzndd[i] - (tstar-pddsig)*deltm;}
2201 else{frzndd[i] = frzndd[i] - tstar*deltm; }
2202 }
2203 } // end of seasonal loop
2204
2205 //******************************************************************
2206 for(i=0;i<NUMVERTICES;i++){
2207 saccu[i] = qm[i];
2208 prect = qmp[i]; // total precipitation during 1 year taking into account des. ef.
2209 Tsum[i]=Tsum[i]/3;
2210
2211 /***** determine PDD factors *****/
2212 if(Tsum[i]< -1.) {
2213 snwmf=2.65*0.001; // ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
2214 smf=17.22*0.001; // ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
2215 }
2216 else if(Tsum[i]< 10){
2217 snwmf = (0.15*Tsum[i] + 2.8)*0.001;
2218 smf = (0.0067*pow((10.-Tsum[i]),3) + 8.3)*0.001;
2219 }
2220 else{
2221 snwmf=4.3*0.001;
2222 smf=8.3*0.001;
2223 }
2224 snwmf=0.95*snwmf;
2225 smf=0.95*smf;
2226
2227 /***** compute PDD ablation and refreezing *****/
2228 pddt = pdd[i] *365;
2229 snwm = snwmf*pddt; // snow that could have been melted in a year
2230 hmx2 = min(h[i],dfrz); // refreeze active layer max depth: dfrz
2231
2232 if(snwm < saccu[i]) {
2233 water=prect-saccu[i] + snwm; //water=rain + snowmelt
2234 // l 2.2= capillary factor
2235 // Should refreezing be controlled by frzndd or by mean annual Tsurf?
2236 // dCovrLm concept is of warming of active layer (thickness =d@=1-
2237 // >2m)
2238 // problem with water seepage into ice: should be sealed after
2239 // refreezing
2240 // so everything needs to be predicated on 1 year scale, except for
2241 // thermal
2242 // conductivity through ice
2243 // also, need to account that melt season has low accum, so what's
2244 // going to
2245 // hold the meltwater around for refreezing? And melt-time will have
2246 // low seasonal frzndd
2247
2248 // Superimposed ice : Pfeffer et al. 1991, Tarasov 2002
2249
2250 supice= min(hmx2*CovrLm*frzndd[i]+2.2*(saccu[i]-snwm), water); // superimposed ice
2251 supcap=min(2.2*(saccu[i]-snwm),water);
2252 runoff=snwm - supice; //meltwater only, does not include rain
2253 }
2254 else { //all snow melted
2255 supice= min(hmx2*CovrLm*frzndd[i], prect );
2256 runoff= saccu[i] + smf*(pddt-saccu[i]/snwmf) - supice;
2257 supcap=0;
2258 }
2259 // pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
2260 // except pdd melt heat source is atmosphere, while refreeze is
2261 // ground/ice stored interim
2262 // assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
2263 // assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
2264 // <0
2265 // assume ndd>pdd, little melt => little supice
2266 // bottom line: compare for Tsurf<0 : supice and no supice case,
2267 // expect Tsurf difference
2268 // except some of cooling flux comes from atmosphere//
2269 // 1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
2270 // does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
2271 // < 0.1
2272
2273 // make more sense to just use residual pdd-ndd except that pdd
2274 // residual not clear yet
2275 // frzndd should not be used up by refreezing in snow, so stick in
2276 // supcap.
2277 diffndd=0;
2278 if (frzndd[i]>0) {
2279 diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd[i]);
2280 frzndd[i]=frzndd[i]-diffndd;
2281 }
2282 if(runoff<0){
2283 saccu[i]= saccu[i] -runoff;
2284 smelt[i] = 0;
2285 precrunoff[i]=prect-saccu[i];
2286 //here assume pdd residual is 0, =>
2287 Tsurf[i]= max(Tsurf[i],-frzndd[i]);
2288 }
2289 else {
2290 smelt[i] = runoff;
2291 precrunoff[i]=prect-max(0.,supice)-saccu[i];}
2292 //here really need pdd balance, try 0.5 fudge factor?
2293 //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
2294 //yet from site plots, can be ice free with Tsurf=-5.5C
2295 if(Tsurf[i]<0) {
2296 Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
2297
2298 agd[i] = -smelt[i]+saccu[i];
2299 pddtj[i]=pddt;
2300
2301 /*Update inputs*/
2302 this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
2303 // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
2304 this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
2305
2306 } //end of the for loop over the vertices
2307}
2308/*}}}*/
2309/*FUNCTION Tria::ProcessResultsUnits{{{1*/
2310void Tria::ProcessResultsUnits(void){
2311
2312 int i;
2313
2314 for(i=0;i<this->results->Size();i++){
2315 ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
2316 elementresult->ProcessUnits(this->parameters);
2317 }
2318}
2319/*}}}*/
2320/*FUNCTION Tria::RequestedOutput{{{1*/
2321void Tria::RequestedOutput(int output_enum,int step,double time){
2322
2323 if(IsInput(output_enum)){
2324 /*just transfer this input to results, and we are done: */
2325 InputToResult(output_enum,step,time);
2326 }
2327 else{
2328 /*this input does not exist, compute it, and then transfer to results: */
2329 switch(output_enum){
2330 case StressTensorEnum:
2331 this->ComputeStressTensor();
2332 InputToResult(StressTensorxxEnum,step,time);
2333 InputToResult(StressTensorxyEnum,step,time);
2334 InputToResult(StressTensorxzEnum,step,time);
2335 InputToResult(StressTensoryyEnum,step,time);
2336 InputToResult(StressTensoryzEnum,step,time);
2337 InputToResult(StressTensorzzEnum,step,time);
2338 break;
2339
2340 default:
2341 /*do nothing, no need to derail the computation because one of the outputs requested cannot be found: */
2342 break;
2343 }
2344 }
2345
2346}
2347/*}}}*/
2348/*FUNCTION Tria::SetClone {{{1*/
2349void Tria::SetClone(int* minranks){
2350
2351 _error_("not implemented yet");
2352}
2353/*}}}1*/
2354/*FUNCTION Tria::SmearFunction {{{1*/
2355void Tria::SmearFunction(Vector* smearedvector,double (*WeightFunction)(double distance,double radius),double radius){
2356 _error_("not implemented yet");
2357
2358}
2359/*}}}1*/
2360/*FUNCTION Tria::SetCurrentConfiguration {{{1*/
2361void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
2362
2363 /*go into parameters and get the analysis_counter: */
2364 int analysis_counter;
2365 parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
2366
2367 /*Get Element type*/
2368 this->element_type=this->element_type_list[analysis_counter];
2369
2370 /*Pick up nodes*/
2371 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
2372 else this->nodes=NULL;
2373
2374}
2375/*}}}*/
2376/*FUNCTION Tria::SurfaceArea {{{1*/
2377double Tria::SurfaceArea(void){
2378
2379 int i;
2380 double S;
2381 double normal[3];
2382 double v13[3],v23[3];
2383 double xyz_list[NUMVERTICES][3];
2384
2385 /*If on water, return 0: */
2386 if(IsOnWater())return 0;
2387
2388 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2389
2390 for (i=0;i<3;i++){
2391 v13[i]=xyz_list[0][i]-xyz_list[2][i];
2392 v23[i]=xyz_list[1][i]-xyz_list[2][i];
2393 }
2394
2395 normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
2396 normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
2397 normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
2398
2399 S = 0.5 * sqrt(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));
2400
2401 /*Return: */
2402 return S;
2403}
2404/*}}}*/
2405/*FUNCTION Tria::SurfaceNormal{{{1*/
2406void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
2407
2408 int i;
2409 double v13[3],v23[3];
2410 double normal[3];
2411 double normal_norm;
2412
2413 for (i=0;i<3;i++){
2414 v13[i]=xyz_list[0][i]-xyz_list[2][i];
2415 v23[i]=xyz_list[1][i]-xyz_list[2][i];
2416 }
2417
2418 normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
2419 normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
2420 normal[2]=v13[0]*v23[1]-v13[1]*v23[0];
2421
2422 normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );
2423
2424 *(surface_normal)=normal[0]/normal_norm;
2425 *(surface_normal+1)=normal[1]/normal_norm;
2426 *(surface_normal+2)=normal[2]/normal_norm;
2427}
2428/*}}}*/
2429/*FUNCTION Tria::TimeAdapt{{{1*/
2430double Tria::TimeAdapt(void){
2431
2432 /*intermediary: */
2433 int i;
2434 double C,dt;
2435 double dx,dy;
2436 double maxx,minx;
2437 double maxy,miny;
2438 double maxabsvx,maxabsvy;
2439 double xyz_list[NUMVERTICES][3];
2440
2441 /*get CFL coefficient:*/
2442 this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum);
2443
2444 /*Get for Vx and Vy, the max of abs value: */
2445 #ifdef _HAVE_RESPONSES_
2446 this->MaxAbsVx(&maxabsvx,false);
2447 this->MaxAbsVy(&maxabsvy,false);
2448 #else
2449 _error_("ISSM was not compiled with responses compiled in, exiting!");
2450 #endif
2451
2452 /* Get node coordinates and dof list: */
2453 GetVerticesCoordinates(&xyz_list[0][0], this->nodes, NUMVERTICES);
2454
2455 minx=xyz_list[0][0];
2456 maxx=xyz_list[0][0];
2457 miny=xyz_list[0][1];
2458 maxy=xyz_list[0][1];
2459
2460 for(i=1;i<NUMVERTICES;i++){
2461 if (xyz_list[i][0]<minx)minx=xyz_list[i][0];
2462 if (xyz_list[i][0]>maxx)maxx=xyz_list[i][0];
2463 if (xyz_list[i][1]<miny)miny=xyz_list[i][1];
2464 if (xyz_list[i][1]>maxy)maxy=xyz_list[i][1];
2465 }
2466 dx=maxx-minx;
2467 dy=maxy-miny;
2468
2469 /*CFL criterion: */
2470 dt=C/(maxabsvy/dx+maxabsvy/dy);
2471
2472 return dt;
2473}
2474/*}}}*/
2475/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/
2476void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index
2477
2478 /*Intermediaries*/
2479 int i,j;
2480 int tria_node_ids[3];
2481 int tria_vertex_ids[3];
2482 int tria_type;
2483 double nodeinputs[3];
2484 double yts;
2485 int progstabilization,balancestabilization;
2486 bool dakota_analysis;
2487
2488 /*Checks if debuging*/
2489 /*{{{2*/
2490 _assert_(iomodel->Data(MeshElementsEnum));
2491 /*}}}*/
2492
2493 /*Fetch parameters: */
2494 iomodel->Constant(&yts,ConstantsYtsEnum);
2495 iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
2496 iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
2497 iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
2498
2499 /*Recover element type*/
2500 if ((analysis_type==PrognosticAnalysisEnum && progstabilization==3) || (analysis_type==BalancethicknessAnalysisEnum && balancestabilization==3)){
2501 /*P1 Discontinuous Galerkin*/
2502 tria_type=P1DGEnum;
2503 }
2504 else{
2505 /*P1 Continuous Galerkin*/
2506 tria_type=P1Enum;
2507 }
2508 this->SetElementType(tria_type,analysis_counter);
2509
2510 /*Recover vertices ids needed to initialize inputs*/
2511 for(i=0;i<3;i++){
2512 tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
2513 }
2514
2515 /*Recover nodes ids needed to initialize the node hook.*/
2516 if (tria_type==P1DGEnum){
2517 /*Discontinuous Galerkin*/
2518 tria_node_ids[0]=iomodel->nodecounter+3*index+1;
2519 tria_node_ids[1]=iomodel->nodecounter+3*index+2;
2520 tria_node_ids[2]=iomodel->nodecounter+3*index+3;
2521 }
2522 else{
2523 /*Continuous Galerkin*/
2524 for(i=0;i<3;i++){
2525 tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->Data(MeshElementsEnum)+3*index+i); //ids for vertices are in the elements array from Matlab
2526 }
2527 }
2528
2529 /*hooks: */
2530 this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
2531
2532 /*Fill with IoModel*/
2533 this->InputUpdateFromIoModel(index,iomodel);
2534
2535 /*Defaults if not provided in iomodel*/
2536 switch(analysis_type){
2537
2538 case DiagnosticHorizAnalysisEnum:
2539
2540 /*default vx,vy and vz: either observation or 0 */
2541 if(!iomodel->Data(VxEnum)){
2542 for(i=0;i<3;i++)nodeinputs[i]=0;
2543 this->inputs->AddInput(new TriaP1Input(VxEnum,nodeinputs));
2544 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVxEnum,nodeinputs));
2545 }
2546 if(!iomodel->Data(VyEnum)){
2547 for(i=0;i<3;i++)nodeinputs[i]=0;
2548 this->inputs->AddInput(new TriaP1Input(VyEnum,nodeinputs));
2549 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVyEnum,nodeinputs));
2550 }
2551 if(!iomodel->Data(VzEnum)){
2552 for(i=0;i<3;i++)nodeinputs[i]=0;
2553 this->inputs->AddInput(new TriaP1Input(VzEnum,nodeinputs));
2554 if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVzEnum,nodeinputs));
2555 }
2556 if(!iomodel->Data(PressureEnum)){
2557 for(i=0;i<3;i++)nodeinputs[i]=0;
2558 if(dakota_analysis){
2559 this->inputs->AddInput(new TriaP1Input(PressureEnum,nodeinputs));
2560 this->inputs->AddInput(new TriaP1Input(QmuPressureEnum,nodeinputs));
2561 }
2562 }
2563 break;
2564
2565 default:
2566 /*No update for other solution types*/
2567 break;
2568
2569 }
2570
2571 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
2572 this->parameters=NULL;
2573}
2574/*}}}*/
2575/*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{1*/
2576int Tria::UpdatePotentialSheetUngrounding(double* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,double* nodes_on_iceshelf){
2577
2578 int i;
2579 int nflipped=0;
2580
2581 /*Go through nodes, and whoever is on the potential_sheet_ungrounding, ends up in nodes_on_iceshelf: */
2582 for(i=0;i<3;i++){
2583 if (vertices_potentially_ungrounding[nodes[i]->Sid()]){
2584 vec_nodes_on_iceshelf->SetValue(nodes[i]->Sid(),1,INS_VAL);
2585
2586 /*If node was not on ice shelf, we flipped*/
2587 if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
2588 nflipped++;
2589 }
2590 }
2591 }
2592 return nflipped;
2593}
2594/*}}}*/
2595
2596#ifdef _HAVE_RESPONSES_
2597/*FUNCTION Tria::IceVolume {{{1*/
2598double Tria::IceVolume(void){
2599
2600 /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
2601 double base,surface,bed;
2602 double xyz_list[NUMVERTICES][3];
2603
2604 if(IsOnWater())return 0;
2605
2606 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2607
2608 /*First calculate the area of the base (cross section triangle)
2609 * http://en.wikipedia.org/wiki/Triangle
2610 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
2611 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]));
2612
2613 /*Now get the average height*/
2614 Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
2615 Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
2616 surface_input->GetInputAverage(&surface);
2617 bed_input->GetInputAverage(&bed);
2618
2619 /*Return: */
2620 return base*(surface-bed);
2621}
2622/*}}}*/
2623/*FUNCTION Tria::MassFlux {{{1*/
2624double Tria::MassFlux( double* segment,bool process_units){
2625
2626 const int numdofs=2;
2627
2628 int i,dim;
2629 double mass_flux=0;
2630 double xyz_list[NUMVERTICES][3];
2631 double normal[2];
2632 double length,rho_ice;
2633 double x1,y1,x2,y2,h1,h2;
2634 double vx1,vx2,vy1,vy2;
2635 GaussTria* gauss_1=NULL;
2636 GaussTria* gauss_2=NULL;
2637
2638 /*Get material parameters :*/
2639 rho_ice=matpar->GetRhoIce();
2640
2641 /*First off, check that this segment belongs to this element: */
2642 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);
2643
2644 /*Recover segment node locations: */
2645 x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
2646
2647 /*Get xyz list: */
2648 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2649
2650 /*get area coordinates of 0 and 1 locations: */
2651 gauss_1=new GaussTria();
2652 gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
2653 gauss_2=new GaussTria();
2654 gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
2655
2656 normal[0]=cos(atan2(x1-x2,y2-y1));
2657 normal[1]=sin(atan2(x1-x2,y2-y1));
2658
2659 length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));
2660
2661 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
2662 this->parameters->FindParam(&dim,MeshDimensionEnum);
2663 Input* vx_input=NULL;
2664 Input* vy_input=NULL;
2665 if(dim==2){
2666 vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2667 vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2668 }
2669 else{
2670 vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
2671 vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
2672 }
2673
2674 thickness_input->GetInputValue(&h1, gauss_1);
2675 thickness_input->GetInputValue(&h2, gauss_2);
2676 vx_input->GetInputValue(&vx1,gauss_1);
2677 vx_input->GetInputValue(&vx2,gauss_2);
2678 vy_input->GetInputValue(&vy1,gauss_1);
2679 vy_input->GetInputValue(&vy2,gauss_2);
2680
2681 mass_flux= rho_ice*length*(
2682 (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
2683 (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
2684 );
2685
2686 /*Process units: */
2687 mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum);
2688
2689 /*clean up and return:*/
2690 delete gauss_1;
2691 delete gauss_2;
2692 return mass_flux;
2693}
2694/*}}}*/
2695/*FUNCTION Tria::MaxAbsVx{{{1*/
2696void Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){
2697
2698 /*Get maximum:*/
2699 double maxabsvx=this->inputs->MaxAbs(VxEnum);
2700
2701 /*process units if requested: */
2702 if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum);
2703
2704 /*Assign output pointers:*/
2705 *pmaxabsvx=maxabsvx;
2706}
2707/*}}}*/
2708/*FUNCTION Tria::MaxAbsVy{{{1*/
2709void Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){
2710
2711 /*Get maximum:*/
2712 double maxabsvy=this->inputs->MaxAbs(VyEnum);
2713
2714 /*process units if requested: */
2715 if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum);
2716
2717 /*Assign output pointers:*/
2718 *pmaxabsvy=maxabsvy;
2719}
2720/*}}}*/
2721/*FUNCTION Tria::MaxAbsVz{{{1*/
2722void Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){
2723
2724 /*Get maximum:*/
2725 double maxabsvz=this->inputs->MaxAbs(VzEnum);
2726
2727 /*process units if requested: */
2728 if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum);
2729
2730 /*Assign output pointers:*/
2731 *pmaxabsvz=maxabsvz;
2732}
2733/*}}}*/
2734/*FUNCTION Tria::MaxVel{{{1*/
2735void Tria::MaxVel(double* pmaxvel, bool process_units){
2736
2737 /*Get maximum:*/
2738 double maxvel=this->inputs->Max(VelEnum);
2739
2740 /*process units if requested: */
2741 if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum);
2742
2743 /*Assign output pointers:*/
2744 *pmaxvel=maxvel;
2745}
2746/*}}}*/
2747/*FUNCTION Tria::MaxVx{{{1*/
2748void Tria::MaxVx(double* pmaxvx, bool process_units){
2749
2750 /*Get maximum:*/
2751 double maxvx=this->inputs->Max(VxEnum);
2752
2753 /*process units if requested: */
2754 if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum);
2755
2756 /*Assign output pointers:*/
2757 *pmaxvx=maxvx;
2758}
2759/*}}}*/
2760/*FUNCTION Tria::MaxVy{{{1*/
2761void Tria::MaxVy(double* pmaxvy, bool process_units){
2762
2763 /*Get maximum:*/
2764 double maxvy=this->inputs->Max(VyEnum);
2765
2766 /*process units if requested: */
2767 if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum);
2768
2769 /*Assign output pointers:*/
2770 *pmaxvy=maxvy;
2771
2772}
2773/*}}}*/
2774/*FUNCTION Tria::MaxVz{{{1*/
2775void Tria::MaxVz(double* pmaxvz, bool process_units){
2776
2777 /*Get maximum:*/
2778 double maxvz=this->inputs->Max(VzEnum);
2779
2780 /*process units if requested: */
2781 if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum);
2782
2783 /*Assign output pointers:*/
2784 *pmaxvz=maxvz;
2785}
2786/*}}}*/
2787/*FUNCTION Tria::MinVel{{{1*/
2788void Tria::MinVel(double* pminvel, bool process_units){
2789
2790 /*Get minimum:*/
2791 double minvel=this->inputs->Min(VelEnum);
2792
2793 /*process units if requested: */
2794 if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum);
2795
2796 /*Assign output pointers:*/
2797 *pminvel=minvel;
2798}
2799/*}}}*/
2800/*FUNCTION Tria::MinVx{{{1*/
2801void Tria::MinVx(double* pminvx, bool process_units){
2802
2803 /*Get minimum:*/
2804 double minvx=this->inputs->Min(VxEnum);
2805
2806 /*process units if requested: */
2807 if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum);
2808
2809 /*Assign output pointers:*/
2810 *pminvx=minvx;
2811}
2812/*}}}*/
2813/*FUNCTION Tria::MinVy{{{1*/
2814void Tria::MinVy(double* pminvy, bool process_units){
2815
2816 /*Get minimum:*/
2817 double minvy=this->inputs->Min(VyEnum);
2818
2819 /*process units if requested: */
2820 if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum);
2821
2822 /*Assign output pointers:*/
2823 *pminvy=minvy;
2824}
2825/*}}}*/
2826/*FUNCTION Tria::MinVz{{{1*/
2827void Tria::MinVz(double* pminvz, bool process_units){
2828
2829 /*Get minimum:*/
2830 double minvz=this->inputs->Min(VzEnum);
2831
2832 /*process units if requested: */
2833 if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum);
2834
2835 /*Assign output pointers:*/
2836 *pminvz=minvz;
2837}
2838/*}}}*/
2839/*FUNCTION Tria::ElementResponse{{{1*/
2840void Tria::ElementResponse(double* presponse,int response_enum,bool process_units){
2841
2842 switch(response_enum){
2843 case MaterialsRheologyBbarEnum:
2844 *presponse=this->matice->GetBbar();
2845 break;
2846 case MaterialsRheologyZbarEnum:
2847 *presponse=this->matice->GetZbar();
2848 break;
2849 case VelEnum:
2850
2851 /*Get input:*/
2852 double vel;
2853 Input* vel_input;
2854
2855 vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
2856 vel_input->GetInputAverage(&vel);
2857
2858 /*process units if requested: */
2859 if(process_units) vel=UnitConversion(vel,IuToExtEnum,VelEnum);
2860
2861 /*Assign output pointers:*/
2862 *presponse=vel;
2863 default:
2864 _error_("Response type %s not supported yet!",EnumToStringx(response_enum));
2865 }
2866
2867}
2868/*}}}*/
2869#endif
2870
2871#ifdef _HAVE_DIAGNOSTIC_
2872/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/
2873ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyeal(void){
2874
2875 /*compute all stiffness matrices for this element*/
2876 ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyealViscous();
2877 ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction();
2878 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
2879
2880 /*clean-up and return*/
2881 delete Ke1;
2882 delete Ke2;
2883 return Ke;
2884}
2885/*}}}*/
2886/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
2887ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){
2888
2889 /*Constants*/
2890 const int numdof=NDOF2*NUMVERTICES;
2891
2892 /*Intermediaries*/
2893 int i,j,ig;
2894 double xyz_list[NUMVERTICES][3];
2895 double viscosity,newviscosity,oldviscosity;
2896 double viscosity_overshoot,thickness,Jdet;
2897 double epsilon[3],oldepsilon[3]; /* epsilon=[exx,eyy,exy]; */
2898 double B[3][numdof];
2899 double Bprime[3][numdof];
2900 double D[3][3] = {0.0};
2901 double D_scalar;
2902 GaussTria *gauss = NULL;
2903
2904 /*Initialize Element matrix*/
2905 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
2906
2907 /*Retrieve all inputs and parameters*/
2908 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2909 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
2910 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2911 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2912 Input* vxold_input=inputs->GetInput(VxPicardEnum); _assert_(vxold_input);
2913 Input* vyold_input=inputs->GetInput(VyPicardEnum); _assert_(vyold_input);
2914 this->parameters->FindParam(&viscosity_overshoot,DiagnosticViscosityOvershootEnum);
2915
2916 /* Start looping on the number of gaussian points: */
2917 gauss=new GaussTria(2);
2918 for (ig=gauss->begin();ig<gauss->end();ig++){
2919
2920 gauss->GaussPoint(ig);
2921
2922 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
2923 GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
2924 GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);
2925
2926 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
2927 this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
2928 matice->GetViscosity2d(&viscosity, &epsilon[0]);
2929 matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
2930 thickness_input->GetInputValue(&thickness, gauss);
2931
2932 newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
2933 D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
2934 for (i=0;i<3;i++) D[i][i]=D_scalar;
2935
2936 TripleMultiply(&B[0][0],3,numdof,1,
2937 &D[0][0],3,3,0,
2938 &Bprime[0][0],3,numdof,0,
2939 &Ke->values[0],1);
2940 }
2941
2942 /*Transform Coordinate System*/
2943 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
2944
2945 /*Clean up and return*/
2946 delete gauss;
2947 return Ke;
2948}
2949/*}}}*/
2950/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
2951ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){
2952
2953 /*Constants*/
2954 const int numdof=NDOF2*NUMVERTICES;
2955
2956 /*Intermediaries*/
2957 int i,j,ig;
2958 int analysis_type;
2959 double MAXSLOPE = .06; // 6 %
2960 double MOUNTAINKEXPONENT = 10;
2961 double slope_magnitude,alpha2;
2962 double Jdet;
2963 double L[2][numdof];
2964 double DL[2][2] = {{ 0,0 },{0,0}};
2965 double DL_scalar;
2966 double slope[2] = {0.0,0.0};
2967 double xyz_list[NUMVERTICES][3];
2968 Friction *friction = NULL;
2969 GaussTria *gauss = NULL;
2970
2971 /*Initialize Element matrix and return if necessary*/
2972 if(IsFloating()) return NULL;
2973 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
2974
2975 /*Retrieve all inputs and parameters*/
2976 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
2977 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
2978 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2979 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2980 Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
2981 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
2982
2983 /*build friction object, used later on: */
2984 friction=new Friction("2d",inputs,matpar,analysis_type);
2985
2986 /* Start looping on the number of gaussian points: */
2987 gauss=new GaussTria(2);
2988 for (ig=gauss->begin();ig<gauss->end();ig++){
2989
2990 gauss->GaussPoint(ig);
2991
2992 // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case,
2993 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
2994 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
2995 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
2996 if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT);
2997 else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
2998
2999 GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
3000 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3001 DL_scalar=alpha2*gauss->weight*Jdet;
3002 for (i=0;i<2;i++) DL[i][i]=DL_scalar;
3003
3004 TripleMultiply( &L[0][0],2,numdof,1,
3005 &DL[0][0],2,2,0,
3006 &L[0][0],2,numdof,0,
3007 &Ke->values[0],1);
3008 }
3009
3010 /*Transform Coordinate System*/
3011 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
3012
3013 /*Clean up and return*/
3014 delete gauss;
3015 delete friction;
3016 return Ke;
3017}
3018/*}}}*/
3019/*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/
3020ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){
3021
3022 /*Intermediaries*/
3023 const int numdof=NUMVERTICES*NDOF2;
3024 int i,connectivity;
3025
3026 /*Initialize Element matrix*/
3027 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
3028
3029 /*Create Element matrix*/
3030 for(i=0;i<NUMVERTICES;i++){
3031 connectivity=nodes[i]->GetConnectivity();
3032 Ke->values[(2*i)*numdof +(2*i) ]=1/(double)connectivity;
3033 Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity;
3034 }
3035
3036 /*Clean up and return*/
3037 return Ke;
3038}
3039/*}}}*/
3040/*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/
3041ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){
3042
3043 /*Constants*/
3044 const int numdof=NDOF2*NUMVERTICES;
3045
3046 /*Intermediaries */
3047 int i,j,ig;
3048 double driving_stress_baseline,thickness;
3049 double Jdet;
3050 double xyz_list[NUMVERTICES][3];
3051 double slope[2];
3052 double basis[3];
3053 double pe_g_gaussian[numdof];
3054 GaussTria* gauss=NULL;
3055
3056 /*Initialize Element vector*/
3057 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
3058
3059 /*Retrieve all inputs and parameters*/
3060 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3061 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3062 Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
3063 Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
3064
3065 /* Start looping on the number of gaussian points: */
3066 gauss=new GaussTria(2);
3067 for(ig=gauss->begin();ig<gauss->end();ig++){
3068
3069 gauss->GaussPoint(ig);
3070
3071 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3072 GetNodalFunctions(basis, gauss);
3073
3074 thickness_input->GetInputValue(&thickness,gauss);
3075 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
3076 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
3077
3078 /*Build pe_g_gaussian vector: */
3079 for (i=0;i<NUMVERTICES;i++){
3080 for (j=0;j<NDOF2;j++){
3081 pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
3082 }
3083 }
3084 }
3085
3086 /*Transform coordinate system*/
3087 TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);
3088
3089 /*Clean up and return*/
3090 delete gauss;
3091 return pe;
3092}
3093/*}}}*/
3094/*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
3095ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
3096
3097 /*Intermediaries */
3098 int i,connectivity;
3099 double constant_part,ub,vb;
3100 double rho_ice,gravity,n,B;
3101 double slope2,thickness;
3102 double slope[2];
3103 GaussTria* gauss=NULL;
3104
3105 /*Initialize Element vector*/
3106 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
3107
3108 /*Retrieve all inputs and parameters*/
3109 rho_ice=matpar->GetRhoIce();
3110 gravity=matpar->GetG();
3111 n=matice->GetN();
3112 B=matice->GetBbar();
3113 Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
3114 Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
3115 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3116
3117 /*Spawn 3 sing elements: */
3118 gauss=new GaussTria();
3119 for(i=0;i<NUMVERTICES;i++){
3120
3121 gauss->GaussVertex(i);
3122
3123 connectivity=nodes[i]->GetConnectivity();
3124
3125 thickness_input->GetInputValue(&thickness,gauss);
3126 slopex_input->GetInputValue(&slope[0],gauss);
3127 slopey_input->GetInputValue(&slope[1],gauss);
3128 slope2=pow(slope[0],2)+pow(slope[1],2);
3129
3130 constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));
3131
3132 ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
3133 vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1];
3134
3135 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;
3136 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;
3137 }
3138
3139 /*Clean up and return*/
3140 delete gauss;
3141 return pe;
3142}
3143/*}}}*/
3144/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{1*/
3145ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
3146
3147 /*Constants*/
3148 const int numdof=NDOF2*NUMVERTICES;
3149
3150 /*Intermediaries */
3151 int i,j,ig;
3152 double xyz_list[NUMVERTICES][3];
3153 double Jdet,thickness;
3154 double eps1dotdphii,eps1dotdphij;
3155 double eps2dotdphii,eps2dotdphij;
3156 double mu_prime;
3157 double epsilon[3];/* epsilon=[exx,eyy,exy];*/
3158 double eps1[2],eps2[2];
3159 double phi[NUMVERTICES];
3160 double dphi[2][NUMVERTICES];
3161 GaussTria *gauss=NULL;
3162
3163 /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
3164 ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
3165
3166 /*Retrieve all inputs and parameters*/
3167 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3168 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3169 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3170 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3171
3172 /* Start looping on the number of gaussian points: */
3173 gauss=new GaussTria(2);
3174 for (ig=gauss->begin();ig<gauss->end();ig++){
3175
3176 gauss->GaussPoint(ig);
3177
3178 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3179 GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
3180
3181 thickness_input->GetInputValue(&thickness, gauss);
3182 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
3183 matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
3184 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
3185 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
3186
3187 for(i=0;i<3;i++){
3188 for(j=0;j<3;j++){
3189 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
3190 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
3191 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
3192 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
3193
3194 Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
3195 Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
3196 Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
3197 Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
3198 }
3199 }
3200 }
3201
3202 /*Transform Coordinate System*/
3203 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
3204
3205 /*Clean up and return*/
3206 delete gauss;
3207 return Ke;
3208}
3209/*}}}*/
3210/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
3211void Tria::GetSolutionFromInputsDiagnosticHoriz(Vector* solution){
3212
3213 const int numdof=NDOF2*NUMVERTICES;
3214
3215 int i;
3216 int* doflist=NULL;
3217 double vx,vy;
3218 double values[numdof];
3219 GaussTria* gauss=NULL;
3220
3221 /*Get dof list: */
3222 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3223
3224 /*Get inputs*/
3225 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3226 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3227
3228 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3229 /*P1 element only for now*/
3230 gauss=new GaussTria();
3231 for(i=0;i<NUMVERTICES;i++){
3232
3233 gauss->GaussVertex(i);
3234
3235 /*Recover vx and vy*/
3236 vx_input->GetInputValue(&vx,gauss);
3237 vy_input->GetInputValue(&vy,gauss);
3238 values[i*NDOF2+0]=vx;
3239 values[i*NDOF2+1]=vy;
3240 }
3241
3242 solution->SetValues(numdof,doflist,values,INS_VAL);
3243
3244 /*Free ressources:*/
3245 delete gauss;
3246 xfree((void**)&doflist);
3247}
3248/*}}}*/
3249/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
3250void Tria::GetSolutionFromInputsDiagnosticHutter(Vector* solution){
3251
3252 const int numdof=NDOF2*NUMVERTICES;
3253
3254 int i;
3255 double vx,vy;
3256 double values[numdof];
3257 int *doflist = NULL;
3258 GaussTria *gauss = NULL;
3259
3260 /*Get dof list: */
3261 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3262
3263 /*Get inputs*/
3264 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3265 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3266
3267 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3268 /*P1 element only for now*/
3269 gauss=new GaussTria();
3270 for(i=0;i<NUMVERTICES;i++){
3271
3272 gauss->GaussVertex(i);
3273
3274 /*Recover vx and vy*/
3275 vx_input->GetInputValue(&vx,gauss);
3276 vy_input->GetInputValue(&vy,gauss);
3277 values[i*NDOF2+0]=vx;
3278 values[i*NDOF2+1]=vy;
3279 }
3280
3281 solution->SetValues(numdof,doflist,values,INS_VAL);
3282
3283 /*Free ressources:*/
3284 delete gauss;
3285 xfree((void**)&doflist);
3286}
3287/*}}}*/
3288/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
3289void Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
3290
3291 const int numdof=NDOF2*NUMVERTICES;
3292
3293 int i;
3294 int* doflist=NULL;
3295 double rho_ice,g;
3296 double values[numdof];
3297 double vx[NUMVERTICES];
3298 double vy[NUMVERTICES];
3299 double vz[NUMVERTICES];
3300 double vel[NUMVERTICES];
3301 double pressure[NUMVERTICES];
3302 double thickness[NUMVERTICES];
3303
3304 /*Get dof list: */
3305 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3306
3307 /*Use the dof list to index into the solution vector: */
3308 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
3309
3310 /*Transform solution in Cartesian Space*/
3311 TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
3312
3313 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3314 for(i=0;i<NUMVERTICES;i++){
3315 vx[i]=values[i*NDOF2+0];
3316 vy[i]=values[i*NDOF2+1];
3317
3318 /*Check solution*/
3319 if(isnan(vx[i])) _error_("NaN found in solution vector");
3320 if(isnan(vy[i])) _error_("NaN found in solution vector");
3321 }
3322
3323 /*Get Vz and compute vel*/
3324 GetInputListOnVertices(&vz[0],VzEnum,0);
3325 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);
3326
3327 /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
3328 *so the pressure is just the pressure at the bedrock: */
3329 rho_ice=matpar->GetRhoIce();
3330 g=matpar->GetG();
3331 GetInputListOnVertices(&thickness[0],ThicknessEnum);
3332 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
3333
3334 /*Now, we have to move the previous Vx and Vy inputs to old
3335 * status, otherwise, we'll wipe them off: */
3336 this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
3337 this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
3338 this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
3339
3340 /*Add vx and vy as inputs to the tria element: */
3341 this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
3342 this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
3343 this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
3344 this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));
3345
3346 /*Free ressources:*/
3347 xfree((void**)&doflist);
3348
3349}
3350/*}}}*/
3351/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{1*/
3352void Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){
3353
3354 const int numdof=NDOF2*NUMVERTICES;
3355
3356 int i;
3357 int* doflist=NULL;
3358 double rho_ice,g;
3359 double values[numdof];
3360 double vx[NUMVERTICES];
3361 double vy[NUMVERTICES];
3362 double vz[NUMVERTICES];
3363 double vel[NUMVERTICES];
3364 double pressure[NUMVERTICES];
3365 double thickness[NUMVERTICES];
3366
3367 /*Get dof list: */
3368 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
3369
3370 /*Use the dof list to index into the solution vector: */
3371 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
3372
3373 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
3374 for(i=0;i<NUMVERTICES;i++){
3375 vx[i]=values[i*NDOF2+0];
3376 vy[i]=values[i*NDOF2+1];
3377
3378 /*Check solution*/
3379 if(isnan(vx[i])) _error_("NaN found in solution vector");
3380 if(isnan(vy[i])) _error_("NaN found in solution vector");
3381 }
3382
3383 /*Now Compute vel*/
3384 GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
3385 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);
3386
3387 /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
3388 *so the pressure is just the pressure at the bedrock: */
3389 rho_ice=matpar->GetRhoIce();
3390 g=matpar->GetG();
3391 GetInputListOnVertices(&thickness[0],ThicknessEnum);
3392 for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];
3393
3394 /*Now, we have to move the previous Vx and Vy inputs to old
3395 * status, otherwise, we'll wipe them off: */
3396 this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
3397 this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
3398 this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
3399
3400 /*Add vx and vy as inputs to the tria element: */
3401 this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
3402 this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
3403 this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
3404 this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));
3405
3406 /*Free ressources:*/
3407 xfree((void**)&doflist);
3408}
3409/*}}}*/
3410#endif
3411
3412#ifdef _HAVE_CONTROL_
3413/*FUNCTION Tria::InputControlUpdate{{{1*/
3414void Tria::InputControlUpdate(double scalar,bool save_parameter){
3415
3416 /*Intermediary*/
3417 int num_controls;
3418 int* control_type=NULL;
3419 Input* input=NULL;
3420
3421 /*retrieve some parameters: */
3422 this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
3423 this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
3424
3425 for(int i=0;i<num_controls;i++){
3426
3427 if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){
3428 input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
3429 }
3430 else{
3431 input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
3432 }
3433
3434 if (input->ObjectEnum()!=ControlInputEnum){
3435 _error_("input %s is not a ControlInput",EnumToStringx(control_type[i]));
3436 }
3437
3438 ((ControlInput*)input)->UpdateValue(scalar);
3439 ((ControlInput*)input)->Constrain();
3440 if (save_parameter) ((ControlInput*)input)->SaveValue();
3441
3442 }
3443
3444 /*Clean up and return*/
3445 xfree((void**)&control_type);
3446}
3447/*}}}*/
3448/*FUNCTION Tria::ControlInputGetGradient{{{1*/
3449void Tria::ControlInputGetGradient(Vector* gradient,int enum_type,int control_index){
3450
3451 int doflist1[NUMVERTICES];
3452 Input* input=NULL;
3453
3454 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
3455 input=(Input*)matice->inputs->GetInput(enum_type);
3456 }
3457 else{
3458 input=inputs->GetInput(enum_type);
3459 }
3460 if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
3461 if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
3462
3463 GradientIndexing(&doflist1[0],control_index);
3464 ((ControlInput*)input)->GetGradient(gradient,&doflist1[0]);
3465
3466}/*}}}*/
3467/*FUNCTION Tria::ControlInputScaleGradient{{{1*/
3468void Tria::ControlInputScaleGradient(int enum_type,double scale){
3469
3470 Input* input=NULL;
3471
3472 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
3473 input=(Input*)matice->inputs->GetInput(enum_type);
3474 }
3475 else{
3476 input=inputs->GetInput(enum_type);
3477 }
3478 if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
3479 if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
3480
3481 ((ControlInput*)input)->ScaleGradient(scale);
3482}/*}}}*/
3483/*FUNCTION Tria::ControlInputSetGradient{{{1*/
3484void Tria::ControlInputSetGradient(double* gradient,int enum_type,int control_index){
3485
3486 int doflist1[NUMVERTICES];
3487 double grad_list[NUMVERTICES];
3488 Input* grad_input=NULL;
3489 Input* input=NULL;
3490
3491 if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
3492 input=(Input*)matice->inputs->GetInput(enum_type);
3493 }
3494 else{
3495 input=inputs->GetInput(enum_type);
3496 }
3497 if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
3498 if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));
3499
3500 GradientIndexing(&doflist1[0],control_index);
3501 for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[doflist1[i]];
3502 grad_input=new TriaP1Input(GradientEnum,grad_list);
3503
3504 ((ControlInput*)input)->SetGradient(grad_input);
3505
3506}/*}}}*/
3507/*FUNCTION Tria::Gradj {{{1*/
3508void Tria::Gradj(Vector* gradient,int control_type,int control_index){
3509 /*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/
3510
3511 /*If on water, grad = 0: */
3512 if(IsOnWater()) return;
3513
3514 /*First deal with ∂/∂alpha(KU-F)*/
3515 switch(control_type){
3516 case FrictionCoefficientEnum:
3517 GradjDragMacAyeal(gradient,control_index);
3518 break;
3519 case MaterialsRheologyBbarEnum:
3520 GradjBMacAyeal(gradient,control_index);
3521 break;
3522 case MaterialsRheologyZbarEnum:
3523 GradjZMacAyeal(gradient,control_index);
3524 break;
3525 case BalancethicknessThickeningRateEnum:
3526 GradjDhDtBalancedthickness(gradient,control_index);
3527 break;
3528 case VxEnum:
3529 GradjVxBalancedthickness(gradient,control_index);
3530 break;
3531 case VyEnum:
3532 GradjVyBalancedthickness(gradient,control_index);
3533 break;
3534 default:
3535 _error_("%s%i","control type not supported yet: ",control_type);
3536 }
3537
3538 /*Now deal with ∂J/∂alpha*/
3539 int *responses = NULL;
3540 int num_responses,resp;
3541 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
3542 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
3543
3544 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
3545 //FIXME: the control type should be checked somewhere (with respect to what variable are we taking the gradient!)
3546
3547 case ThicknessAbsMisfitEnum:
3548 case ThicknessAbsGradientEnum:
3549 case SurfaceAbsVelMisfitEnum:
3550 case SurfaceRelVelMisfitEnum:
3551 case SurfaceLogVelMisfitEnum:
3552 case SurfaceLogVxVyMisfitEnum:
3553 case SurfaceAverageVelMisfitEnum:
3554 /*Nothing, J does not depends on the parameter being inverted for*/
3555 break;
3556 case DragCoefficientAbsGradientEnum:
3557 GradjDragGradient(gradient,resp,control_index);
3558 break;
3559 case RheologyBbarAbsGradientEnum:
3560 GradjBGradient(gradient,resp,control_index);
3561 break;
3562 default:
3563 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
3564 }
3565
3566 xfree((void**)&responses);
3567}
3568/*}}}*/
3569/*FUNCTION Tria::GradjBGradient{{{1*/
3570void Tria::GradjBGradient(Vector* gradient,int weight_index,int control_index){
3571
3572 int i,ig;
3573 int doflist1[NUMVERTICES];
3574 double Jdet,weight;
3575 double xyz_list[NUMVERTICES][3];
3576 double dbasis[NDOF2][NUMVERTICES];
3577 double dk[NDOF2];
3578 double grade_g[NUMVERTICES]={0.0};
3579 GaussTria *gauss=NULL;
3580
3581 /*Retrieve all inputs we will be needing: */
3582 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3583 GradientIndexing(&doflist1[0],control_index);
3584 Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
3585 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
3586
3587 /* Start looping on the number of gaussian points: */
3588 gauss=new GaussTria(2);
3589 for (ig=gauss->begin();ig<gauss->end();ig++){
3590
3591 gauss->GaussPoint(ig);
3592
3593 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3594 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
3595 weights_input->GetInputValue(&weight,gauss,weight_index);
3596
3597 /*Build alpha_complement_list: */
3598 rheologyb_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3599
3600 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3601 for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
3602 }
3603 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3604
3605 /*Clean up and return*/
3606 delete gauss;
3607}
3608/*}}}*/
3609/*FUNCTION Tria::GradjZGradient{{{1*/
3610void Tria::GradjZGradient(Vector* gradient,int weight_index,int control_index){
3611
3612 int i,ig;
3613 int doflist1[NUMVERTICES];
3614 double Jdet,weight;
3615 double xyz_list[NUMVERTICES][3];
3616 double dbasis[NDOF2][NUMVERTICES];
3617 double dk[NDOF2];
3618 double grade_g[NUMVERTICES]={0.0};
3619 GaussTria *gauss=NULL;
3620
3621 /*Retrieve all inputs we will be needing: */
3622 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3623 GradientIndexing(&doflist1[0],control_index);
3624 Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
3625 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
3626
3627 /* Start looping on the number of gaussian points: */
3628 gauss=new GaussTria(2);
3629 for (ig=gauss->begin();ig<gauss->end();ig++){
3630
3631 gauss->GaussPoint(ig);
3632
3633 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3634 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
3635 weights_input->GetInputValue(&weight,gauss,weight_index);
3636
3637 /*Build alpha_complement_list: */
3638 rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3639
3640 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3641 for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
3642 }
3643 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3644
3645 /*Clean up and return*/
3646 delete gauss;
3647}
3648/*}}}*/
3649/*FUNCTION Tria::GradjBMacAyeal{{{1*/
3650void Tria::GradjBMacAyeal(Vector* gradient,int control_index){
3651
3652 /*Intermediaries*/
3653 int i,ig;
3654 int doflist[NUMVERTICES];
3655 double vx,vy,lambda,mu,thickness,Jdet;
3656 double viscosity_complement;
3657 double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2];
3658 double xyz_list[NUMVERTICES][3];
3659 double basis[3],epsilon[3];
3660 double grad[NUMVERTICES]={0.0};
3661 GaussTria *gauss = NULL;
3662
3663 /* Get node coordinates and dof list: */
3664 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3665 GradientIndexing(&doflist[0],control_index);
3666
3667 /*Retrieve all inputs*/
3668 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3669 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3670 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3671 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input);
3672 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input);
3673 Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
3674
3675 /* Start looping on the number of gaussian points: */
3676 gauss=new GaussTria(4);
3677 for (ig=gauss->begin();ig<gauss->end();ig++){
3678
3679 gauss->GaussPoint(ig);
3680
3681 thickness_input->GetInputValue(&thickness,gauss);
3682 rheologyb_input->GetInputDerivativeValue(&dB[0],&xyz_list[0][0],gauss);
3683 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
3684 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
3685 adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
3686 adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
3687
3688 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
3689 matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);
3690
3691 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3692 GetNodalFunctions(basis,gauss);
3693
3694 /*standard gradient dJ/dki*/
3695 for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
3696 (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
3697 )*Jdet*gauss->weight*basis[i];
3698 }
3699
3700 gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
3701
3702 /*clean-up*/
3703 delete gauss;
3704}
3705/*}}}*/
3706/*FUNCTION Tria::GradjZMacAyeal{{{1*/
3707void Tria::GradjZMacAyeal(Vector* gradient,int control_index){
3708
3709 /*Intermediaries*/
3710 int i,ig;
3711 int doflist[NUMVERTICES];
3712 double vx,vy,lambda,mu,thickness,Jdet;
3713 double viscosity_complement;
3714 double dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2];
3715 double xyz_list[NUMVERTICES][3];
3716 double basis[3],epsilon[3];
3717 double grad[NUMVERTICES]={0.0};
3718 GaussTria *gauss = NULL;
3719
3720 /* Get node coordinates and dof list: */
3721 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3722 GradientIndexing(&doflist[0],control_index);
3723
3724 /*Retrieve all inputs*/
3725 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3726 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3727 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3728 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input);
3729 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input);
3730 Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
3731
3732 /* Start looping on the number of gaussian points: */
3733 gauss=new GaussTria(4);
3734 for (ig=gauss->begin();ig<gauss->end();ig++){
3735
3736 gauss->GaussPoint(ig);
3737
3738 thickness_input->GetInputValue(&thickness,gauss);
3739 rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss);
3740 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
3741 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
3742 adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
3743 adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
3744
3745 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
3746 matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]);
3747
3748 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3749 GetNodalFunctions(basis,gauss);
3750
3751 /*standard gradient dJ/dki*/
3752 for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
3753 (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
3754 )*Jdet*gauss->weight*basis[i];
3755 }
3756
3757 gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
3758
3759 /*clean-up*/
3760 delete gauss;
3761}
3762/*}}}*/
3763/*FUNCTION Tria::GradjDragMacAyeal {{{1*/
3764void Tria::GradjDragMacAyeal(Vector* gradient,int control_index){
3765
3766 int i,ig;
3767 int analysis_type;
3768 int doflist1[NUMVERTICES];
3769 int connectivity[NUMVERTICES];
3770 double vx,vy,lambda,mu,alpha_complement,Jdet;
3771 double bed,thickness,Neff,drag;
3772 double xyz_list[NUMVERTICES][3];
3773 double dk[NDOF2];
3774 double grade_g[NUMVERTICES]={0.0};
3775 double grade_g_gaussian[NUMVERTICES];
3776 double basis[3];
3777 double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
3778 Friction* friction=NULL;
3779 GaussTria *gauss=NULL;
3780
3781 if(IsFloating())return;
3782
3783 /*retrive parameters: */
3784 parameters->FindParam(&analysis_type,AnalysisTypeEnum);
3785 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3786 GradientIndexing(&doflist1[0],control_index);
3787 this->GetConnectivityList(&connectivity[0]);
3788
3789 /*Build friction element, needed later: */
3790 friction=new Friction("2d",inputs,matpar,analysis_type);
3791
3792 /*Retrieve all inputs we will be needing: */
3793 Input* adjointx_input=inputs->GetInput(AdjointxEnum); _assert_(adjointx_input);
3794 Input* adjointy_input=inputs->GetInput(AdjointyEnum); _assert_(adjointy_input);
3795 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
3796 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
3797 Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
3798
3799 /* Start looping on the number of gaussian points: */
3800 gauss=new GaussTria(4);
3801 for (ig=gauss->begin();ig<gauss->end();ig++){
3802
3803 gauss->GaussPoint(ig);
3804
3805 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3806 GetNodalFunctions(basis, gauss);
3807
3808 /*Build alpha_complement_list: */
3809 friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
3810
3811 dragcoefficient_input->GetInputValue(&drag, gauss);
3812 adjointx_input->GetInputValue(&lambda, gauss);
3813 adjointy_input->GetInputValue(&mu, gauss);
3814 vx_input->GetInputValue(&vx,gauss);
3815 vy_input->GetInputValue(&vy,gauss);
3816 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3817
3818 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3819 for (i=0;i<NUMVERTICES;i++){
3820 grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
3821 }
3822
3823 /*Add gradje_g_gaussian vector to gradje_g: */
3824 for(i=0;i<NUMVERTICES;i++){
3825 _assert_(!isnan(grade_g[i]));
3826 grade_g[i]+=grade_g_gaussian[i];
3827 }
3828 }
3829 /*Analytical gradient*/
3830 //delete gauss;
3831 //gauss=new GaussTria();
3832 //for (int iv=0;iv<NUMVERTICES;iv++){
3833 // gauss->GaussVertex(iv);
3834 // friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
3835 // dragcoefficient_input->GetInputValue(&drag, gauss);
3836 // adjointx_input->GetInputValue(&lambda, gauss);
3837 // adjointy_input->GetInputValue(&mu, gauss);
3838 // vx_input->GetInputValue(&vx,gauss);
3839 // vy_input->GetInputValue(&vy,gauss);
3840 // grade_g[iv] = -2*1.e+7*drag*alpha_complement*(lambda*vx+mu*vy)/((double)connectivity[iv]);
3841 //}
3842 /*End Analytical gradient*/
3843
3844 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3845
3846 /*Clean up and return*/
3847 delete gauss;
3848 delete friction;
3849}
3850/*}}}*/
3851/*FUNCTION Tria::GradjDragGradient{{{1*/
3852void Tria::GradjDragGradient(Vector* gradient, int weight_index,int control_index){
3853
3854 int i,ig;
3855 int doflist1[NUMVERTICES];
3856 double Jdet,weight;
3857 double xyz_list[NUMVERTICES][3];
3858 double dbasis[NDOF2][NUMVERTICES];
3859 double dk[NDOF2];
3860 double grade_g[NUMVERTICES]={0.0};
3861 GaussTria *gauss=NULL;
3862
3863 /*Retrieve all inputs we will be needing: */
3864 if(IsFloating())return;
3865 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3866 GradientIndexing(&doflist1[0],control_index);
3867 Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
3868 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
3869
3870 /* Start looping on the number of gaussian points: */
3871 gauss=new GaussTria(2);
3872 for (ig=gauss->begin();ig<gauss->end();ig++){
3873
3874 gauss->GaussPoint(ig);
3875
3876 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3877 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
3878 weights_input->GetInputValue(&weight,gauss,weight_index);
3879
3880 /*Build alpha_complement_list: */
3881 dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
3882
3883 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
3884 for (i=0;i<NUMVERTICES;i++){
3885 grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
3886 _assert_(!isnan(grade_g[i]));
3887 }
3888 }
3889 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3890
3891 /*Clean up and return*/
3892 delete gauss;
3893}
3894/*}}}*/
3895/*FUNCTION Tria::GradjDhDtBalancedthickness{{{1*/
3896void Tria::GradjDhDtBalancedthickness(Vector* gradient,int control_index){
3897
3898 /*Intermediaries*/
3899 int doflist1[NUMVERTICES];
3900 double lambda[NUMVERTICES];
3901 double gradient_g[NUMVERTICES];
3902
3903 /*Compute Gradient*/
3904 GradientIndexing(&doflist1[0],control_index);
3905 GetInputListOnVertices(&lambda[0],AdjointEnum);
3906 for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i];
3907
3908 gradient->SetValues(NUMVERTICES,doflist1,gradient_g,INS_VAL);
3909}
3910/*}}}*/
3911/*FUNCTION Tria::GradjVxBalancedthickness{{{1*/
3912void Tria::GradjVxBalancedthickness(Vector* gradient,int control_index){
3913
3914 /*Intermediaries*/
3915 int i,ig;
3916 int doflist1[NUMVERTICES];
3917 double thickness,Jdet;
3918 double basis[3];
3919 double Dlambda[2],dp[2];
3920 double xyz_list[NUMVERTICES][3];
3921 double grade_g[NUMVERTICES] = {0.0};
3922 GaussTria *gauss = NULL;
3923
3924 /* Get node coordinates and dof list: */
3925 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3926 GradientIndexing(&doflist1[0],control_index);
3927
3928 /*Retrieve all inputs we will be needing: */
3929 Input* adjoint_input=inputs->GetInput(AdjointEnum); _assert_(adjoint_input);
3930 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3931
3932 /* Start looping on the number of gaussian points: */
3933 gauss=new GaussTria(2);
3934 for (ig=gauss->begin();ig<gauss->end();ig++){
3935
3936 gauss->GaussPoint(ig);
3937
3938 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3939 GetNodalFunctions(basis, gauss);
3940
3941 adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
3942 thickness_input->GetInputValue(&thickness, gauss);
3943 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
3944
3945 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
3946 }
3947
3948 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3949
3950 /*Clean up and return*/
3951 delete gauss;
3952}
3953/*}}}*/
3954/*FUNCTION Tria::GradjVyBalancedthickness{{{1*/
3955void Tria::GradjVyBalancedthickness(Vector* gradient,int control_index){
3956
3957 /*Intermediaries*/
3958 int i,ig;
3959 int doflist1[NUMVERTICES];
3960 double thickness,Jdet;
3961 double basis[3];
3962 double Dlambda[2],dp[2];
3963 double xyz_list[NUMVERTICES][3];
3964 double grade_g[NUMVERTICES] = {0.0};
3965 GaussTria *gauss = NULL;
3966
3967 /* Get node coordinates and dof list: */
3968 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
3969 GradientIndexing(&doflist1[0],control_index);
3970
3971 /*Retrieve all inputs we will be needing: */
3972 Input* adjoint_input=inputs->GetInput(AdjointEnum); _assert_(adjoint_input);
3973 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
3974
3975 /* Start looping on the number of gaussian points: */
3976 gauss=new GaussTria(2);
3977 for (ig=gauss->begin();ig<gauss->end();ig++){
3978
3979 gauss->GaussPoint(ig);
3980
3981 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
3982 GetNodalFunctions(basis, gauss);
3983
3984 adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
3985 thickness_input->GetInputValue(&thickness, gauss);
3986 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
3987
3988 for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
3989 }
3990 gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
3991
3992 /*Clean up and return*/
3993 delete gauss;
3994}
3995/*}}}*/
3996/*FUNCTION Tria::GradientIndexing{{{1*/
3997void Tria::GradientIndexing(int* indexing,int control_index){
3998
3999 /*Get some parameters*/
4000 int num_controls;
4001 parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
4002
4003 /*get gradient indices*/
4004 for(int i=0;i<NUMVERTICES;i++){
4005 indexing[i]=num_controls*this->nodes[i]->GetVertexDof() + control_index;
4006 }
4007
4008}
4009/*}}}*/
4010/*FUNCTION Tria::RheologyBbarAbsGradient{{{1*/
4011double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){
4012
4013 /* Intermediaries */
4014 int ig;
4015 double Jelem = 0;
4016 double weight;
4017 double Jdet;
4018 double xyz_list[NUMVERTICES][3];
4019 double dp[NDOF2];
4020 GaussTria *gauss = NULL;
4021
4022 /*retrieve parameters and inputs*/
4023
4024 /*If on water, return 0: */
4025 if(IsOnWater()) return 0;
4026
4027 /*Retrieve all inputs we will be needing: */
4028 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4029 Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4030 Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
4031
4032 /* Start looping on the number of gaussian points: */
4033 gauss=new GaussTria(2);
4034 for (ig=gauss->begin();ig<gauss->end();ig++){
4035
4036 gauss->GaussPoint(ig);
4037
4038 /* Get Jacobian determinant: */
4039 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4040
4041 /*Get all parameters at gaussian point*/
4042 weights_input->GetInputValue(&weight,gauss,weight_index);
4043 rheologyb_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
4044
4045 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
4046 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
4047 }
4048
4049 /*Clean up and return*/
4050 delete gauss;
4051 return Jelem;
4052}
4053/*}}}*/
4054/*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/
4055double Tria::SurfaceAverageVelMisfit(bool process_units,int weight_index){
4056
4057 const int numdof=2*NUMVERTICES;
4058
4059 int i,ig;
4060 double Jelem=0,S,Jdet;
4061 double misfit;
4062 double vx,vy,vxobs,vyobs,weight;
4063 double xyz_list[NUMVERTICES][3];
4064 GaussTria *gauss=NULL;
4065
4066 /*If on water, return 0: */
4067 if(IsOnWater())return 0;
4068
4069 /* Get node coordinates and dof list: */
4070 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4071
4072 /*Retrieve all inputs we will be needing: */
4073 inputs->GetInputValue(&S,SurfaceAreaEnum);
4074 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4075 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4076 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4077 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4078 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4079
4080 /* Start looping on the number of gaussian points: */
4081 gauss=new GaussTria(3);
4082 for (ig=gauss->begin();ig<gauss->end();ig++){
4083
4084 gauss->GaussPoint(ig);
4085
4086 /* Get Jacobian determinant: */
4087 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4088
4089 /*Get all parameters at gaussian point*/
4090 weights_input->GetInputValue(&weight,gauss,weight_index);
4091 vx_input->GetInputValue(&vx,gauss);
4092 vy_input->GetInputValue(&vy,gauss);
4093 vxobs_input->GetInputValue(&vxobs,gauss);
4094 vyobs_input->GetInputValue(&vyobs,gauss);
4095
4096 /*Compute SurfaceAverageVelMisfitEnum:
4097 *
4098 * 1 2 2
4099 * J = --- sqrt( (u - u ) + (v - v ) )
4100 * S obs obs
4101 */
4102 misfit=1/S*pow( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) ,0.5);
4103
4104 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
4105
4106 /*Add to cost function*/
4107 Jelem+=misfit*weight*Jdet*gauss->weight;
4108 }
4109
4110 /*clean-up and Return: */
4111 delete gauss;
4112 return Jelem;
4113}
4114/*}}}*/
4115/*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/
4116double Tria::SurfaceLogVelMisfit(bool process_units,int weight_index){
4117
4118 const int numdof=NDOF2*NUMVERTICES;
4119
4120 int i,ig;
4121 double Jelem=0;
4122 double misfit,Jdet;
4123 double epsvel=2.220446049250313e-16;
4124 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4125 double velocity_mag,obs_velocity_mag;
4126 double xyz_list[NUMVERTICES][3];
4127 double vx,vy,vxobs,vyobs,weight;
4128 GaussTria *gauss=NULL;
4129
4130 /*If on water, return 0: */
4131 if(IsOnWater())return 0;
4132
4133 /* Get node coordinates and dof list: */
4134 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4135
4136 /*Retrieve all inputs we will be needing: */
4137 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4138 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4139 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4140 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4141 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4142
4143 /* Start looping on the number of gaussian points: */
4144 gauss=new GaussTria(4);
4145 for (ig=gauss->begin();ig<gauss->end();ig++){
4146
4147 gauss->GaussPoint(ig);
4148
4149 /* Get Jacobian determinant: */
4150 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4151
4152 /*Get all parameters at gaussian point*/
4153 weights_input->GetInputValue(&weight,gauss,weight_index);
4154 vx_input->GetInputValue(&vx,gauss);
4155 vy_input->GetInputValue(&vy,gauss);
4156 vxobs_input->GetInputValue(&vxobs,gauss);
4157 vyobs_input->GetInputValue(&vyobs,gauss);
4158
4159 /*Compute SurfaceLogVelMisfit:
4160 * [ vel + eps ] 2
4161 * J = 4 \bar{v}^2 | log ( ----------- ) |
4162 * [ vel + eps ]
4163 * obs
4164 */
4165 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel;
4166 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
4167 misfit=4*pow(meanvel,2.)*pow(log(velocity_mag/obs_velocity_mag),2.);
4168
4169 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVelMisfitEnum);
4170
4171 /*Add to cost function*/
4172 Jelem+=misfit*weight*Jdet*gauss->weight;
4173 }
4174
4175 /*clean-up and Return: */
4176 delete gauss;
4177 return Jelem;
4178}
4179/*}}}*/
4180/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
4181double Tria::SurfaceLogVxVyMisfit(bool process_units,int weight_index){
4182
4183 const int numdof=NDOF2*NUMVERTICES;
4184
4185 int i,ig;
4186 int fit=-1;
4187 double Jelem=0, S=0;
4188 double epsvel=2.220446049250313e-16;
4189 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4190 double misfit, Jdet;
4191 double vx,vy,vxobs,vyobs,weight;
4192 double xyz_list[NUMVERTICES][3];
4193 GaussTria *gauss=NULL;
4194
4195 /*If on water, return 0: */
4196 if(IsOnWater())return 0;
4197
4198 /* Get node coordinates and dof list: */
4199 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4200
4201 /*Retrieve all inputs we will be needing: */
4202 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4203 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4204 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4205 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4206 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4207
4208 /* Start looping on the number of gaussian points: */
4209 gauss=new GaussTria(4);
4210 for (ig=gauss->begin();ig<gauss->end();ig++){
4211
4212 gauss->GaussPoint(ig);
4213
4214 /* Get Jacobian determinant: */
4215 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4216
4217 /*Get all parameters at gaussian point*/
4218 weights_input->GetInputValue(&weight,gauss,weight_index);
4219 vx_input->GetInputValue(&vx,gauss);
4220 vy_input->GetInputValue(&vy,gauss);
4221 vxobs_input->GetInputValue(&vxobs,gauss);
4222 vyobs_input->GetInputValue(&vyobs,gauss);
4223
4224 /*Compute SurfaceRelVelMisfit:
4225 *
4226 * 1 [ |u| + eps 2 |v| + eps 2 ]
4227 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
4228 * 2 [ |u |+ eps |v |+ eps ]
4229 * obs obs
4230 */
4231 misfit=0.5*pow(meanvel,2.)*(
4232 pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2.) +
4233 pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2.) );
4234
4235 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVxVyMisfitEnum);
4236
4237 /*Add to cost function*/
4238 Jelem+=misfit*weight*Jdet*gauss->weight;
4239 }
4240
4241 /*clean-up and Return: */
4242 delete gauss;
4243 return Jelem;
4244}
4245/*}}}*/
4246/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/
4247double Tria::SurfaceAbsVelMisfit(bool process_units,int weight_index){
4248
4249 const int numdof=NDOF2*NUMVERTICES;
4250
4251 int i,ig;
4252 double Jelem=0;
4253 double misfit,Jdet;
4254 double vx,vy,vxobs,vyobs,weight;
4255 double xyz_list[NUMVERTICES][3];
4256 GaussTria *gauss=NULL;
4257
4258 /*If on water, return 0: */
4259 if(IsOnWater())return 0;
4260
4261 /* Get node coordinates and dof list: */
4262 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4263
4264 /*Retrieve all inputs we will be needing: */
4265 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4266 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4267 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4268 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4269 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4270
4271 /* Start looping on the number of gaussian points: */
4272 gauss=new GaussTria(2);
4273 for (ig=gauss->begin();ig<gauss->end();ig++){
4274
4275 gauss->GaussPoint(ig);
4276
4277 /* Get Jacobian determinant: */
4278 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4279
4280 /*Get all parameters at gaussian point*/
4281 weights_input->GetInputValue(&weight,gauss,weight_index);
4282 vx_input->GetInputValue(&vx,gauss);
4283 vy_input->GetInputValue(&vy,gauss);
4284 vxobs_input->GetInputValue(&vxobs,gauss);
4285 vyobs_input->GetInputValue(&vyobs,gauss);
4286
4287 /*Compute SurfaceAbsVelMisfitEnum:
4288 *
4289 * 1 [ 2 2 ]
4290 * J = --- | (u - u ) + (v - v ) |
4291 * 2 [ obs obs ]
4292 *
4293 */
4294 misfit=0.5*( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) );
4295
4296 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);
4297
4298 /*Add to cost function*/
4299 Jelem+=misfit*weight*Jdet*gauss->weight;
4300 }
4301
4302 /*clean up and Return: */
4303 delete gauss;
4304 return Jelem;
4305}
4306/*}}}*/
4307/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
4308double Tria::SurfaceRelVelMisfit(bool process_units,int weight_index){
4309 const int numdof=2*NUMVERTICES;
4310
4311 int i,ig;
4312 double Jelem=0;
4313 double scalex=1,scaley=1;
4314 double misfit,Jdet;
4315 double epsvel=2.220446049250313e-16;
4316 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4317 double vx,vy,vxobs,vyobs,weight;
4318 double xyz_list[NUMVERTICES][3];
4319 GaussTria *gauss=NULL;
4320
4321 /*If on water, return 0: */
4322 if(IsOnWater())return 0;
4323
4324 /* Get node coordinates and dof list: */
4325 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4326
4327 /*Retrieve all inputs we will be needing: */
4328 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4329 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4330 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4331 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4332 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4333
4334 /* Start looping on the number of gaussian points: */
4335 gauss=new GaussTria(4);
4336 for (ig=gauss->begin();ig<gauss->end();ig++){
4337
4338 gauss->GaussPoint(ig);
4339
4340 /* Get Jacobian determinant: */
4341 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4342
4343 /*Get all parameters at gaussian point*/
4344 weights_input->GetInputValue(&weight,gauss,weight_index);
4345 vx_input->GetInputValue(&vx,gauss);
4346 vy_input->GetInputValue(&vy,gauss);
4347 vxobs_input->GetInputValue(&vxobs,gauss);
4348 vyobs_input->GetInputValue(&vyobs,gauss);
4349
4350 /*Compute SurfaceRelVelMisfit:
4351 *
4352 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
4353 * J = --- | ------------- (u - u ) + ------------- (v - v ) |
4354 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
4355 * obs obs
4356 */
4357 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
4358 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
4359 misfit=0.5*(scalex*pow((vx-vxobs),2.)+scaley*pow((vy-vyobs),2.));
4360 if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceRelVelMisfitEnum);
4361
4362 /*Add to cost function*/
4363 Jelem+=misfit*weight*Jdet*gauss->weight;
4364 }
4365
4366 /*clean up and Return: */
4367 delete gauss;
4368 return Jelem;
4369}
4370/*}}}*/
4371/*FUNCTION Tria::ThicknessAbsGradient{{{1*/
4372double Tria::ThicknessAbsGradient(bool process_units,int weight_index){
4373
4374 /* Intermediaries */
4375 int ig;
4376 double Jelem = 0;
4377 double weight;
4378 double Jdet;
4379 double xyz_list[NUMVERTICES][3];
4380 double dp[NDOF2];
4381 GaussTria *gauss = NULL;
4382
4383 /*retrieve parameters and inputs*/
4384
4385 /*If on water, return 0: */
4386 if(IsOnWater()) return 0;
4387
4388 /*Retrieve all inputs we will be needing: */
4389 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4390 Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4391 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4392
4393 /* Start looping on the number of gaussian points: */
4394 gauss=new GaussTria(2);
4395 for (ig=gauss->begin();ig<gauss->end();ig++){
4396
4397 gauss->GaussPoint(ig);
4398
4399 /* Get Jacobian determinant: */
4400 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4401
4402 /*Get all parameters at gaussian point*/
4403 weights_input->GetInputValue(&weight,gauss,weight_index);
4404 thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
4405
4406 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
4407 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
4408 }
4409
4410 /*Clean up and return*/
4411 delete gauss;
4412 return Jelem;
4413}
4414/*}}}*/
4415/*FUNCTION Tria::ThicknessAbsMisfit {{{1*/
4416double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){
4417
4418 /*Intermediaries*/
4419 int i,ig;
4420 double thickness,thicknessobs,weight;
4421 double Jdet;
4422 double Jelem = 0;
4423 double xyz_list[NUMVERTICES][3];
4424 GaussTria *gauss = NULL;
4425 double dH[2];
4426
4427 /*If on water, return 0: */
4428 if(IsOnWater())return 0;
4429
4430 /*Retrieve all inputs we will be needing: */
4431 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4432 Input* thickness_input =inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4433 Input* thicknessobs_input=inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
4434 Input* weights_input =inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4435
4436 /* Start looping on the number of gaussian points: */
4437 gauss=new GaussTria(2);
4438 for (ig=gauss->begin();ig<gauss->end();ig++){
4439
4440 gauss->GaussPoint(ig);
4441
4442 /* Get Jacobian determinant: */
4443 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4444
4445 /*Get parameters at gauss point*/
4446 thickness_input->GetInputValue(&thickness,gauss);
4447 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
4448 thicknessobs_input->GetInputValue(&thicknessobs,gauss);
4449 weights_input->GetInputValue(&weight,gauss,weight_index);
4450
4451 /*compute ThicknessAbsMisfit*/
4452 Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight;
4453 }
4454
4455 /* clean up and Return: */
4456 delete gauss;
4457 return Jelem;
4458}
4459/*}}}*/
4460/*FUNCTION Tria::CreatePVectorAdjointBalancethickness{{{1*/
4461ElementVector* Tria::CreatePVectorAdjointBalancethickness(void){
4462
4463 /*Constants*/
4464 const int numdof=1*NUMVERTICES;
4465
4466 /*Intermediaries */
4467 int i,ig,resp;
4468 double Jdet;
4469 double thickness,thicknessobs,weight;
4470 int *responses = NULL;
4471 int num_responses;
4472 double xyz_list[NUMVERTICES][3];
4473 double basis[3];
4474 double dbasis[NDOF2][NUMVERTICES];
4475 double dH[2];
4476 GaussTria* gauss=NULL;
4477
4478 /*Initialize Element vector*/
4479 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
4480
4481 /*Retrieve all inputs and parameters*/
4482 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4483 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
4484 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
4485 Input* thickness_input = inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4486 Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
4487 Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4488
4489 /* Start looping on the number of gaussian points: */
4490 gauss=new GaussTria(2);
4491 for(ig=gauss->begin();ig<gauss->end();ig++){
4492
4493 gauss->GaussPoint(ig);
4494
4495 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4496 GetNodalFunctions(basis, gauss);
4497 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
4498
4499 thickness_input->GetInputValue(&thickness, gauss);
4500 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
4501 thicknessobs_input->GetInputValue(&thicknessobs, gauss);
4502
4503 /*Loop over all requested responses*/
4504 for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
4505
4506 case ThicknessAbsMisfitEnum:
4507 weights_input->GetInputValue(&weight, gauss,resp);
4508 for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
4509 break;
4510 case ThicknessAbsGradientEnum:
4511 weights_input->GetInputValue(&weight, gauss,resp);
4512 for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
4513 for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
4514 break;
4515 default:
4516 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
4517 }
4518 }
4519
4520 /*Clean up and return*/
4521 delete gauss;
4522 xfree((void**)&responses);
4523 return pe;
4524}
4525/*}}}*/
4526/*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/
4527ElementVector* Tria::CreatePVectorAdjointHoriz(void){
4528
4529 /*Constants*/
4530 const int numdof=NDOF2*NUMVERTICES;
4531
4532 /*Intermediaries */
4533 int i,resp,ig;
4534 int *responses=NULL;
4535 int num_responses;
4536 double Jdet;
4537 double obs_velocity_mag,velocity_mag;
4538 double dux,duy;
4539 double epsvel=2.220446049250313e-16;
4540 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4541 double scalex=0,scaley=0,scale=0,S=0;
4542 double vx,vy,vxobs,vyobs,weight;
4543 double xyz_list[NUMVERTICES][3];
4544 double basis[3];
4545 GaussTria* gauss=NULL;
4546
4547 /*Initialize Element vector*/
4548 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
4549
4550 /*Retrieve all inputs and parameters*/
4551 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4552 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
4553 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
4554 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4555 Input* vx_input =inputs->GetInput(VxEnum); _assert_(vx_input);
4556 Input* vy_input =inputs->GetInput(VyEnum); _assert_(vy_input);
4557 Input* vxobs_input =inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4558 Input* vyobs_input =inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4559
4560 /*Get Surface if required by one response*/
4561 for(resp=0;resp<num_responses;resp++){
4562 if(responses[resp]==SurfaceAverageVelMisfitEnum){
4563 inputs->GetInputValue(&S,SurfaceAreaEnum); break;
4564 }
4565 }
4566
4567 /* Start looping on the number of gaussian points: */
4568 gauss=new GaussTria(4);
4569 for (ig=gauss->begin();ig<gauss->end();ig++){
4570
4571 gauss->GaussPoint(ig);
4572
4573 /* Get Jacobian determinant: */
4574 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4575
4576 /*Get all parameters at gaussian point*/
4577 vx_input->GetInputValue(&vx,gauss);
4578 vy_input->GetInputValue(&vy,gauss);
4579 vxobs_input->GetInputValue(&vxobs,gauss);
4580 vyobs_input->GetInputValue(&vyobs,gauss);
4581 GetNodalFunctions(basis, gauss);
4582
4583 /*Loop over all requested responses*/
4584 for(resp=0;resp<num_responses;resp++){
4585
4586 weights_input->GetInputValue(&weight,gauss,resp);
4587
4588 switch(responses[resp]){
4589 case SurfaceAbsVelMisfitEnum:
4590 /*
4591 * 1 [ 2 2 ]
4592 * J = --- | (u - u ) + (v - v ) |
4593 * 2 [ obs obs ]
4594 *
4595 * dJ
4596 * DU = - -- = (u - u )
4597 * du obs
4598 */
4599 for (i=0;i<NUMVERTICES;i++){
4600 dux=vxobs-vx;
4601 duy=vyobs-vy;
4602 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4603 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4604 }
4605 break;
4606 case SurfaceRelVelMisfitEnum:
4607 /*
4608 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
4609 * J = --- | ------------- (u - u ) + ------------- (v - v ) |
4610 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
4611 * obs obs
4612 *
4613 * dJ \bar{v}^2
4614 * DU = - -- = ------------- (u - u )
4615 * du (u + eps)^2 obs
4616 * obs
4617 */
4618 for (i=0;i<NUMVERTICES;i++){
4619 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
4620 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
4621 dux=scalex*(vxobs-vx);
4622 duy=scaley*(vyobs-vy);
4623 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4624 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4625 }
4626 break;
4627 case SurfaceLogVelMisfitEnum:
4628 /*
4629 * [ vel + eps ] 2
4630 * J = 4 \bar{v}^2 | log ( ----------- ) |
4631 * [ vel + eps ]
4632 * obs
4633 *
4634 * dJ 2 * log(...)
4635 * DU = - -- = - 4 \bar{v}^2 ------------- u
4636 * du vel^2 + eps
4637 *
4638 */
4639 for (i=0;i<NUMVERTICES;i++){
4640 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel;
4641 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
4642 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
4643 dux=scale*vx;
4644 duy=scale*vy;
4645 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4646 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4647 }
4648 break;
4649 case SurfaceAverageVelMisfitEnum:
4650 /*
4651 * 1 2 2
4652 * J = --- sqrt( (u - u ) + (v - v ) )
4653 * S obs obs
4654 *
4655 * dJ 1 1
4656 * DU = - -- = - --- ----------- * 2 (u - u )
4657 * du S 2 sqrt(...) obs
4658 */
4659 for (i=0;i<NUMVERTICES;i++){
4660 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
4661 dux=scale*(vxobs-vx);
4662 duy=scale*(vyobs-vy);
4663 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4664 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4665 }
4666 break;
4667 case SurfaceLogVxVyMisfitEnum:
4668 /*
4669 * 1 [ |u| + eps 2 |v| + eps 2 ]
4670 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
4671 * 2 [ |u |+ eps |v |+ eps ]
4672 * obs obs
4673 * dJ 1 u 1
4674 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------
4675 * du |u| + eps |u| u + eps
4676 */
4677 for (i=0;i<NUMVERTICES;i++){
4678 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
4679 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
4680 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4681 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4682 }
4683 break;
4684 case DragCoefficientAbsGradientEnum:
4685 /*Nothing in P vector*/
4686 break;
4687 case ThicknessAbsGradientEnum:
4688 /*Nothing in P vector*/
4689 break;
4690 case RheologyBbarAbsGradientEnum:
4691 /*Nothing in P vector*/
4692 break;
4693 default:
4694 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
4695 }
4696 }
4697 }
4698
4699 /*Clean up and return*/
4700 delete gauss;
4701 xfree((void**)&responses);
4702 return pe;
4703}
4704/*}}}*/
4705/*FUNCTION Tria::CreatePVectorAdjointStokes{{{1*/
4706ElementVector* Tria::CreatePVectorAdjointStokes(void){
4707
4708 /*Intermediaries */
4709 int i,resp,ig;
4710 int *responses=NULL;
4711 int num_responses;
4712 double Jdet;
4713 double obs_velocity_mag,velocity_mag;
4714 double dux,duy;
4715 double epsvel=2.220446049250313e-16;
4716 double meanvel=3.170979198376458e-05; /*1000 m/yr*/
4717 double scalex=0,scaley=0,scale=0,S=0;
4718 double vx,vy,vxobs,vyobs,weight;
4719 double xyz_list[NUMVERTICES][3];
4720 double basis[3];
4721 GaussTria* gauss=NULL;
4722
4723 /*Initialize Element vector*/
4724 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
4725
4726 /*Retrieve all inputs and parameters*/
4727 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4728 this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
4729 this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
4730 Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4731 Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
4732 Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
4733 Input* vxobs_input = inputs->GetInput(InversionVxObsEnum); _assert_(vxobs_input);
4734 Input* vyobs_input = inputs->GetInput(InversionVyObsEnum); _assert_(vyobs_input);
4735
4736 /*Get Surface if required by one response*/
4737 for(resp=0;resp<num_responses;resp++){
4738 if(responses[resp]==SurfaceAverageVelMisfitEnum){
4739 inputs->GetInputValue(&S,SurfaceAreaEnum); break;
4740 }
4741 }
4742
4743 /* Start looping on the number of gaussian points: */
4744 gauss=new GaussTria(4);
4745 for (ig=gauss->begin();ig<gauss->end();ig++){
4746
4747 gauss->GaussPoint(ig);
4748
4749 /* Get Jacobian determinant: */
4750 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4751
4752 /*Get all parameters at gaussian point*/
4753 vx_input->GetInputValue(&vx,gauss);
4754 vy_input->GetInputValue(&vy,gauss);
4755 vxobs_input->GetInputValue(&vxobs,gauss);
4756 vyobs_input->GetInputValue(&vyobs,gauss);
4757 GetNodalFunctions(basis, gauss);
4758
4759 /*Loop over all requested responses*/
4760 for(resp=0;resp<num_responses;resp++){
4761
4762 weights_input->GetInputValue(&weight,gauss,resp);
4763
4764 switch(responses[resp]){
4765
4766 case SurfaceAbsVelMisfitEnum:
4767 /*
4768 * 1 [ 2 2 ]
4769 * J = --- | (u - u ) + (v - v ) |
4770 * 2 [ obs obs ]
4771 *
4772 * dJ
4773 * DU = - -- = (u - u )
4774 * du obs
4775 */
4776 for (i=0;i<NUMVERTICES;i++){
4777 dux=vxobs-vx;
4778 duy=vyobs-vy;
4779 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4780 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4781 }
4782 break;
4783 case SurfaceRelVelMisfitEnum:
4784 /*
4785 * 1 [ \bar{v}^2 2 \bar{v}^2 2 ]
4786 * J = --- | ------------- (u - u ) + ------------- (v - v ) |
4787 * 2 [ (u + eps)^2 obs (v + eps)^2 obs ]
4788 * obs obs
4789 *
4790 * dJ \bar{v}^2
4791 * DU = - -- = ------------- (u - u )
4792 * du (u + eps)^2 obs
4793 * obs
4794 */
4795 for (i=0;i<NUMVERTICES;i++){
4796 scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
4797 scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
4798 dux=scalex*(vxobs-vx);
4799 duy=scaley*(vyobs-vy);
4800 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4801 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4802 }
4803 break;
4804 case SurfaceLogVelMisfitEnum:
4805 /*
4806 * [ vel + eps ] 2
4807 * J = 4 \bar{v}^2 | log ( ----------- ) |
4808 * [ vel + eps ]
4809 * obs
4810 *
4811 * dJ 2 * log(...)
4812 * DU = - -- = - 4 \bar{v}^2 ------------- u
4813 * du vel^2 + eps
4814 *
4815 */
4816 for (i=0;i<NUMVERTICES;i++){
4817 velocity_mag =sqrt(pow(vx, 2.)+pow(vy, 2.))+epsvel;
4818 obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
4819 scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
4820 dux=scale*vx;
4821 duy=scale*vy;
4822 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4823 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4824 }
4825 break;
4826 case SurfaceAverageVelMisfitEnum:
4827 /*
4828 * 1 2 2
4829 * J = --- sqrt( (u - u ) + (v - v ) )
4830 * S obs obs
4831 *
4832 * dJ 1 1
4833 * DU = - -- = - --- ----------- * 2 (u - u )
4834 * du S 2 sqrt(...) obs
4835 */
4836 for (i=0;i<NUMVERTICES;i++){
4837 scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
4838 dux=scale*(vxobs-vx);
4839 duy=scale*(vyobs-vy);
4840 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4841 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4842 }
4843 break;
4844 case SurfaceLogVxVyMisfitEnum:
4845 /*
4846 * 1 [ |u| + eps 2 |v| + eps 2 ]
4847 * J = --- \bar{v}^2 | log ( ----------- ) + log ( ----------- ) |
4848 * 2 [ |u |+ eps |v |+ eps ]
4849 * obs obs
4850 * dJ 1 u 1
4851 * DU = - -- = - \bar{v}^2 log(u...) --------- ---- ~ - \bar{v}^2 log(u...) ------
4852 * du |u| + eps |u| u + eps
4853 */
4854 for (i=0;i<NUMVERTICES;i++){
4855 dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
4856 duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
4857 pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i];
4858 pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i];
4859 }
4860 break;
4861 case DragCoefficientAbsGradientEnum:
4862 /*Nothing in P vector*/
4863 break;
4864 case ThicknessAbsGradientEnum:
4865 /*Nothing in P vector*/
4866 break;
4867 case RheologyBbarAbsGradientEnum:
4868 /*Nothing in P vector*/
4869 break;
4870 default:
4871 _error_("response %s not supported yet",EnumToStringx(responses[resp]));
4872 }
4873 }
4874 }
4875
4876 /*Clean up and return*/
4877 delete gauss;
4878 xfree((void**)&responses);
4879 return pe;
4880}
4881/*}}}*/
4882/*FUNCTION Tria::DragCoefficientAbsGradient{{{1*/
4883double Tria::DragCoefficientAbsGradient(bool process_units,int weight_index){
4884
4885 /* Intermediaries */
4886 int ig;
4887 double Jelem = 0;
4888 double weight;
4889 double Jdet;
4890 double xyz_list[NUMVERTICES][3];
4891 double dp[NDOF2];
4892 GaussTria *gauss = NULL;
4893
4894 /*retrieve parameters and inputs*/
4895
4896 /*If on water, return 0: */
4897 if(IsOnWater()) return 0;
4898
4899 /*Retrieve all inputs we will be needing: */
4900 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4901 Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
4902 Input* drag_input =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input);
4903
4904 /* Start looping on the number of gaussian points: */
4905 gauss=new GaussTria(2);
4906 for (ig=gauss->begin();ig<gauss->end();ig++){
4907
4908 gauss->GaussPoint(ig);
4909
4910 /* Get Jacobian determinant: */
4911 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4912
4913 /*Get all parameters at gaussian point*/
4914 weights_input->GetInputValue(&weight,gauss,weight_index);
4915 drag_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
4916
4917 /*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */
4918 Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
4919 }
4920
4921 /*Clean up and return*/
4922 delete gauss;
4923 return Jelem;
4924}
4925/*}}}*/
4926/*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{1*/
4927ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){
4928
4929 ElementMatrix* Ke=NULL;
4930
4931 /*Get Element Matrix of the forward model*/
4932 switch(GetElementType()){
4933 case P1Enum:
4934 Ke=CreateKMatrixBalancethickness_CG();
4935 break;
4936 case P1DGEnum:
4937 Ke=CreateKMatrixBalancethickness_DG();
4938 break;
4939 default:
4940 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
4941 }
4942
4943 /*Transpose and return Ke*/
4944 Ke->Transpose();
4945 return Ke;
4946}
4947/*}}}*/
4948/*FUNCTION Tria::CreateKMatrixAdjointMacAyeal{{{1*/
4949ElementMatrix* Tria::CreateKMatrixAdjointMacAyeal(void){
4950
4951 /*Constants*/
4952 const int numdof=NDOF2*NUMVERTICES;
4953
4954 /*Intermediaries */
4955 int i,j,ig;
4956 bool incomplete_adjoint;
4957 double xyz_list[NUMVERTICES][3];
4958 double Jdet,thickness;
4959 double eps1dotdphii,eps1dotdphij;
4960 double eps2dotdphii,eps2dotdphij;
4961 double mu_prime;
4962 double epsilon[3];/* epsilon=[exx,eyy,exy];*/
4963 double eps1[2],eps2[2];
4964 double phi[NUMVERTICES];
4965 double dphi[2][NUMVERTICES];
4966 GaussTria *gauss=NULL;
4967
4968 /*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
4969 parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
4970 ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
4971 if(incomplete_adjoint) return Ke;
4972
4973 /*Retrieve all inputs and parameters*/
4974 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
4975 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
4976 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
4977 Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
4978
4979 /* Start looping on the number of gaussian points: */
4980 gauss=new GaussTria(2);
4981 for (ig=gauss->begin();ig<gauss->end();ig++){
4982
4983 gauss->GaussPoint(ig);
4984
4985 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
4986 GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);
4987
4988 thickness_input->GetInputValue(&thickness, gauss);
4989 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
4990 matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
4991 eps1[0]=2*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
4992 eps1[1]=epsilon[2]; eps2[1]=epsilon[0]+2*epsilon[1];
4993
4994 for(i=0;i<3;i++){
4995 for(j=0;j<3;j++){
4996 eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
4997 eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
4998 eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
4999 eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];
5000
5001 Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
5002 Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
5003 Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
5004 Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
5005 }
5006 }
5007 }
5008
5009 /*Transform Coordinate System*/
5010 TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
5011
5012 /*Clean up and return*/
5013 delete gauss;
5014 //Ke->Transpose();
5015 return Ke;
5016}
5017/*}}}*/
5018/*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/
5019void Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){
5020
5021 const int numdof=NDOF2*NUMVERTICES;
5022
5023 int i;
5024 int* doflist=NULL;
5025 double values[numdof];
5026 double lambdax[NUMVERTICES];
5027 double lambday[NUMVERTICES];
5028
5029 /*Get dof list: */
5030 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5031
5032 /*Use the dof list to index into the solution vector: */
5033 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
5034
5035 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
5036 for(i=0;i<NUMVERTICES;i++){
5037 lambdax[i]=values[i*NDOF2+0];
5038 lambday[i]=values[i*NDOF2+1];
5039
5040 /*Check solution*/
5041 if(isnan(lambdax[i])) _error_("NaN found in solution vector");
5042 if(isnan(lambday[i])) _error_("NaN found in solution vector");
5043 }
5044
5045 /*Add vx and vy as inputs to the tria element: */
5046 this->inputs->AddInput(new TriaP1Input(AdjointxEnum,lambdax));
5047 this->inputs->AddInput(new TriaP1Input(AdjointyEnum,lambday));
5048
5049 /*Free ressources:*/
5050 xfree((void**)&doflist);
5051}
5052/*}}}*/
5053
5054/*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{1*/
5055void Tria::InputUpdateFromSolutionAdjointBalancethickness(double* solution){
5056
5057 const int numdof=NDOF1*NUMVERTICES;
5058
5059 int i;
5060 int* doflist=NULL;
5061 double values[numdof];
5062 double lambda[NUMVERTICES];
5063
5064 /*Get dof list: */
5065 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5066
5067 /*Use the dof list to index into the solution vector: */
5068 for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
5069
5070 /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
5071 for(i=0;i<numdof;i++){
5072 lambda[i]=values[i];
5073 if(isnan(lambda[i])) _error_("NaN found in solution vector");
5074 }
5075
5076 /*Add vx and vy as inputs to the tria element: */
5077 this->inputs->AddInput(new TriaP1Input(AdjointEnum,lambda));
5078
5079 /*Free ressources:*/
5080 xfree((void**)&doflist);
5081}
5082/*}}}*/
5083/*FUNCTION Tria::GetVectorFromControlInputs{{{1*/
5084void Tria::GetVectorFromControlInputs(Vector* vector,int control_enum,int control_index,const char* data){
5085
5086 int doflist1[NUMVERTICES];
5087 Input *input=NULL;
5088
5089 /*Get out if this is not an element input*/
5090 if(!IsInput(control_enum)) return;
5091
5092 /*Prepare index list*/
5093 GradientIndexing(&doflist1[0],control_index);
5094
5095 /*Get input (either in element or material)*/
5096 if(control_enum==MaterialsRheologyBbarEnum){
5097 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
5098 }
5099 else{
5100 input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
5101 }
5102
5103 /*Check that it is a ControlInput*/
5104 if (input->ObjectEnum()!=ControlInputEnum){
5105 _error_("input %s is not a ControlInput",EnumToStringx(control_enum));
5106 }
5107
5108 ((ControlInput*)input)->GetVectorFromInputs(vector,&doflist1[0],data);
5109}
5110/*}}}*/
5111/*FUNCTION Tria::SetControlInputsFromVector{{{1*/
5112void Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){
5113
5114 double values[NUMVERTICES];
5115 int doflist1[NUMVERTICES];
5116 Input *input = NULL;
5117 Input *new_input = NULL;
5118
5119 /*Get out if this is not an element input*/
5120 if(!IsInput(control_enum)) return;
5121
5122 /*Prepare index list*/
5123 GradientIndexing(&doflist1[0],control_index);
5124
5125 /*Get values on vertices*/
5126 for (int i=0;i<NUMVERTICES;i++){
5127 values[i]=vector[doflist1[i]];
5128 }
5129 new_input = new TriaP1Input(control_enum,values);
5130
5131 if(control_enum==MaterialsRheologyBbarEnum){
5132 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
5133 }
5134 else{
5135 input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
5136 }
5137
5138 if (input->ObjectEnum()!=ControlInputEnum){
5139 _error_("input %s is not a ControlInput",EnumToStringx(control_enum));
5140 }
5141
5142 ((ControlInput*)input)->SetInput(new_input);
5143}
5144/*}}}*/
5145#endif
5146
5147#ifdef _HAVE_HYDROLOGY_
5148/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{1*/
5149void Tria::CreateHydrologyWaterVelocityInput(void){
5150
5151 /*material parameters: */
5152 double mu_water;
5153 double VelocityFactor; // This factor represents the number 12 in laminar flow velocity which can vary by differnt hydrology.CR
5154 double n_man,CR;
5155 double w;
5156 double rho_ice, rho_water, g;
5157 double dsdx,dsdy,dbdx,dbdy;
5158 double vx[NUMVERTICES];
5159 double vy[NUMVERTICES];
5160 GaussTria *gauss = NULL;
5161
5162 /*Retrieve all inputs and parameters*/
5163 rho_ice=matpar->GetRhoIce();
5164 rho_water=matpar->GetRhoWater();
5165 g=matpar->GetG();
5166 CR=matpar->GetHydrologyCR(); // To have Lebrocq equavalent equation: CR=0.01,n_man=0.02
5167 n_man=matpar->GetHydrologyN();
5168 mu_water=matpar->GetMuWater();
5169 Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
5170 Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
5171 Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(bedslopex_input);
5172 Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(bedslopey_input);
5173 Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
5174
5175 /* compute VelocityFactor */
5176 VelocityFactor= n_man*pow(CR,2)*rho_water*g/mu_water;
5177
5178 gauss=new GaussTria();
5179 for (int iv=0;iv<NUMVERTICES;iv++){
5180 gauss->GaussVertex(iv);
5181 surfaceslopex_input->GetInputValue(&dsdx,gauss);
5182 surfaceslopey_input->GetInputValue(&dsdy,gauss);
5183 bedslopex_input->GetInputValue(&dbdx,gauss);
5184 bedslopey_input->GetInputValue(&dbdy,gauss);
5185 watercolumn_input->GetInputValue(&w,gauss);
5186
5187 /* Water velocity x and y components */
5188 // vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
5189 // vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
5190
5191 vx[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
5192 vy[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
5193 }
5194
5195 /*clean-up*/
5196 delete gauss;
5197
5198 /*Add to inputs*/
5199 this->inputs->AddInput(new TriaP1Input(HydrologyWaterVxEnum,vx));
5200 this->inputs->AddInput(new TriaP1Input(HydrologyWaterVyEnum,vy));
5201}
5202/*}}}*/
5203/*FUNCTION Tria::CreateKMatrixHydrology{{{1*/
5204ElementMatrix* Tria::CreateKMatrixHydrology(void){
5205
5206 /*Constants*/
5207 const int numdof=NDOF1*NUMVERTICES;
5208
5209 /*Intermediaries */
5210 double diffusivity;
5211 int i,j,ig;
5212 double Jdettria,DL_scalar,dt,h;
5213 double vx,vy,vel,dvxdx,dvydy;
5214 double dvx[2],dvy[2];
5215 double v_gauss[2]={0.0};
5216 double xyz_list[NUMVERTICES][3];
5217 double L[NUMVERTICES];
5218 double B[2][NUMVERTICES];
5219 double Bprime[2][NUMVERTICES];
5220 double K[2][2] ={0.0};
5221 double KDL[2][2] ={0.0};
5222 double DL[2][2] ={0.0};
5223 double DLprime[2][2] ={0.0};
5224 GaussTria *gauss=NULL;
5225
5226 /*Skip if water or ice shelf element*/
5227 if(IsOnWater() | IsFloating()) return NULL;
5228
5229 /*Initialize Element matrix*/
5230 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
5231
5232 /*Create water velocity vx and vy from current inputs*/
5233 CreateHydrologyWaterVelocityInput();
5234
5235 /*Retrieve all inputs and parameters*/
5236 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5237 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
5238 this->parameters->FindParam(&diffusivity,HydrologyStabilizationEnum);
5239 Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
5240 Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
5241 h=sqrt(2*this->GetArea());
5242
5243 /* Start looping on the number of gaussian points: */
5244 gauss=new GaussTria(2);
5245 for (ig=gauss->begin();ig<gauss->end();ig++){
5246
5247 gauss->GaussPoint(ig);
5248
5249 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5250 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
5251
5252 vx_input->GetInputValue(&vx,gauss);
5253 vy_input->GetInputValue(&vy,gauss);
5254 vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
5255 vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
5256
5257 DL_scalar=gauss->weight*Jdettria;
5258
5259 TripleMultiply( &L[0],1,numdof,1,
5260 &DL_scalar,1,1,0,
5261 &L[0],1,numdof,0,
5262 &Ke->values[0],1);
5263
5264 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
5265 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
5266
5267 dvxdx=dvx[0];
5268 dvydy=dvy[1];
5269 DL_scalar=dt*gauss->weight*Jdettria;
5270
5271 DL[0][0]=DL_scalar*dvxdx;
5272 DL[1][1]=DL_scalar*dvydy;
5273 DLprime[0][0]=DL_scalar*vx;
5274 DLprime[1][1]=DL_scalar*vy;
5275
5276 TripleMultiply( &B[0][0],2,numdof,1,
5277 &DL[0][0],2,2,0,
5278 &B[0][0],2,numdof,0,
5279 &Ke->values[0],1);
5280
5281 TripleMultiply( &B[0][0],2,numdof,1,
5282 &DLprime[0][0],2,2,0,
5283 &Bprime[0][0],2,numdof,0,
5284 &Ke->values[0],1);
5285
5286 /*Artificial diffusivity*/
5287 vel=sqrt(pow(vx,2.)+pow(vy,2.));
5288 K[0][0]=diffusivity*h/(2*vel)*vx*vx;
5289 K[1][0]=diffusivity*h/(2*vel)*vy*vx;
5290 K[0][1]=diffusivity*h/(2*vel)*vx*vy;
5291 K[1][1]=diffusivity*h/(2*vel)*vy*vy;
5292 KDL[0][0]=DL_scalar*K[0][0];
5293 KDL[1][0]=DL_scalar*K[1][0];
5294 KDL[0][1]=DL_scalar*K[0][1];
5295 KDL[1][1]=DL_scalar*K[1][1];
5296
5297 TripleMultiply( &Bprime[0][0],2,numdof,1,
5298 &KDL[0][0],2,2,0,
5299 &Bprime[0][0],2,numdof,0,
5300 &Ke->values[0],1);
5301 }
5302
5303 /*Clean up and return*/
5304 delete gauss;
5305 return Ke;
5306}
5307/*}}}*/
5308/*FUNCTION Tria::CreatePVectorHydrology {{{1*/
5309ElementVector* Tria::CreatePVectorHydrology(void){
5310
5311 /*Constants*/
5312 const int numdof=NDOF1*NUMVERTICES;
5313
5314 /*Intermediaries */
5315 int i,j,ig;
5316 double Jdettria,dt;
5317 double basal_melting_g;
5318 double old_watercolumn_g;
5319 double xyz_list[NUMVERTICES][3];
5320 double basis[numdof];
5321 GaussTria* gauss=NULL;
5322
5323 /*Skip if water or ice shelf element*/
5324 if(IsOnWater() | IsFloating()) return NULL;
5325
5326 /*Initialize Element vector*/
5327 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
5328
5329 /*Retrieve all inputs and parameters*/
5330 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5331 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
5332 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
5333 Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum); _assert_(old_watercolumn_input);
5334
5335 /*Initialize basal_melting_correction_g to 0, do not forget!:*/
5336 /* Start looping on the number of gaussian points: */
5337 gauss=new GaussTria(2);
5338 for(ig=gauss->begin();ig<gauss->end();ig++){
5339
5340 gauss->GaussPoint(ig);
5341
5342 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5343 GetNodalFunctions(basis, gauss);
5344
5345 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
5346 old_watercolumn_input->GetInputValue(&old_watercolumn_g,gauss);
5347
5348 if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
5349 else for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
5350 }
5351
5352 /*Clean up and return*/
5353 delete gauss;
5354 return pe;
5355}
5356/*}}}*/
5357/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
5358void Tria::GetSolutionFromInputsHydrology(Vector* solution){
5359
5360 const int numdof=NDOF1*NUMVERTICES;
5361
5362 int i;
5363 int* doflist=NULL;
5364 double watercolumn;
5365 double values[numdof];
5366 GaussTria* gauss=NULL;
5367
5368 /*Get dof list: */
5369 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5370
5371 /*Get inputs*/
5372 Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
5373
5374 /*Ok, we have watercolumn values, fill in watercolumn array: */
5375 /*P1 element only for now*/
5376 gauss=new GaussTria();
5377 for(i=0;i<NUMVERTICES;i++){
5378
5379 gauss->GaussVertex(i);
5380
5381 /*Recover watercolumn*/
5382 watercolumn_input->GetInputValue(&watercolumn,gauss);
5383 values[i]=watercolumn;
5384 }
5385
5386 solution->SetValues(numdof,doflist,values,INS_VAL);
5387
5388 /*Free ressources:*/
5389 delete gauss;
5390 xfree((void**)&doflist);
5391}
5392/*}}}*/
5393/*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{1*/
5394void Tria::InputUpdateFromSolutionHydrology(double* solution){
5395
5396 /*Intermediaries*/
5397 const int numdof = NDOF1*NUMVERTICES;
5398
5399 int i;
5400 int* doflist=NULL;
5401 double values[numdof];
5402
5403 /*Get dof list: */
5404 GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
5405
5406 /*Use the dof list to index into the solution vector: */
5407 for(i=0;i<numdof;i++){
5408 values[i]=solution[doflist[i]];
5409 if(isnan(values[i])) _error_("NaN found in solution vector");
5410 if (values[i]<pow((double)10,(double)-10))values[i]=pow((double)10,(double)-10); //correcting the water column to positive values
5411
5412 }
5413
5414 /*Add input to the element: */
5415 this->inputs->AddInput(new TriaP1Input(WatercolumnEnum,values));
5416
5417 /*Free ressources:*/
5418 xfree((void**)&doflist);
5419}
5420/*}}}*/
5421#endif
5422
5423#ifdef _HAVE_DAKOTA_
5424/*FUNCTION Tria::InputUpdateFromVectorDakota(double* vector, int name, int type);{{{1*/
5425void Tria::InputUpdateFromVectorDakota(double* vector, int name, int type){
5426
5427 int i,j;
5428
5429 /*Check that name is an element input*/
5430 if (!IsInput(name)) return;
5431
5432 switch(type){
5433
5434 case VertexEnum:
5435
5436 /*New TriaP1Input*/
5437 double values[3];
5438
5439 /*Get values on the 3 vertices*/
5440 for (i=0;i<3;i++){
5441 values[i]=vector[this->nodes[i]->GetSidList()]; //careful, vector of values here is not parallel distributed, but serial distributed (from a serial Dakota core!)
5442 }
5443
5444 /*Branch on the specified type of update: */
5445 switch(name){
5446 case ThicknessEnum:
5447 /*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium {{{2*/
5448 double thickness[3];
5449 double thickness_init[3];
5450 double hydrostatic_ratio[3];
5451 double surface[3];
5452 double bed[3];
5453
5454 /*retrieve inputs: */
5455 GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
5456 GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
5457 GetInputListOnVertices(&bed[0],BedEnum);
5458 GetInputListOnVertices(&surface[0],SurfaceEnum);
5459
5460 /*build new thickness: */
5461// for(j=0;j<3;j++)thickness[j]=values[j];
5462
5463 /*build new bed and surface: */
5464 if (this->IsFloating()){
5465 /*hydrostatic equilibrium: */
5466 double rho_ice,rho_water,di;
5467 rho_ice=this->matpar->GetRhoIce();
5468 rho_water=this->matpar->GetRhoWater();
5469
5470 di=rho_ice/rho_water;
5471
5472 /*build new thickness: */
5473 for (j=0; j<3; j++) {
5474 /* for observed/interpolated/hydrostatic thickness, remove scaling from any hydrostatic thickness */
5475 if (hydrostatic_ratio[j] >= 0.)
5476 thickness[j]=values[j]-(values[j]/thickness_init[j]-1.)*hydrostatic_ratio[j]*surface[j]/(1.-di);
5477 /* for minimum thickness, don't scale */
5478 else
5479 thickness[j]=thickness_init[j];
5480
5481 /* check the computed thickness and update bed */
5482 if (thickness[j] < 0.)
5483 thickness[j]=1./(1.-di);
5484 bed[j]=surface[j]-thickness[j];
5485 }
5486
5487// for(j=0;j<3;j++){
5488// surface[j]=(1-di)*thickness[j];
5489// bed[j]=-di*thickness[j];
5490// }
5491 }
5492 else{
5493 /*build new thickness: */
5494 for (j=0; j<3; j++) {
5495 /* for observed thickness, use scaled value */
5496 if (hydrostatic_ratio[j] >= 0.)
5497 thickness[j]=values[j];
5498 /* for minimum thickness, don't scale */
5499 else
5500 thickness[j]=thickness_init[j];
5501 }
5502
5503 /*update bed on grounded ice: */
5504// for(j=0;j<3;j++)surface[j]=bed[j]+thickness[j];
5505 for(j=0;j<3;j++)bed[j]=surface[j]-thickness[j];
5506 }
5507
5508 /*Add new inputs: */
5509 this->inputs->AddInput(new TriaP1Input(ThicknessEnum,thickness));
5510 this->inputs->AddInput(new TriaP1Input(BedEnum,bed));
5511 this->inputs->AddInput(new TriaP1Input(SurfaceEnum,surface));
5512
5513 /*}}}*/
5514 break;
5515 default:
5516 this->inputs->AddInput(new TriaP1Input(name,values));
5517 }
5518 break;
5519
5520 default:
5521 _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
5522 }
5523
5524}
5525/*}}}*/
5526/*FUNCTION Tria::InputUpdateFromVectorDakota(int* vector, int name, int type);{{{1*/
5527void Tria::InputUpdateFromVectorDakota(int* vector, int name, int type){
5528 _error_(" not supported yet!");
5529}
5530/*}}}*/
5531/*FUNCTION Tria::InputUpdateFromVectorDakota(bool* vector, int name, int type);{{{1*/
5532void Tria::InputUpdateFromVectorDakota(bool* vector, int name, int type){
5533 _error_(" not supported yet!");
5534}
5535/*}}}*/
5536/*FUNCTION Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type);{{{1*/
5537void Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){
5538
5539 int i,j,t;
5540 TransientInput* transientinput=NULL;
5541 double values[3];
5542 double time;
5543 int row;
5544 double yts;
5545
5546 /*Check that name is an element input*/
5547 if (!IsInput(name)) return;
5548
5549 switch(type){
5550
5551 case VertexEnum:
5552
5553 /*Create transient input: */
5554
5555 parameters->FindParam(&yts,ConstantsYtsEnum);
5556 for(t=0;t<ncols;t++){ //ncols is the number of times
5557
5558 /*create input values: */
5559 for(i=0;i<3;i++){
5560 row=this->nodes[i]->GetSidList();
5561 values[i]=(double)matrix[ncols*row+t];
5562 }
5563
5564 /*time? :*/
5565 time=(double)matrix[(nrows-1)*ncols+t]*yts;
5566
5567 if(t==0) transientinput=new TransientInput(name);
5568 transientinput->AddTimeInput(new TriaP1Input(name,values),time);
5569 transientinput->Configure(parameters);
5570 }
5571 this->inputs->AddInput(transientinput);
5572 break;
5573
5574 default:
5575 _error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
5576 }
5577
5578}
5579/*}}}*/
5580#endif
5581
5582#ifdef _HAVE_BALANCED_
5583/*FUNCTION Tria::CreateKMatrixBalancethickness {{{1*/
5584ElementMatrix* Tria::CreateKMatrixBalancethickness(void){
5585
5586 switch(GetElementType()){
5587 case P1Enum:
5588 return CreateKMatrixBalancethickness_CG();
5589 case P1DGEnum:
5590 return CreateKMatrixBalancethickness_DG();
5591 default:
5592 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
5593 }
5594
5595}
5596/*}}}*/
5597/*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{1*/
5598ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){
5599
5600 /*Constants*/
5601 const int numdof=NDOF1*NUMVERTICES;
5602
5603 /*Intermediaries */
5604 int stabilization;
5605 int i,j,ig,dim;
5606 double Jdettria,vx,vy,dvxdx,dvydy,vel,h;
5607 double dvx[2],dvy[2];
5608 double xyz_list[NUMVERTICES][3];
5609 double L[NUMVERTICES];
5610 double B[2][NUMVERTICES];
5611 double Bprime[2][NUMVERTICES];
5612 double K[2][2] = {0.0};
5613 double KDL[2][2] = {0.0};
5614 double DL[2][2] = {0.0};
5615 double DLprime[2][2] = {0.0};
5616 double DL_scalar;
5617 GaussTria *gauss = NULL;
5618
5619 /*Initialize Element matrix*/
5620 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
5621
5622 /*Retrieve all Inputs and parameters: */
5623 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5624 this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum);
5625 this->parameters->FindParam(&dim,MeshDimensionEnum);
5626 Input* vxaverage_input=NULL;
5627 Input* vyaverage_input=NULL;
5628 if(dim==2){
5629 vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
5630 vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
5631 }
5632 else{
5633 vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
5634 vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
5635 }
5636 h=sqrt(2*this->GetArea());
5637
5638 /*Start looping on the number of gaussian points:*/
5639 gauss=new GaussTria(2);
5640 for (ig=gauss->begin();ig<gauss->end();ig++){
5641
5642 gauss->GaussPoint(ig);
5643
5644 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5645 GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
5646 GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
5647
5648 vxaverage_input->GetInputValue(&vx,gauss);
5649 vyaverage_input->GetInputValue(&vy,gauss);
5650 vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
5651 vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
5652
5653 dvxdx=dvx[0];
5654 dvydy=dvy[1];
5655 DL_scalar=gauss->weight*Jdettria;
5656
5657 DL[0][0]=DL_scalar*dvxdx;
5658 DL[1][1]=DL_scalar*dvydy;
5659
5660 DLprime[0][0]=DL_scalar*vx;
5661 DLprime[1][1]=DL_scalar*vy;
5662
5663 TripleMultiply( &B[0][0],2,numdof,1,
5664 &DL[0][0],2,2,0,
5665 &B[0][0],2,numdof,0,
5666 &Ke->values[0],1);
5667
5668 TripleMultiply( &B[0][0],2,numdof,1,
5669 &DLprime[0][0],2,2,0,
5670 &Bprime[0][0],2,numdof,0,
5671 &Ke->values[0],1);
5672
5673 if(stabilization==1){
5674 /*Streamline upwinding*/
5675 vel=sqrt(pow(vx,2.)+pow(vy,2.));
5676 K[0][0]=h/(2*vel)*vx*vx;
5677 K[1][0]=h/(2*vel)*vy*vx;
5678 K[0][1]=h/(2*vel)*vx*vy;
5679 K[1][1]=h/(2*vel)*vy*vy;
5680 }
5681 else if(stabilization==2){
5682 /*MacAyeal*/
5683 vxaverage_input->GetInputAverage(&vx);
5684 vyaverage_input->GetInputAverage(&vy);
5685 K[0][0]=h/2.0*fabs(vx);
5686 K[0][1]=0.;
5687 K[1][0]=0.;
5688 K[1][1]=h/2.0*fabs(vy);
5689 }
5690 if(stabilization==1 || stabilization==2){
5691 KDL[0][0]=DL_scalar*K[0][0];
5692 KDL[1][0]=DL_scalar*K[1][0];
5693 KDL[0][1]=DL_scalar*K[0][1];
5694 KDL[1][1]=DL_scalar*K[1][1];
5695 TripleMultiply( &Bprime[0][0],2,numdof,1,
5696 &KDL[0][0],2,2,0,
5697 &Bprime[0][0],2,numdof,0,
5698 &Ke->values[0],1);
5699 }
5700 }
5701
5702 /*Clean up and return*/
5703 delete gauss;
5704 return Ke;
5705}
5706/*}}}*/
5707/*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{1*/
5708ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){
5709
5710 /*Constants*/
5711 const int numdof=NDOF1*NUMVERTICES;
5712
5713 /*Intermediaries*/
5714 int i,j,ig,dim;
5715 double vx,vy,Jdettria;
5716 double xyz_list[NUMVERTICES][3];
5717 double B[2][NUMVERTICES];
5718 double Bprime[2][NUMVERTICES];
5719 double DL[2][2]={0.0};
5720 double DL_scalar;
5721 GaussTria *gauss=NULL;
5722
5723 /*Initialize Element matrix*/
5724 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
5725
5726 /*Retrieve all inputs and parameters*/
5727 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5728 this->parameters->FindParam(&dim,MeshDimensionEnum);
5729 Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
5730 Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
5731
5732 /*Start looping on the number of gaussian points:*/
5733 gauss=new GaussTria(2);
5734 for (ig=gauss->begin();ig<gauss->end();ig++){
5735
5736 gauss->GaussPoint(ig);
5737
5738 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5739 /*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
5740 GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
5741 GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);
5742
5743 vx_input->GetInputValue(&vx,gauss);
5744 vy_input->GetInputValue(&vy,gauss);
5745
5746 DL_scalar=-gauss->weight*Jdettria;
5747 DL[0][0]=DL_scalar*vx;
5748 DL[1][1]=DL_scalar*vy;
5749
5750 TripleMultiply( &B[0][0],2,numdof,1,
5751 &DL[0][0],2,2,0,
5752 &Bprime[0][0],2,numdof,0,
5753 &Ke->values[0],1);
5754 }
5755
5756 /*Clean up and return*/
5757 delete gauss;
5758 return Ke;
5759}
5760/*}}}*/
5761/*FUNCTION Tria::CreatePVectorBalancethickness{{{1*/
5762ElementVector* Tria::CreatePVectorBalancethickness(void){
5763
5764 switch(GetElementType()){
5765 case P1Enum:
5766 return CreatePVectorBalancethickness_CG();
5767 break;
5768 case P1DGEnum:
5769 return CreatePVectorBalancethickness_DG();
5770 default:
5771 _error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
5772 }
5773}
5774/*}}}*/
5775/*FUNCTION Tria::CreatePVectorBalancethickness_CG{{{1*/
5776ElementVector* Tria::CreatePVectorBalancethickness_CG(void){
5777
5778 /*Constants*/
5779 const int numdof=NDOF1*NUMVERTICES;
5780
5781 /*Intermediaries */
5782 int i,j,ig;
5783 double xyz_list[NUMVERTICES][3];
5784 double dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
5785 double L[NUMVERTICES];
5786 GaussTria* gauss=NULL;
5787
5788 /*Initialize Element vector*/
5789 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
5790
5791 /*Retrieve all inputs and parameters*/
5792 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5793 Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
5794 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
5795 Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
5796
5797 /* Start looping on the number of gaussian points: */
5798 gauss=new GaussTria(2);
5799 for(ig=gauss->begin();ig<gauss->end();ig++){
5800
5801 gauss->GaussPoint(ig);
5802
5803 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
5804 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
5805 dhdt_input->GetInputValue(&dhdt_g,gauss);
5806
5807 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5808 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
5809
5810 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
5811 }
5812
5813 /*Clean up and return*/
5814 delete gauss;
5815 return pe;
5816}
5817/*}}}*/
5818/*FUNCTION Tria::CreatePVectorBalancethickness_DG {{{1*/
5819ElementVector* Tria::CreatePVectorBalancethickness_DG(void){
5820
5821 /*Constants*/
5822 const int numdof=NDOF1*NUMVERTICES;
5823
5824 /*Intermediaries */
5825 int i,j,ig;
5826 double xyz_list[NUMVERTICES][3];
5827 double basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
5828 double L[NUMVERTICES];
5829 GaussTria* gauss=NULL;
5830
5831 /*Initialize Element vector*/
5832 ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
5833
5834 /*Retrieve all inputs and parameters*/
5835 GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
5836 Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
5837 Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
5838 Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
5839
5840 /* Start looping on the number of gaussian points: */
5841 gauss=new GaussTria(2);
5842 for(ig=gauss->begin();ig<gauss->end();ig++){
5843
5844 gauss->GaussPoint(ig);
5845
5846 surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
5847 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
5848 dhdt_input->GetInputValue(&dhdt_g,gauss);
5849
5850 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
5851 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
5852
5853 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
5854 }
5855
5856 /*Clean up and return*/
5857 delete gauss;
5858 return pe;
5859}
5860/*}}}*/
5861#endif
Note: See TracBrowser for help on using the repository browser.