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

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

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