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

Last change on this file since 11427 was 11332, checked in by Mathieu Morlighem, 13 years ago

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