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

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

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