Context Navigation

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

Last change on this file since 12004 was 12004, checked in by cborstad, 13 years ago
merged trunk-jpl into trunk-jpl-damage through revision 11990
File size: 180.2 KB

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

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

Download in other formats: