Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

source: issm/trunk/src/c/classes/objects/Elements/Tria.cpp@ 13395

Last change on this file since 13395 was 13395, checked in by Mathieu Morlighem, 12 years ago
merged trunk-jpl and trunk for revision 13393
File size: 190.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: