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

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

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