source: issm/oecreview/Archive/18296-19100/ISSM-18910-18911.diff@ 19102

Last change on this file since 19102 was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 109.9 KB
RevLine 
[19102]1Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18910)
4+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 18911)
5@@ -154,6 +154,99 @@
6 this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
7 }
8 /*}}}*/
9+void Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
10+
11+ bool already = false;
12+ int i,j;
13+ int partition[NUMVERTICES];
14+ int offsetsid[NUMVERTICES];
15+ int offsetdof[NUMVERTICES];
16+ IssmDouble area;
17+ IssmDouble mean;
18+
19+ /*First, get the area: */
20+ area=this->GetArea();
21+
22+ /*Figure out the average for this element: */
23+ this->GetVerticesSidList(&offsetsid[0]);
24+ this->GetVertexPidList(&offsetdof[0]);
25+ mean=0;
26+ for(i=0;i<NUMVERTICES;i++){
27+ partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
28+ mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
29+ }
30+
31+ /*Add contribution: */
32+ for(i=0;i<NUMVERTICES;i++){
33+ already=false;
34+ for(j=0;j<i;j++){
35+ if (partition[i]==partition[j]){
36+ already=true;
37+ break;
38+ }
39+ }
40+ if(!already){
41+ partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
42+ partition_areas->SetValue(partition[i],area,ADD_VAL);
43+ };
44+ }
45+}
46+/*}}}*/
47+void Tria::CalvingRateLevermann(){/*{{{*/
48+
49+ IssmDouble xyz_list[NUMVERTICES][3];
50+ GaussTria* gauss=NULL;
51+ IssmDouble vx,vy,vel;
52+ IssmDouble strainparallel;
53+ IssmDouble propcoeff;
54+ IssmDouble strainperpendicular;
55+ IssmDouble calvingratex[NUMVERTICES];
56+ IssmDouble calvingratey[NUMVERTICES];
57+ IssmDouble calvingrate[NUMVERTICES];
58+
59+
60+ /* Get node coordinates and dof list: */
61+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
62+
63+ /*Retrieve all inputs and parameters we will need*/
64+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
65+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
66+ Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input);
67+ Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input);
68+ Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input);
69+
70+ /* Start looping on the number of vertices: */
71+ gauss=new GaussTria();
72+ for (int iv=0;iv<NUMVERTICES;iv++){
73+ gauss->GaussVertex(iv);
74+
75+ /* Get the value we need*/
76+ vx_input->GetInputValue(&vx,gauss);
77+ vy_input->GetInputValue(&vy,gauss);
78+ vel=vx*vx+vy*vy;
79+ strainparallel_input->GetInputValue(&strainparallel,gauss);
80+ strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
81+ levermanncoeff_input->GetInputValue(&propcoeff,gauss);
82+
83+ /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
84+ calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
85+ if(calvingrate[iv]<0){
86+ calvingrate[iv]=0;
87+ }
88+ calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
89+ calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
90+ }
91+
92+ /*Add input*/
93+ this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
94+ this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
95+ this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
96+
97+ /*Clean up and return*/
98+ delete gauss;
99+
100+}
101+/*}}}*/
102 IssmDouble Tria::CharacteristicLength(void){/*{{{*/
103
104 return sqrt(2*this->GetArea());
105@@ -163,6 +256,56 @@
106 _error_("Not Implemented yet");
107 }
108 /*}}}*/
109+void Tria::ComputeDeviatoricStressTensor(){/*{{{*/
110+
111+ IssmDouble xyz_list[NUMVERTICES][3];
112+ IssmDouble viscosity;
113+ IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
114+ IssmDouble tau_xx[NUMVERTICES];
115+ IssmDouble tau_yy[NUMVERTICES];
116+ IssmDouble tau_zz[NUMVERTICES]={0,0,0};
117+ IssmDouble tau_xy[NUMVERTICES];
118+ IssmDouble tau_xz[NUMVERTICES]={0,0,0};
119+ IssmDouble tau_yz[NUMVERTICES]={0,0,0};
120+ GaussTria* gauss=NULL;
121+ int domaintype,dim=2;
122+
123+ /* Get node coordinates and dof list: */
124+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
125+
126+ /*Retrieve all inputs we will be needing: */
127+ this->FindParam(&domaintype,DomainTypeEnum);
128+ if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
129+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
130+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
131+
132+ /* Start looping on the number of vertices: */
133+ gauss=new GaussTria();
134+ for (int iv=0;iv<NUMVERTICES;iv++){
135+ gauss->GaussVertex(iv);
136+
137+ /*Compute strain rate and viscosity: */
138+ this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
139+ this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
140+
141+ /*Compute Stress*/
142+ tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
143+ tau_yy[iv]=2*viscosity*epsilon[1];
144+ tau_xy[iv]=2*viscosity*epsilon[2];
145+ }
146+
147+ /*Add Stress tensor components into inputs*/
148+ this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
149+ this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
150+ this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
151+ this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
152+ this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
153+ this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
154+
155+ /*Clean up and return*/
156+ delete gauss;
157+}
158+/*}}}*/
159 void Tria::ComputeSigmaNN(){/*{{{*/
160
161 if(!IsOnBase()){
162@@ -275,203 +418,6 @@
163 delete gauss;
164 }
165 /*}}}*/
166-void Tria::ComputeDeviatoricStressTensor(){/*{{{*/
167-
168- IssmDouble xyz_list[NUMVERTICES][3];
169- IssmDouble viscosity;
170- IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
171- IssmDouble tau_xx[NUMVERTICES];
172- IssmDouble tau_yy[NUMVERTICES];
173- IssmDouble tau_zz[NUMVERTICES]={0,0,0};
174- IssmDouble tau_xy[NUMVERTICES];
175- IssmDouble tau_xz[NUMVERTICES]={0,0,0};
176- IssmDouble tau_yz[NUMVERTICES]={0,0,0};
177- GaussTria* gauss=NULL;
178- int domaintype,dim=2;
179-
180- /* Get node coordinates and dof list: */
181- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
182-
183- /*Retrieve all inputs we will be needing: */
184- this->FindParam(&domaintype,DomainTypeEnum);
185- if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
186- Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
187- Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
188-
189- /* Start looping on the number of vertices: */
190- gauss=new GaussTria();
191- for (int iv=0;iv<NUMVERTICES;iv++){
192- gauss->GaussVertex(iv);
193-
194- /*Compute strain rate and viscosity: */
195- this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
196- this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
197-
198- /*Compute Stress*/
199- tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
200- tau_yy[iv]=2*viscosity*epsilon[1];
201- tau_xy[iv]=2*viscosity*epsilon[2];
202- }
203-
204- /*Add Stress tensor components into inputs*/
205- this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
206- this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
207- this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
208- this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
209- this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
210- this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
211-
212- /*Clean up and return*/
213- delete gauss;
214-}
215-/*}}}*/
216-void Tria::StrainRateparallel(){/*{{{*/
217-
218- IssmDouble *xyz_list = NULL;
219- IssmDouble epsilon[3];
220- GaussTria* gauss=NULL;
221- IssmDouble vx,vy,vel;
222- IssmDouble strainxx;
223- IssmDouble strainxy;
224- IssmDouble strainyy;
225- IssmDouble strainparallel[NUMVERTICES];
226-
227- /* Get node coordinates and dof list: */
228- this->GetVerticesCoordinates(&xyz_list);
229-
230- /*Retrieve all inputs we will need*/
231- Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
232- Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
233-
234- /* Start looping on the number of vertices: */
235- gauss=new GaussTria();
236- for (int iv=0;iv<NUMVERTICES;iv++){
237- gauss->GaussVertex(iv);
238-
239- /* Get the value we need*/
240- vx_input->GetInputValue(&vx,gauss);
241- vy_input->GetInputValue(&vy,gauss);
242- vel=vx*vx+vy*vy;
243-
244- /*Compute strain rate viscosity and pressure: */
245- this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
246- strainxx=epsilon[0];
247- strainyy=epsilon[1];
248- strainxy=epsilon[2];
249-
250- /*strainparallel= Strain rate along the ice flow direction */
251- strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
252- }
253-
254- /*Add input*/
255- this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
256-
257- /*Clean up and return*/
258- delete gauss;
259- xDelete<IssmDouble>(xyz_list);
260-}
261-/*}}}*/
262-void Tria::StrainRateperpendicular(){/*{{{*/
263-
264- IssmDouble *xyz_list = NULL;
265- GaussTria* gauss=NULL;
266- IssmDouble epsilon[3];
267- IssmDouble vx,vy,vel;
268- IssmDouble strainxx;
269- IssmDouble strainxy;
270- IssmDouble strainyy;
271- IssmDouble strainperpendicular[NUMVERTICES];
272-
273- /* Get node coordinates and dof list: */
274- this->GetVerticesCoordinates(&xyz_list);
275-
276- /*Retrieve all inputs we will need*/
277- Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
278- Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
279-
280- /* Start looping on the number of vertices: */
281- gauss=new GaussTria();
282- for (int iv=0;iv<NUMVERTICES;iv++){
283- gauss->GaussVertex(iv);
284-
285- /* Get the value we need*/
286- vx_input->GetInputValue(&vx,gauss);
287- vy_input->GetInputValue(&vy,gauss);
288- vel=vx*vx+vy*vy;
289-
290- /*Compute strain rate viscosity and pressure: */
291- this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
292- strainxx=epsilon[0];
293- strainyy=epsilon[1];
294- strainxy=epsilon[2];
295-
296- /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
297- strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
298- }
299-
300- /*Add input*/
301- this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
302-
303- /*Clean up and return*/
304- delete gauss;
305- xDelete<IssmDouble>(xyz_list);
306-}
307-/*}}}*/
308-void Tria::CalvingRateLevermann(){/*{{{*/
309-
310- IssmDouble xyz_list[NUMVERTICES][3];
311- GaussTria* gauss=NULL;
312- IssmDouble vx,vy,vel;
313- IssmDouble strainparallel;
314- IssmDouble propcoeff;
315- IssmDouble strainperpendicular;
316- IssmDouble calvingratex[NUMVERTICES];
317- IssmDouble calvingratey[NUMVERTICES];
318- IssmDouble calvingrate[NUMVERTICES];
319-
320-
321- /* Get node coordinates and dof list: */
322- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
323-
324- /*Retrieve all inputs and parameters we will need*/
325- Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
326- Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
327- Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum); _assert_(strainparallel_input);
328- Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum); _assert_(strainperpendicular_input);
329- Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum); _assert_(levermanncoeff_input);
330-
331- /* Start looping on the number of vertices: */
332- gauss=new GaussTria();
333- for (int iv=0;iv<NUMVERTICES;iv++){
334- gauss->GaussVertex(iv);
335-
336- /* Get the value we need*/
337- vx_input->GetInputValue(&vx,gauss);
338- vy_input->GetInputValue(&vy,gauss);
339- vel=vx*vx+vy*vy;
340- strainparallel_input->GetInputValue(&strainparallel,gauss);
341- strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
342- levermanncoeff_input->GetInputValue(&propcoeff,gauss);
343-
344- /*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
345- calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
346- if(calvingrate[iv]<0){
347- calvingrate[iv]=0;
348- }
349- calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
350- calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
351- }
352-
353- /*Add input*/
354- this->inputs->AddInput(new TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
355- this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
356- this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
357-
358- /*Clean up and return*/
359- delete gauss;
360-
361-}
362-/*}}}*/
363 void Tria::Configure(Elements* elementsin, Loads* loadsin,Nodes* nodesin,Vertices *verticesin,Materials* materialsin, Parameters* parametersin){/*{{{*/
364
365 /*go into parameters and get the analysis_counter: */
366@@ -507,6 +453,54 @@
367
368 }
369 /*}}}*/
370+void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
371+
372+ int vertexpidlist[NUMVERTICES];
373+ IssmDouble grad_list[NUMVERTICES];
374+ Input* grad_input=NULL;
375+
376+ Input* input=inputs->GetInput(enum_type);
377+ if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
378+ if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
379+
380+ GradientIndexing(&vertexpidlist[0],control_index);
381+ for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
382+ grad_input=new TriaInput(GradientEnum,grad_list,P1Enum);
383+
384+ ((ControlInput*)input)->SetGradient(grad_input);
385+
386+}/*}}}*/
387+void Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
388+
389+ Input* input=inputs->GetInput(control_enum);
390+ if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
391+ if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
392+
393+ int sidlist[NUMVERTICES];
394+ int connectivity[NUMVERTICES];
395+ IssmPDouble values[NUMVERTICES];
396+ IssmPDouble gradients[NUMVERTICES];
397+ IssmDouble value,gradient;
398+
399+ this->GetVerticesConnectivityList(&connectivity[0]);
400+ this->GetVerticesSidList(&sidlist[0]);
401+
402+ GaussTria* gauss=new GaussTria();
403+ for (int iv=0;iv<NUMVERTICES;iv++){
404+ gauss->GaussVertex(iv);
405+
406+ ((ControlInput*)input)->GetInputValue(&value,gauss);
407+ ((ControlInput*)input)->GetGradientValue(&gradient,gauss);
408+
409+ values[iv] = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
410+ gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
411+ }
412+ delete gauss;
413+
414+ vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
415+ vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
416+
417+}/*}}}*/
418 void Tria::Delta18oParameterization(void){/*{{{*/
419
420 int i;
421@@ -577,6 +571,107 @@
422 delete gauss;
423 }
424 /*}}}*/
425+int Tria::EdgeOnBaseIndex(void){/*{{{*/
426+
427+ IssmDouble values[NUMVERTICES];
428+ int indices[3][2] = {{1,2},{2,0},{0,1}};
429+
430+ /*Retrieve all inputs and parameters*/
431+ GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
432+
433+ for(int i=0;i<3;i++){
434+ if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
435+ return i;
436+ }
437+ }
438+
439+ _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
440+ _error_("Could not find 2 vertices on bed");
441+}
442+/*}}}*/
443+void Tria::EdgeOnBaseIndices(int* pindex1,int* pindex2){/*{{{*/
444+
445+ IssmDouble values[NUMVERTICES];
446+ int indices[3][2] = {{1,2},{2,0},{0,1}};
447+
448+ /*Retrieve all inputs and parameters*/
449+ GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
450+
451+ for(int i=0;i<3;i++){
452+ if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
453+ *pindex1 = indices[i][0];
454+ *pindex2 = indices[i][1];
455+ return;
456+ }
457+ }
458+
459+ _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
460+ _error_("Could not find 2 vertices on bed");
461+}
462+/*}}}*/
463+int Tria::EdgeOnSurfaceIndex(void){/*{{{*/
464+
465+ IssmDouble values[NUMVERTICES];
466+ int indices[3][2] = {{1,2},{2,0},{0,1}};
467+
468+ /*Retrieve all inputs and parameters*/
469+ GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
470+
471+ for(int i=0;i<3;i++){
472+ if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
473+ return i;
474+ }
475+ }
476+
477+ _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
478+ _error_("Could not find 2 vertices on surface");
479+}
480+/*}}}*/
481+void Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){/*{{{*/
482+
483+ IssmDouble values[NUMVERTICES];
484+ int indices[3][2] = {{1,2},{2,0},{0,1}};
485+
486+ /*Retrieve all inputs and parameters*/
487+ GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
488+
489+ for(int i=0;i<3;i++){
490+ if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
491+ *pindex1 = indices[i][0];
492+ *pindex2 = indices[i][1];
493+ return;
494+ }
495+ }
496+
497+ _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
498+ _error_("Could not find 2 vertices on surface");
499+}
500+/*}}}*/
501+void Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
502+
503+ switch(response_enum){
504+ case MaterialsRheologyBbarEnum:
505+ *presponse=this->material->GetBbar();
506+ break;
507+
508+ case VelEnum:{
509+
510+ /*Get input:*/
511+ IssmDouble vel;
512+ Input* vel_input;
513+
514+ vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
515+ vel_input->GetInputAverage(&vel);
516+
517+ /*Assign output pointers:*/
518+ *presponse=vel;}
519+ break;
520+ default:
521+ _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
522+ }
523+
524+}
525+/*}}}*/
526 void Tria::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
527
528 IssmDouble xyz_list[NUMVERTICES][3];
529@@ -604,10 +699,90 @@
530 return this->element_type;
531 }
532 /*}}}*/
533-int Tria::ObjectEnum(void){/*{{{*/
534+void Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
535
536- return TriaEnum;
537+ if(!IsOnBase()) return;
538
539+ int approximation;
540+ inputs->GetInputValue(&approximation,ApproximationEnum);
541+
542+ if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
543+ for(int i=0;i<NUMVERTICES;i++){
544+ vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
545+ vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
546+ }
547+ }
548+ else{
549+ /*Intermediaries*/
550+ IssmDouble* xyz_list = NULL;
551+ IssmDouble* xyz_list_base = NULL;
552+ IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base;
553+ IssmDouble bed_normal[2];
554+ IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
555+ IssmDouble surface=0,value=0;
556+ bool grounded;
557+
558+ /* Get node coordinates and dof list: */
559+ GetVerticesCoordinates(&xyz_list);
560+ GetVerticesCoordinatesBase(&xyz_list_base);
561+
562+ /*Retrieve all inputs we will be needing: */
563+ Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
564+ Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
565+ Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
566+ Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
567+ Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
568+
569+ /*Create gauss point in the middle of the basal edge*/
570+ Gauss* gauss=NewGaussBase(1);
571+ gauss->GaussPoint(0);
572+
573+ if(!IsFloating()){
574+ /*Check for basal force only if grounded and touching GL*/
575+ // if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
576+ this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
577+ this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
578+ pressure_input->GetInputValue(&pressure, gauss);
579+ base_input->GetInputValue(&base, gauss);
580+
581+ /*Compute Stress*/
582+ IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
583+ IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
584+ IssmDouble sigma_xy=2.*viscosity*epsilon[2];
585+
586+ /*Get normal vector to the bed */
587+ NormalBase(&bed_normal[0],xyz_list_base);
588+
589+ /*basalforce*/
590+ sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
591+
592+ /*Compute water pressure*/
593+ IssmDouble rho_ice = matpar->GetRhoIce();
594+ IssmDouble rho_water = matpar->GetRhoWater();
595+ IssmDouble gravity = matpar->GetG();
596+ water_pressure=gravity*rho_water*base;
597+
598+ /*Compare basal stress to water pressure and determine whether it should ground*/
599+ if (sigma_nn<water_pressure) grounded=true;
600+ else grounded=false;
601+ }
602+ else{
603+ /*Check for basal elevation if floating*/
604+ base_input->GetInputValue(&base, gauss);
605+ bed_input->GetInputValue(&bed, gauss);
606+ if(base<bed) grounded=true;
607+ else grounded=false;
608+ }
609+ for(int i=0;i<NUMVERTICES;i++){
610+ if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
611+ else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
612+ }
613+
614+ /*clean up*/
615+ delete gauss;
616+ xDelete<IssmDouble>(xyz_list);
617+ xDelete<IssmDouble>(xyz_list_base);
618+ }
619 }
620 /*}}}*/
621 IssmDouble Tria::GetArea(void){/*{{{*/
622@@ -667,56 +842,6 @@
623
624 }
625 /*}}}*/
626-void Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
627-
628- /*Computeportion of the element that has a positive levelset*/
629-
630- bool negative=true;
631- int point;
632- const IssmPDouble epsilon= 1.e-15;
633- IssmDouble f1,f2;
634-
635- /*Be sure that values are not zero*/
636- if(gl[0]==0.) gl[0]=gl[0]+epsilon;
637- if(gl[1]==0.) gl[1]=gl[1]+epsilon;
638- if(gl[2]==0.) gl[2]=gl[2]+epsilon;
639-
640- /*Check that not all nodes are positive or negative*/
641- if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
642- point=0;
643- f1=1.;
644- f2=1.;
645- }
646- else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
647- point=0;
648- f1=0.;
649- f2=0.;
650- }
651- else{
652- if(gl[0]*gl[1]*gl[2]<0) negative=false;
653-
654- if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
655- point=2;
656- f1=gl[2]/(gl[2]-gl[0]);
657- f2=gl[2]/(gl[2]-gl[1]);
658- }
659- else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
660- point=0;
661- f1=gl[0]/(gl[0]-gl[1]);
662- f2=gl[0]/(gl[0]-gl[2]);
663- }
664- else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
665- point=1;
666- f1=gl[1]/(gl[1]-gl[2]);
667- f2=gl[1]/(gl[1]-gl[0]);
668- }
669- }
670- *point1=point;
671- *fraction1=f1;
672- *fraction2=f2;
673- *mainlynegative=negative;
674-}
675-/*}}}*/
676 void Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/
677 /*Computeportion of the element that is grounded*/
678
679@@ -888,149 +1013,6 @@
680 return phi;
681 }
682 /*}}}*/
683-void Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
684-
685- int indices[2];
686- IssmDouble xyz_list[NUMVERTICES][3];
687-
688- /*Element XYZ list*/
689- ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
690-
691- /*Allocate Output*/
692- IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
693- this->EdgeOnBaseIndices(&indices[0],&indices[1]);
694- for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
695-
696- /*Assign output pointer*/
697- *pxyz_list = xyz_list_edge;
698-
699-}/*}}}*/
700-void Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
701-
702- int indices[2];
703- IssmDouble xyz_list[NUMVERTICES][3];
704-
705- /*Element XYZ list*/
706- ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
707-
708- /*Allocate Output*/
709- IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
710- this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
711- for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
712-
713- /*Assign output pointer*/
714- *pxyz_list = xyz_list_edge;
715-
716-}/*}}}*/
717-void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
718-
719- /*Build unit outward pointing vector*/
720- IssmDouble vector[2];
721- IssmDouble norm;
722-
723- vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
724- vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
725-
726- norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
727-
728- normal[0]= + vector[1]/norm;
729- normal[1]= - vector[0]/norm;
730-}
731-/*}}}*/
732-void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
733-
734- int normal_orientation=0;
735- IssmDouble s1,s2;
736- IssmDouble levelset[NUMVERTICES];
737-
738- /*Recover parameters and values*/
739- IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
740- GetInputListOnVertices(&levelset[0],levelsetenum);
741-
742- if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
743- /*Portion of the segments*/
744- s1=levelset[2]/(levelset[2]-levelset[1]);
745- s2=levelset[2]/(levelset[2]-levelset[0]);
746-
747- if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle depending on distribution of levelsetfunction
748- /*New point 1*/
749- xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
750- xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
751- xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
752-
753- /*New point 0*/
754- xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
755- xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
756- xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
757- }
758- else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
759- /*Portion of the segments*/
760- s1=levelset[0]/(levelset[0]-levelset[2]);
761- s2=levelset[0]/(levelset[0]-levelset[1]);
762-
763- if(levelset[0]<0.) normal_orientation=1;
764- /*New point 1*/
765- xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
766- xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
767- xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
768-
769- /*New point 2*/
770- xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
771- xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
772- xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
773- }
774- else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
775- /*Portion of the segments*/
776- s1=levelset[1]/(levelset[1]-levelset[0]);
777- s2=levelset[1]/(levelset[1]-levelset[2]);
778-
779- if(levelset[1]<0.) normal_orientation=1;
780- /*New point 0*/
781- xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
782- xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
783- xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
784-
785- /*New point 2*/
786- xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
787- xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
788- xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
789- }
790- else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
791- xyz_zero[3*0+0]=xyz_list[0*3+0];
792- xyz_zero[3*0+1]=xyz_list[0*3+1];
793- xyz_zero[3*0+2]=xyz_list[0*3+2];
794-
795- /*New point 2*/
796- xyz_zero[3*1+0]=xyz_list[1*3+0];
797- xyz_zero[3*1+1]=xyz_list[1*3+1];
798- xyz_zero[3*1+2]=xyz_list[1*3+2];
799- }
800- else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
801- xyz_zero[3*0+0]=xyz_list[2*3+0];
802- xyz_zero[3*0+1]=xyz_list[2*3+1];
803- xyz_zero[3*0+2]=xyz_list[2*3+2];
804-
805- /*New point 2*/
806- xyz_zero[3*1+0]=xyz_list[0*3+0];
807- xyz_zero[3*1+1]=xyz_list[0*3+1];
808- xyz_zero[3*1+2]=xyz_list[0*3+2];
809- }
810- else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
811- xyz_zero[3*0+0]=xyz_list[1*3+0];
812- xyz_zero[3*0+1]=xyz_list[1*3+1];
813- xyz_zero[3*0+2]=xyz_list[1*3+2];
814-
815- /*New point 2*/
816- xyz_zero[3*1+0]=xyz_list[2*3+0];
817- xyz_zero[3*1+1]=xyz_list[2*3+1];
818- xyz_zero[3*1+2]=xyz_list[2*3+2];
819- }
820- else _error_("Case not covered");
821-
822- /*Assign output pointer*/
823- *pxyz_zero= xyz_zero;
824-}
825-/*}}}*/
826 void Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
827
828 /* Intermediaries */
829@@ -1071,6 +1053,18 @@
830
831 xDelete<int>(indicesfront);
832 }/*}}}*/
833+void Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
834+
835+ Input* input=inputs->GetInput(enumtype);
836+ if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
837+
838+ GaussTria* gauss=new GaussTria();
839+ gauss->GaussVertex(this->GetNodeIndex(node));
840+
841+ input->GetInputValue(pvalue,gauss);
842+ delete gauss;
843+}
844+/*}}}*/
845 void Tria::GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){/*{{{*/
846
847 /* Intermediaries */
848@@ -1111,6 +1105,62 @@
849
850 xDelete<int>(indicesfront);
851 }/*}}}*/
852+void Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
853+
854+ /*Computeportion of the element that has a positive levelset*/
855+
856+ bool negative=true;
857+ int point;
858+ const IssmPDouble epsilon= 1.e-15;
859+ IssmDouble f1,f2;
860+
861+ /*Be sure that values are not zero*/
862+ if(gl[0]==0.) gl[0]=gl[0]+epsilon;
863+ if(gl[1]==0.) gl[1]=gl[1]+epsilon;
864+ if(gl[2]==0.) gl[2]=gl[2]+epsilon;
865+
866+ /*Check that not all nodes are positive or negative*/
867+ if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
868+ point=0;
869+ f1=1.;
870+ f2=1.;
871+ }
872+ else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
873+ point=0;
874+ f1=0.;
875+ f2=0.;
876+ }
877+ else{
878+ if(gl[0]*gl[1]*gl[2]<0) negative=false;
879+
880+ if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
881+ point=2;
882+ f1=gl[2]/(gl[2]-gl[0]);
883+ f2=gl[2]/(gl[2]-gl[1]);
884+ }
885+ else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
886+ point=0;
887+ f1=gl[0]/(gl[0]-gl[1]);
888+ f2=gl[0]/(gl[0]-gl[2]);
889+ }
890+ else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
891+ point=1;
892+ f1=gl[1]/(gl[1]-gl[2]);
893+ f2=gl[1]/(gl[1]-gl[0]);
894+ }
895+ }
896+ *point1=point;
897+ *fraction1=f1;
898+ *fraction2=f2;
899+ *mainlynegative=negative;
900+}
901+/*}}}*/
902+Node* Tria::GetNode(int node_number){/*{{{*/
903+ _assert_(node_number>=0);
904+ _assert_(node_number<this->NumberofNodes(this->element_type));
905+ return this->nodes[node_number];
906+
907+}/*}}}*/
908 int Tria::GetNodeIndex(Node* node){/*{{{*/
909
910 _assert_(nodes);
911@@ -1134,24 +1184,221 @@
912 return NUMVERTICES;
913 }
914 /*}}}*/
915-void Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
916+void Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
917
918- Input* input=inputs->GetInput(enumtype);
919- if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
920+ int *doflist = NULL;
921+ IssmDouble value;
922
923+ /*Fetch number of nodes for this finite element*/
924+ int numnodes = this->NumberofNodes(this->element_type);
925+
926+ /*Fetch dof list and allocate solution vector*/
927+ GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
928+ IssmDouble* values = xNew<IssmDouble>(numnodes);
929+
930+ /*Get inputs*/
931+ Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
932+
933+ /*Ok, we have the values, fill in the array: */
934 GaussTria* gauss=new GaussTria();
935- gauss->GaussVertex(this->GetNodeIndex(node));
936+ for(int i=0;i<numnodes;i++){
937+ gauss->GaussNode(this->element_type,i);
938
939- input->GetInputValue(pvalue,gauss);
940+ enum_input->GetInputValue(&value,gauss);
941+ values[i]=value;
942+ }
943+
944+ solution->SetValues(numnodes,doflist,values,INS_VAL);
945+
946+ /*Free ressources:*/
947+ xDelete<int>(doflist);
948+ xDelete<IssmDouble>(values);
949 delete gauss;
950 }
951 /*}}}*/
952-Node* Tria::GetNode(int node_number){/*{{{*/
953- _assert_(node_number>=0);
954- _assert_(node_number<this->NumberofNodes(this->element_type));
955- return this->nodes[node_number];
956+void Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
957
958+ int vertexidlist[NUMVERTICES];
959+ Input *input=NULL;
960+
961+ /*Get out if this is not an element input*/
962+ if(!IsInput(control_enum)) return;
963+
964+ /*Prepare index list*/
965+ GradientIndexing(&vertexidlist[0],control_index,onsid);
966+
967+ /*Get input (either in element or material)*/
968+ input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
969+
970+ /*Check that it is a ControlInput*/
971+ if (input->ObjectEnum()!=ControlInputEnum){
972+ _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
973+ }
974+
975+ ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
976+}
977+/*}}}*/
978+void Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
979+
980+ int indices[2];
981+ IssmDouble xyz_list[NUMVERTICES][3];
982+
983+ /*Element XYZ list*/
984+ ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
985+
986+ /*Allocate Output*/
987+ IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
988+ this->EdgeOnBaseIndices(&indices[0],&indices[1]);
989+ for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
990+
991+ /*Assign output pointer*/
992+ *pxyz_list = xyz_list_edge;
993+
994 }/*}}}*/
995+void Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
996+
997+ int indices[2];
998+ IssmDouble xyz_list[NUMVERTICES][3];
999+
1000+ /*Element XYZ list*/
1001+ ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
1002+
1003+ /*Allocate Output*/
1004+ IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
1005+ this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
1006+ for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
1007+
1008+ /*Assign output pointer*/
1009+ *pxyz_list = xyz_list_edge;
1010+
1011+}/*}}}*/
1012+bool Tria::HasEdgeOnBase(){/*{{{*/
1013+
1014+ IssmDouble values[NUMVERTICES];
1015+ IssmDouble sum;
1016+
1017+ /*Retrieve all inputs and parameters*/
1018+ GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
1019+ sum = values[0]+values[1]+values[2];
1020+
1021+ _assert_(sum==0. || sum==1. || sum==2.);
1022+
1023+ if(sum==3.) _error_("Two edges on bed not supported yet...");
1024+
1025+ if(sum>1.){
1026+ return true;
1027+ }
1028+ else{
1029+ return false;
1030+ }
1031+}
1032+/*}}}*/
1033+bool Tria::HasEdgeOnSurface(){/*{{{*/
1034+
1035+ IssmDouble values[NUMVERTICES];
1036+ IssmDouble sum;
1037+
1038+ /*Retrieve all inputs and parameters*/
1039+ GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
1040+ sum = values[0]+values[1]+values[2];
1041+
1042+ _assert_(sum==0. || sum==1. || sum==2.);
1043+
1044+ if(sum==3.) _error_("Two edges on surface not supported yet...");
1045+
1046+ if(sum>1.){
1047+ return true;
1048+ }
1049+ else{
1050+ return false;
1051+ }
1052+}
1053+/*}}}*/
1054+IssmDouble Tria::IceVolume(void){/*{{{*/
1055+
1056+ /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
1057+ IssmDouble base,surface,bed;
1058+ IssmDouble xyz_list[NUMVERTICES][3];
1059+
1060+ if(!IsIceInElement())return 0;
1061+
1062+ /*First get back the area of the base*/
1063+ base=this->GetArea();
1064+
1065+ /*Now get the average height*/
1066+ Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
1067+ Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
1068+ surface_input->GetInputAverage(&surface);
1069+ base_input->GetInputAverage(&bed);
1070+
1071+ /*Return: */
1072+ int domaintype;
1073+ parameters->FindParam(&domaintype,DomainTypeEnum);
1074+ if(domaintype==Domain2DverticalEnum){
1075+ return base;
1076+ }
1077+ else{
1078+ return base*(surface-bed);
1079+ }
1080+}
1081+/*}}}*/
1082+IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/
1083+
1084+ /*The volume above floatation: H + rho_water/rho_ice * bathymetry */
1085+ IssmDouble rho_ice,rho_water;
1086+ IssmDouble base,surface,bed,bathymetry;
1087+ IssmDouble xyz_list[NUMVERTICES][3];
1088+
1089+ if(!IsIceInElement() || IsFloating())return 0;
1090+
1091+ rho_ice=matpar->GetRhoIce();
1092+ rho_water=matpar->GetRhoWater();
1093+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
1094+
1095+ /*First calculate the area of the base (cross section triangle)
1096+ * http://en.wikipedia.org/wiki/Triangle
1097+ * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
1098+ 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]));
1099+
1100+ /*Now get the average height and bathymetry*/
1101+ Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
1102+ Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
1103+ Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
1104+ surface_input->GetInputAverage(&surface);
1105+ base_input->GetInputAverage(&bed);
1106+ bed_input->GetInputAverage(&bathymetry);
1107+
1108+ /*Return: */
1109+ return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
1110+}
1111+/*}}}*/
1112+void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
1113+
1114+ /*Intermediary*/
1115+ int num_controls;
1116+ int* control_type=NULL;
1117+ Input* input=NULL;
1118+
1119+ /*retrieve some parameters: */
1120+ this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
1121+ this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
1122+
1123+ for(int i=0;i<num_controls;i++){
1124+ input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
1125+ if (input->ObjectEnum()!=ControlInputEnum){
1126+ _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
1127+ }
1128+
1129+ ((ControlInput*)input)->UpdateValue(scalar);
1130+ ((ControlInput*)input)->Constrain();
1131+ if (save_parameter) ((ControlInput*)input)->SaveValue();
1132+
1133+ }
1134+
1135+ /*Clean up and return*/
1136+ xDelete<int>(control_type);
1137+}
1138+/*}}}*/
1139 void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
1140
1141 /*New input*/
1142@@ -1370,6 +1617,59 @@
1143
1144 }
1145 /*}}}*/
1146+bool Tria::IsFaceOnBoundary(void){/*{{{*/
1147+
1148+ IssmDouble values[NUMVERTICES];
1149+ IssmDouble sum;
1150+
1151+ /*Retrieve all inputs and parameters*/
1152+ GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
1153+ sum = values[0]+values[1]+values[2];
1154+
1155+ _assert_(sum==0. || sum==1. || sum==2.);
1156+
1157+ if(sum==3.) _error_("Two edges on boundary not supported yet...");
1158+
1159+ if(sum>1.){
1160+ return true;
1161+ }
1162+ else{
1163+ return false;
1164+ }
1165+}/*}}}*/
1166+bool Tria::IsIcefront(void){/*{{{*/
1167+
1168+ bool isicefront;
1169+ int i,nrice;
1170+ IssmDouble ls[NUMVERTICES];
1171+
1172+ /*Retrieve all inputs and parameters*/
1173+ GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
1174+
1175+ /* If only one vertex has ice, there is an ice front here */
1176+ isicefront=false;
1177+ if(IsIceInElement()){
1178+ nrice=0;
1179+ for(i=0;i<NUMVERTICES;i++)
1180+ if(ls[i]<0.) nrice++;
1181+ if(nrice==1) isicefront= true;
1182+ }
1183+ return isicefront;
1184+}/*}}}*/
1185+bool Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
1186+
1187+ int i;
1188+ bool shelf=false;
1189+
1190+ for(i=0;i<NUMVERTICES;i++){
1191+ if (flags[vertices[i]->Pid()]<0.){
1192+ shelf=true;
1193+ break;
1194+ }
1195+ }
1196+ return shelf;
1197+}
1198+/*}}}*/
1199 bool Tria::IsOnBase(){/*{{{*/
1200
1201 int domaintype;
1202@@ -1396,6 +1696,25 @@
1203 }
1204 }
1205 /*}}}*/
1206+bool Tria::IsZeroLevelset(int levelset_enum){/*{{{*/
1207+
1208+ bool iszerols;
1209+ IssmDouble ls[NUMVERTICES];
1210+
1211+ /*Retrieve all inputs and parameters*/
1212+ GetInputListOnVertices(&ls[0],levelset_enum);
1213+
1214+ /*If the level set is awlays <0, there is no ice front here*/
1215+ iszerols= false;
1216+ if(IsIceInElement()){
1217+ if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
1218+ iszerols = true;
1219+ }
1220+ }
1221+
1222+ return iszerols;
1223+}
1224+/*}}}*/
1225 void Tria::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1226
1227 _assert_(gauss->Enum()==GaussTriaEnum);
1228@@ -1424,222 +1743,267 @@
1229
1230 }
1231 /*}}}*/
1232-bool Tria::HasEdgeOnBase(){/*{{{*/
1233+IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
1234
1235- IssmDouble values[NUMVERTICES];
1236- IssmDouble sum;
1237
1238- /*Retrieve all inputs and parameters*/
1239- GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
1240- sum = values[0]+values[1]+values[2];
1241+ /*intermediary: */
1242+ IssmDouble* values=NULL;
1243+ Input* thickness_input=NULL;
1244+ IssmDouble thickness;
1245+ IssmDouble weight;
1246+ IssmDouble Jdet;
1247+ IssmDouble volume;
1248+ IssmDouble rho_ice;
1249+ IssmDouble* xyz_list=NULL;
1250+ int point1;
1251+ IssmDouble fraction1,fraction2;
1252+ bool mainlynegative=true;
1253+
1254+ /*Output:*/
1255+ volume=0;
1256
1257- _assert_(sum==0. || sum==1. || sum==2.);
1258+ /* Get node coordinates and dof list: */
1259+ GetVerticesCoordinates(&xyz_list);
1260
1261- if(sum==3.) _error_("Two edges on bed not supported yet...");
1262+ /*Retrieve inputs required:*/
1263+ thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
1264+
1265+ /*Retrieve material parameters: */
1266+ rho_ice=matpar->GetRhoIce();
1267
1268- if(sum>1.){
1269- return true;
1270+ /*Retrieve values of the levelset defining the masscon: */
1271+ values = xNew<IssmDouble>(NUMVERTICES);
1272+ for(int i=0;i<NUMVERTICES;i++){
1273+ values[i]=levelset[this->vertices[i]->Sid()];
1274 }
1275- else{
1276- return false;
1277+
1278+ /*Ok, use the level set values to figure out where we put our gaussian points:*/
1279+ this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
1280+ Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
1281+
1282+ volume=0;
1283+
1284+ for(int ig=gauss->begin();ig<gauss->end();ig++){
1285+ gauss->GaussPoint(ig);
1286+
1287+ this->JacobianDeterminant(&Jdet,xyz_list,gauss);
1288+ thickness_input->GetInputValue(&thickness, gauss);
1289+
1290+ volume+=thickness*gauss->weight*Jdet;
1291 }
1292+
1293+ /* clean up and Return: */
1294+ xDelete<IssmDouble>(xyz_list);
1295+ xDelete<IssmDouble>(values);
1296+ delete gauss;
1297+ return rho_ice*volume;
1298 }
1299 /*}}}*/
1300-bool Tria::HasEdgeOnSurface(){/*{{{*/
1301+IssmDouble Tria::MassFlux(IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
1302
1303- IssmDouble values[NUMVERTICES];
1304- IssmDouble sum;
1305+ int domaintype;
1306+ IssmDouble mass_flux=0.;
1307+ IssmDouble xyz_list[NUMVERTICES][3];
1308+ IssmDouble normal[2];
1309+ IssmDouble length,rho_ice;
1310+ IssmDouble h1,h2;
1311+ IssmDouble vx1,vx2,vy1,vy2;
1312+ GaussTria* gauss_1=NULL;
1313+ GaussTria* gauss_2=NULL;
1314
1315- /*Retrieve all inputs and parameters*/
1316- GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
1317- sum = values[0]+values[1]+values[2];
1318+ /*Get material parameters :*/
1319+ rho_ice=matpar->GetRhoIce();
1320
1321- _assert_(sum==0. || sum==1. || sum==2.);
1322+ /*First off, check that this segment belongs to this element: */
1323+ if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id);
1324
1325- if(sum==3.) _error_("Two edges on surface not supported yet...");
1326+ /*Get xyz list: */
1327+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
1328
1329- if(sum>1.){
1330- return true;
1331+ /*get area coordinates of 0 and 1 locations: */
1332+ gauss_1=new GaussTria();
1333+ gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
1334+ gauss_2=new GaussTria();
1335+ gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
1336+
1337+ normal[0]=cos(atan2(x1-x2,y2-y1));
1338+ normal[1]=sin(atan2(x1-x2,y2-y1));
1339+
1340+ length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
1341+
1342+ Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
1343+ this->parameters->FindParam(&domaintype,DomainTypeEnum);
1344+ Input* vx_input=NULL;
1345+ Input* vy_input=NULL;
1346+ if(domaintype==Domain2DhorizontalEnum){
1347+ vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
1348+ vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
1349 }
1350 else{
1351- return false;
1352+ vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
1353+ vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
1354 }
1355-}
1356-/*}}}*/
1357-void Tria::EdgeOnBaseIndices(int* pindex1,int* pindex2){/*{{{*/
1358
1359- IssmDouble values[NUMVERTICES];
1360- int indices[3][2] = {{1,2},{2,0},{0,1}};
1361+ thickness_input->GetInputValue(&h1, gauss_1);
1362+ thickness_input->GetInputValue(&h2, gauss_2);
1363+ vx_input->GetInputValue(&vx1,gauss_1);
1364+ vx_input->GetInputValue(&vx2,gauss_2);
1365+ vy_input->GetInputValue(&vy1,gauss_1);
1366+ vy_input->GetInputValue(&vy2,gauss_2);
1367
1368- /*Retrieve all inputs and parameters*/
1369- GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
1370+ mass_flux= rho_ice*length*(
1371+ (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
1372+ (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
1373+ );
1374
1375- for(int i=0;i<3;i++){
1376- if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1377- *pindex1 = indices[i][0];
1378- *pindex2 = indices[i][1];
1379- return;
1380- }
1381- }
1382-
1383- _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
1384- _error_("Could not find 2 vertices on bed");
1385+ /*clean up and return:*/
1386+ delete gauss_1;
1387+ delete gauss_2;
1388+ return mass_flux;
1389 }
1390 /*}}}*/
1391-void Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){/*{{{*/
1392+IssmDouble Tria::MassFlux(IssmDouble* segment){/*{{{*/
1393
1394- IssmDouble values[NUMVERTICES];
1395- int indices[3][2] = {{1,2},{2,0},{0,1}};
1396+ int domaintype;
1397+ IssmDouble mass_flux=0.;
1398+ IssmDouble xyz_list[NUMVERTICES][3];
1399+ IssmDouble normal[2];
1400+ IssmDouble length,rho_ice;
1401+ IssmDouble x1,y1,x2,y2,h1,h2;
1402+ IssmDouble vx1,vx2,vy1,vy2;
1403+ GaussTria* gauss_1=NULL;
1404+ GaussTria* gauss_2=NULL;
1405
1406- /*Retrieve all inputs and parameters*/
1407- GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
1408+ /*Get material parameters :*/
1409+ rho_ice=matpar->GetRhoIce();
1410
1411- for(int i=0;i<3;i++){
1412- if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1413- *pindex1 = indices[i][0];
1414- *pindex2 = indices[i][1];
1415- return;
1416- }
1417- }
1418+ /*First off, check that this segment belongs to this element: */
1419+ if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
1420
1421- _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
1422- _error_("Could not find 2 vertices on surface");
1423-}
1424-/*}}}*/
1425-int Tria::EdgeOnBaseIndex(void){/*{{{*/
1426+ /*Recover segment node locations: */
1427+ x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
1428
1429- IssmDouble values[NUMVERTICES];
1430- int indices[3][2] = {{1,2},{2,0},{0,1}};
1431+ /*Get xyz list: */
1432+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
1433
1434- /*Retrieve all inputs and parameters*/
1435- GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
1436+ /*get area coordinates of 0 and 1 locations: */
1437+ gauss_1=new GaussTria();
1438+ gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
1439+ gauss_2=new GaussTria();
1440+ gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
1441
1442- for(int i=0;i<3;i++){
1443- if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1444- return i;
1445- }
1446+ normal[0]=cos(atan2(x1-x2,y2-y1));
1447+ normal[1]=sin(atan2(x1-x2,y2-y1));
1448+
1449+ length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
1450+
1451+ Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
1452+ this->parameters->FindParam(&domaintype,DomainTypeEnum);
1453+ Input* vx_input=NULL;
1454+ Input* vy_input=NULL;
1455+ if(domaintype==Domain2DhorizontalEnum){
1456+ vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
1457+ vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
1458 }
1459+ else{
1460+ vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
1461+ vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
1462+ }
1463
1464- _printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
1465- _error_("Could not find 2 vertices on bed");
1466+ thickness_input->GetInputValue(&h1, gauss_1);
1467+ thickness_input->GetInputValue(&h2, gauss_2);
1468+ vx_input->GetInputValue(&vx1,gauss_1);
1469+ vx_input->GetInputValue(&vx2,gauss_2);
1470+ vy_input->GetInputValue(&vy1,gauss_1);
1471+ vy_input->GetInputValue(&vy2,gauss_2);
1472+
1473+ mass_flux= rho_ice*length*(
1474+ (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
1475+ (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
1476+ );
1477+
1478+ /*clean up and return:*/
1479+ delete gauss_1;
1480+ delete gauss_2;
1481+ return mass_flux;
1482 }
1483 /*}}}*/
1484-int Tria::EdgeOnSurfaceIndex(void){/*{{{*/
1485+IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/
1486
1487- IssmDouble values[NUMVERTICES];
1488- int indices[3][2] = {{1,2},{2,0},{0,1}};
1489+ /*Intermediaries*/
1490+ IssmDouble model,observation,weight;
1491+ IssmDouble Jdet;
1492+ IssmDouble Jelem = 0;
1493+ IssmDouble xyz_list[NUMVERTICES][3];
1494+ GaussTria *gauss = NULL;
1495
1496- /*Retrieve all inputs and parameters*/
1497- GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
1498+ /*If on water, return 0: */
1499+ if(!IsIceInElement())return 0;
1500
1501- for(int i=0;i<3;i++){
1502- if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
1503- return i;
1504- }
1505- }
1506+ /*Retrieve all inputs we will be needing: */
1507+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
1508+ Input* model_input=inputs->GetInput(modelenum); _assert_(model_input);
1509+ Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input);
1510+ Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input);
1511
1512- _printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
1513- _error_("Could not find 2 vertices on surface");
1514-}
1515-/*}}}*/
1516-void Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
1517+ /* Start looping on the number of gaussian points: */
1518+ gauss=new GaussTria(2);
1519+ for(int ig=gauss->begin();ig<gauss->end();ig++){
1520
1521- if(!IsOnBase()) return;
1522+ gauss->GaussPoint(ig);
1523
1524- int approximation;
1525- inputs->GetInputValue(&approximation,ApproximationEnum);
1526+ /* Get Jacobian determinant: */
1527+ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
1528
1529- if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
1530- for(int i=0;i<NUMVERTICES;i++){
1531- vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
1532- vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
1533- }
1534+ /*Get parameters at gauss point*/
1535+ model_input->GetInputValue(&model,gauss);
1536+ observation_input->GetInputValue(&observation,gauss);
1537+ weights_input->GetInputValue(&weight,gauss);
1538+
1539+ /*compute misfit between model and observation */
1540+ Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight;
1541 }
1542- else{
1543- /*Intermediaries*/
1544- IssmDouble* xyz_list = NULL;
1545- IssmDouble* xyz_list_base = NULL;
1546- IssmDouble pressure,water_pressure,sigma_nn,viscosity,bed,base;
1547- IssmDouble bed_normal[2];
1548- IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
1549- IssmDouble surface=0,value=0;
1550- bool grounded;
1551
1552- /* Get node coordinates and dof list: */
1553- GetVerticesCoordinates(&xyz_list);
1554- GetVerticesCoordinatesBase(&xyz_list_base);
1555+ /* clean up and Return: */
1556+ delete gauss;
1557+ return Jelem;
1558+}
1559+/*}}}*/
1560+IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/
1561
1562- /*Retrieve all inputs we will be needing: */
1563- Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
1564- Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
1565- Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
1566- Input* vx_input = inputs->GetInput(VxEnum); _assert_(vx_input);
1567- Input* vy_input = inputs->GetInput(VyEnum); _assert_(vy_input);
1568+ /*Intermediaries*/
1569+ IssmDouble weight;
1570+ IssmDouble Jdet;
1571+ IssmDouble Jelem = 0;
1572+ IssmDouble xyz_list[NUMVERTICES][3];
1573+ GaussTria *gauss = NULL;
1574
1575- /*Create gauss point in the middle of the basal edge*/
1576- Gauss* gauss=NewGaussBase(1);
1577- gauss->GaussPoint(0);
1578+ /*If on water, return 0: */
1579+ if(!IsIceInElement())return 0;
1580
1581- if(!IsFloating()){
1582- /*Check for basal force only if grounded and touching GL*/
1583- // if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
1584- this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
1585- this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
1586- pressure_input->GetInputValue(&pressure, gauss);
1587- base_input->GetInputValue(&base, gauss);
1588+ /*Retrieve all inputs we will be needing: */
1589+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
1590+ Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input);
1591
1592- /*Compute Stress*/
1593- IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
1594- IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
1595- IssmDouble sigma_xy=2.*viscosity*epsilon[2];
1596+ /* Start looping on the number of gaussian points: */
1597+ gauss=new GaussTria(2);
1598+ for(int ig=gauss->begin();ig<gauss->end();ig++){
1599
1600- /*Get normal vector to the bed */
1601- NormalBase(&bed_normal[0],xyz_list_base);
1602+ gauss->GaussPoint(ig);
1603
1604- /*basalforce*/
1605- sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
1606+ /* Get Jacobian determinant: */
1607+ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
1608
1609- /*Compute water pressure*/
1610- IssmDouble rho_ice = matpar->GetRhoIce();
1611- IssmDouble rho_water = matpar->GetRhoWater();
1612- IssmDouble gravity = matpar->GetG();
1613- water_pressure=gravity*rho_water*base;
1614+ /*Get parameters at gauss point*/
1615+ weights_input->GetInputValue(&weight,gauss);
1616
1617- /*Compare basal stress to water pressure and determine whether it should ground*/
1618- if (sigma_nn<water_pressure) grounded=true;
1619- else grounded=false;
1620- }
1621- else{
1622- /*Check for basal elevation if floating*/
1623- base_input->GetInputValue(&base, gauss);
1624- bed_input->GetInputValue(&bed, gauss);
1625- if(base<bed) grounded=true;
1626- else grounded=false;
1627- }
1628- for(int i=0;i<NUMVERTICES;i++){
1629- if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
1630- else vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
1631- }
1632-
1633- /*clean up*/
1634- delete gauss;
1635- xDelete<IssmDouble>(xyz_list);
1636- xDelete<IssmDouble>(xyz_list_base);
1637+ /*compute misfit between model and observation */
1638+ Jelem+=Jdet*weight*gauss->weight;
1639 }
1640-}
1641-/*}}}*/
1642-bool Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
1643
1644- int i;
1645- bool shelf=false;
1646-
1647- for(i=0;i<NUMVERTICES;i++){
1648- if (flags[vertices[i]->Pid()]<0.){
1649- shelf=true;
1650- break;
1651- }
1652- }
1653- return shelf;
1654+ /* clean up and Return: */
1655+ delete gauss;
1656+ return Jelem;
1657 }
1658 /*}}}*/
1659 Gauss* Tria::NewGauss(void){/*{{{*/
1660@@ -1690,59 +2054,59 @@
1661
1662 }
1663 /*}}}*/
1664-void Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
1665+void Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1666
1667 _assert_(gauss->Enum()==GaussTriaEnum);
1668- this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
1669+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
1670
1671 }
1672 /*}}}*/
1673-void Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
1674+void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1675
1676 _assert_(gauss->Enum()==GaussTriaEnum);
1677- this->GetNodalFunctions(basis,(GaussTria*)gauss,P2Enum);
1678+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());
1679
1680 }
1681 /*}}}*/
1682-void Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1683+void Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
1684
1685 _assert_(gauss->Enum()==GaussTriaEnum);
1686- this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
1687+ this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());
1688
1689 }
1690 /*}}}*/
1691-void Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1692+void Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
1693
1694 _assert_(gauss->Enum()==GaussTriaEnum);
1695- this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
1696+ this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
1697
1698 }
1699 /*}}}*/
1700-void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1701+void Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
1702
1703 _assert_(gauss->Enum()==GaussTriaEnum);
1704- this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());
1705+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
1706
1707 }
1708 /*}}}*/
1709-void Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
1710+void Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
1711
1712 _assert_(gauss->Enum()==GaussTriaEnum);
1713- this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());
1714+ this->GetNodalFunctions(basis,(GaussTria*)gauss,P2Enum);
1715
1716 }
1717 /*}}}*/
1718-void Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
1719+void Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
1720
1721 _assert_(gauss->Enum()==GaussTriaEnum);
1722- this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());
1723+ this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());
1724
1725 }
1726 /*}}}*/
1727-void Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
1728+void Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
1729
1730 _assert_(gauss->Enum()==GaussTriaEnum);
1731- this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());
1732+ this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());
1733
1734 }
1735 /*}}}*/
1736@@ -1794,6 +2158,21 @@
1737 _assert_(bed_normal[1]<0);
1738 }
1739 /*}}}*/
1740+void Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
1741+
1742+ /*Build unit outward pointing vector*/
1743+ IssmDouble vector[2];
1744+ IssmDouble norm;
1745+
1746+ vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
1747+ vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
1748+
1749+ norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
1750+
1751+ normal[0]= + vector[1]/norm;
1752+ normal[1]= - vector[0]/norm;
1753+}
1754+/*}}}*/
1755 void Tria::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){/*{{{*/
1756
1757 /*Build unit outward pointing vector*/
1758@@ -1812,18 +2191,16 @@
1759 _assert_(top_normal[1]>0);
1760 }
1761 /*}}}*/
1762-int Tria::VelocityInterpolation(void){/*{{{*/
1763- return TriaRef::VelocityInterpolation(this->element_type);
1764+int Tria::ObjectEnum(void){/*{{{*/
1765+
1766+ return TriaEnum;
1767+
1768 }
1769 /*}}}*/
1770 int Tria::PressureInterpolation(void){/*{{{*/
1771 return TriaRef::PressureInterpolation(this->element_type);
1772 }
1773 /*}}}*/
1774-int Tria::TensorInterpolation(void){/*{{{*/
1775- return TriaRef::TensorInterpolation(this->element_type);
1776-}
1777-/*}}}*/
1778 int Tria::NumberofNodesPressure(void){/*{{{*/
1779 return TriaRef::NumberofNodes(this->PressureInterpolation());
1780 }
1781@@ -1900,6 +2277,33 @@
1782 delete gauss;
1783 }
1784 /*}}}*/
1785+void Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
1786+
1787+ IssmDouble h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
1788+ IssmDouble bed_hydro;
1789+ IssmDouble rho_water,rho_ice,density;
1790+
1791+ /*material parameters: */
1792+ rho_water=matpar->GetRhoWater();
1793+ rho_ice=matpar->GetRhoIce();
1794+ density=rho_ice/rho_water;
1795+ GetInputListOnVertices(&h[0],ThicknessEnum);
1796+ GetInputListOnVertices(&r[0],BedEnum);
1797+ GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
1798+
1799+ /*go through vertices, and figure out which ones are grounded and want to unground: */
1800+ for(int i=0;i<NUMVERTICES;i++){
1801+ /*Find if grounded vertices want to start floating*/
1802+ if (gl[i]>0.){
1803+ bed_hydro=-density*h[i];
1804+ if(bed_hydro>r[i]){
1805+ /*Vertex that could potentially unground, flag it*/
1806+ potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
1807+ }
1808+ }
1809+ }
1810+}
1811+/*}}}*/
1812 void Tria::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/
1813
1814 /*Static condensation if requested*/
1815@@ -2010,6 +2414,82 @@
1816 _error_("not implemented yet");
1817 }
1818 /*}}}*/
1819+void Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
1820+
1821+ IssmDouble values[NUMVERTICES];
1822+ int vertexpidlist[NUMVERTICES],control_init;
1823+
1824+
1825+ /*Get Domain type*/
1826+ int domaintype;
1827+ parameters->FindParam(&domaintype,DomainTypeEnum);
1828+
1829+ /*Specific case for depth averaged quantities*/
1830+ control_init=control_enum;
1831+ if(domaintype==Domain2DverticalEnum){
1832+ if(control_enum==MaterialsRheologyBbarEnum){
1833+ control_enum=MaterialsRheologyBEnum;
1834+ if(!IsOnBase()) return;
1835+ }
1836+ if(control_enum==DamageDbarEnum){
1837+ control_enum=DamageDEnum;
1838+ if(!IsOnBase()) return;
1839+ }
1840+ }
1841+
1842+ /*Get out if this is not an element input*/
1843+ if(!IsInput(control_enum)) return;
1844+
1845+ /*Prepare index list*/
1846+ GradientIndexing(&vertexpidlist[0],control_index);
1847+
1848+ /*Get values on vertices*/
1849+ for(int i=0;i<NUMVERTICES;i++){
1850+ values[i]=vector[vertexpidlist[i]];
1851+ }
1852+ Input* new_input = new TriaInput(control_enum,values,P1Enum);
1853+ Input* input = (Input*)this->inputs->GetInput(control_enum); _assert_(input);
1854+ if(input->ObjectEnum()!=ControlInputEnum){
1855+ _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
1856+ }
1857+
1858+ ((ControlInput*)input)->SetInput(new_input);
1859+}
1860+/*}}}*/
1861+void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
1862+
1863+ /*go into parameters and get the analysis_counter: */
1864+ int analysis_counter;
1865+ parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
1866+
1867+ /*Get Element type*/
1868+ if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter];
1869+
1870+ /*Pick up nodes*/
1871+ if(this->hnodes && this->hnodes[analysis_counter]){
1872+ this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
1873+ }
1874+
1875+}
1876+/*}}}*/
1877+Element* Tria::SpawnBasalElement(void){/*{{{*/
1878+
1879+ int index1,index2;
1880+ int domaintype;
1881+
1882+ this->parameters->FindParam(&domaintype,DomainTypeEnum);
1883+ switch(domaintype){
1884+ case Domain2DhorizontalEnum:
1885+ return this;
1886+ case Domain2DverticalEnum:
1887+ _assert_(HasEdgeOnBase());
1888+ this->EdgeOnBaseIndices(&index1,&index2);
1889+ return SpawnSeg(index1,index2);
1890+ default:
1891+ _error_("not implemented yet");
1892+ }
1893+}
1894+/*}}}*/
1895 Seg* Tria::SpawnSeg(int index1,int index2){/*{{{*/
1896
1897 int analysis_counter;
1898@@ -2037,7 +2517,7 @@
1899 return seg;
1900 }
1901 /*}}}*/
1902-Element* Tria::SpawnBasalElement(void){/*{{{*/
1903+Element* Tria::SpawnTopElement(void){/*{{{*/
1904
1905 int index1,index2;
1906 int domaintype;
1907@@ -2047,46 +2527,104 @@
1908 case Domain2DhorizontalEnum:
1909 return this;
1910 case Domain2DverticalEnum:
1911- _assert_(HasEdgeOnBase());
1912- this->EdgeOnBaseIndices(&index1,&index2);
1913- return SpawnSeg(index1,index2);
1914+ _assert_(HasEdgeOnSurface());
1915+ this->EdgeOnSurfaceIndices(&index1,&index2);
1916+ return SpawnSeg(index2,index1); //reverse order
1917 default:
1918 _error_("not implemented yet");
1919 }
1920 }
1921 /*}}}*/
1922-Element* Tria::SpawnTopElement(void){/*{{{*/
1923+void Tria::StrainRateparallel(){/*{{{*/
1924
1925- int index1,index2;
1926- int domaintype;
1927+ IssmDouble *xyz_list = NULL;
1928+ IssmDouble epsilon[3];
1929+ GaussTria* gauss=NULL;
1930+ IssmDouble vx,vy,vel;
1931+ IssmDouble strainxx;
1932+ IssmDouble strainxy;
1933+ IssmDouble strainyy;
1934+ IssmDouble strainparallel[NUMVERTICES];
1935
1936- this->parameters->FindParam(&domaintype,DomainTypeEnum);
1937- switch(domaintype){
1938- case Domain2DhorizontalEnum:
1939- return this;
1940- case Domain2DverticalEnum:
1941- _assert_(HasEdgeOnSurface());
1942- this->EdgeOnSurfaceIndices(&index1,&index2);
1943- return SpawnSeg(index2,index1); //reverse order
1944- default:
1945- _error_("not implemented yet");
1946+ /* Get node coordinates and dof list: */
1947+ this->GetVerticesCoordinates(&xyz_list);
1948+
1949+ /*Retrieve all inputs we will need*/
1950+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
1951+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
1952+
1953+ /* Start looping on the number of vertices: */
1954+ gauss=new GaussTria();
1955+ for (int iv=0;iv<NUMVERTICES;iv++){
1956+ gauss->GaussVertex(iv);
1957+
1958+ /* Get the value we need*/
1959+ vx_input->GetInputValue(&vx,gauss);
1960+ vy_input->GetInputValue(&vy,gauss);
1961+ vel=vx*vx+vy*vy;
1962+
1963+ /*Compute strain rate viscosity and pressure: */
1964+ this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
1965+ strainxx=epsilon[0];
1966+ strainyy=epsilon[1];
1967+ strainxy=epsilon[2];
1968+
1969+ /*strainparallel= Strain rate along the ice flow direction */
1970+ strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
1971 }
1972+
1973+ /*Add input*/
1974+ this->inputs->AddInput(new TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
1975+
1976+ /*Clean up and return*/
1977+ delete gauss;
1978+ xDelete<IssmDouble>(xyz_list);
1979 }
1980 /*}}}*/
1981-void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
1982+void Tria::StrainRateperpendicular(){/*{{{*/
1983
1984- /*go into parameters and get the analysis_counter: */
1985- int analysis_counter;
1986- parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
1987+ IssmDouble *xyz_list = NULL;
1988+ GaussTria* gauss=NULL;
1989+ IssmDouble epsilon[3];
1990+ IssmDouble vx,vy,vel;
1991+ IssmDouble strainxx;
1992+ IssmDouble strainxy;
1993+ IssmDouble strainyy;
1994+ IssmDouble strainperpendicular[NUMVERTICES];
1995
1996- /*Get Element type*/
1997- if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter];
1998+ /* Get node coordinates and dof list: */
1999+ this->GetVerticesCoordinates(&xyz_list);
2000
2001- /*Pick up nodes*/
2002- if(this->hnodes && this->hnodes[analysis_counter]){
2003- this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
2004+ /*Retrieve all inputs we will need*/
2005+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2006+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2007+
2008+ /* Start looping on the number of vertices: */
2009+ gauss=new GaussTria();
2010+ for (int iv=0;iv<NUMVERTICES;iv++){
2011+ gauss->GaussVertex(iv);
2012+
2013+ /* Get the value we need*/
2014+ vx_input->GetInputValue(&vx,gauss);
2015+ vy_input->GetInputValue(&vy,gauss);
2016+ vel=vx*vx+vy*vy;
2017+
2018+ /*Compute strain rate viscosity and pressure: */
2019+ this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
2020+ strainxx=epsilon[0];
2021+ strainyy=epsilon[1];
2022+ strainxy=epsilon[2];
2023+
2024+ /*strainperpendicular= Strain rate perpendicular to the ice flow direction */
2025+ strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
2026 }
2027
2028+ /*Add input*/
2029+ this->inputs->AddInput(new TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
2030+
2031+ /*Clean up and return*/
2032+ delete gauss;
2033+ xDelete<IssmDouble>(xyz_list);
2034 }
2035 /*}}}*/
2036 IssmDouble Tria::SurfaceArea(void){/*{{{*/
2037@@ -2116,6 +2654,10 @@
2038 return S;
2039 }
2040 /*}}}*/
2041+int Tria::TensorInterpolation(void){/*{{{*/
2042+ return TriaRef::TensorInterpolation(this->element_type);
2043+}
2044+/*}}}*/
2045 IssmDouble Tria::TimeAdapt(void){/*{{{*/
2046
2047 /*intermediary: */
2048@@ -2157,6 +2699,34 @@
2049 return dt;
2050 }
2051 /*}}}*/
2052+IssmDouble Tria::TotalSmb(void){/*{{{*/
2053+
2054+ /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
2055+ IssmDouble base,smb,rho_ice;
2056+ IssmDouble Total_Smb=0;
2057+ IssmDouble xyz_list[NUMVERTICES][3];
2058+
2059+ /*Get material parameters :*/
2060+ rho_ice=matpar->GetRhoIce();
2061+
2062+ if(!IsIceInElement())return 0;
2063+
2064+ ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
2065+
2066+ /*First calculate the area of the base (cross section triangle)
2067+ * http://en.wikipedia.org/wiki/Triangle
2068+ * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
2069+ base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
2070+
2071+ /*Now get the average SMB over the element*/
2072+ Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
2073+ smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1
2074+ Total_Smb=rho_ice*base*smb; // smb on element in kg s-1
2075+
2076+ /*Return: */
2077+ return Total_Smb;
2078+}
2079+/*}}}*/
2080 void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){/*{{{*/
2081
2082 /*Intermediaries*/
2083@@ -2346,487 +2916,134 @@
2084
2085 }
2086 /*}}}*/
2087-void Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
2088- TriaRef::GetInputValue(pvalue,values,gauss,P1Enum);
2089-}
2090-/*}}}*/
2091-void Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
2092- TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
2093-}
2094-/*}}}*/
2095-int Tria::VertexConnectivity(int vertexindex){/*{{{*/
2096- _assert_(this->vertices);
2097- return this->vertices[vertexindex]->Connectivity();
2098-}
2099-/*}}}*/
2100-bool Tria::IsZeroLevelset(int levelset_enum){/*{{{*/
2101+int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
2102
2103- bool iszerols;
2104- IssmDouble ls[NUMVERTICES];
2105+ int i;
2106+ int nflipped=0;
2107
2108- /*Retrieve all inputs and parameters*/
2109- GetInputListOnVertices(&ls[0],levelset_enum);
2110+ /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
2111+ for(i=0;i<3;i++){
2112+ if (reCast<bool>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
2113+ vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
2114
2115- /*If the level set is awlays <0, there is no ice front here*/
2116- iszerols= false;
2117- if(IsIceInElement()){
2118- if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
2119- iszerols = true;
2120+ /*If node was not on ice shelf, we flipped*/
2121+ if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
2122+ nflipped++;
2123+ }
2124 }
2125 }
2126-
2127- return iszerols;
2128+ return nflipped;
2129 }
2130 /*}}}*/
2131-bool Tria::IsIcefront(void){/*{{{*/
2132-
2133- bool isicefront;
2134- int i,nrice;
2135- IssmDouble ls[NUMVERTICES];
2136-
2137- /*Retrieve all inputs and parameters*/
2138- GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
2139-
2140- /* If only one vertex has ice, there is an ice front here */
2141- isicefront=false;
2142- if(IsIceInElement()){
2143- nrice=0;
2144- for(i=0;i<NUMVERTICES;i++)
2145- if(ls[i]<0.) nrice++;
2146- if(nrice==1) isicefront= true;
2147- }
2148- return isicefront;
2149-}/*}}}*/
2150-bool Tria::IsFaceOnBoundary(void){/*{{{*/
2151-
2152- IssmDouble values[NUMVERTICES];
2153- IssmDouble sum;
2154-
2155- /*Retrieve all inputs and parameters*/
2156- GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
2157- sum = values[0]+values[1]+values[2];
2158-
2159- _assert_(sum==0. || sum==1. || sum==2.);
2160-
2161- if(sum==3.) _error_("Two edges on boundary not supported yet...");
2162-
2163- if(sum>1.){
2164- return true;
2165- }
2166- else{
2167- return false;
2168- }
2169-}/*}}}*/
2170-void Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
2171-
2172- bool already = false;
2173- int i,j;
2174- int partition[NUMVERTICES];
2175- int offsetsid[NUMVERTICES];
2176- int offsetdof[NUMVERTICES];
2177- IssmDouble area;
2178- IssmDouble mean;
2179-
2180- /*First, get the area: */
2181- area=this->GetArea();
2182-
2183- /*Figure out the average for this element: */
2184- this->GetVerticesSidList(&offsetsid[0]);
2185- this->GetVertexPidList(&offsetdof[0]);
2186- mean=0;
2187- for(i=0;i<NUMVERTICES;i++){
2188- partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
2189- mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
2190- }
2191-
2192- /*Add contribution: */
2193- for(i=0;i<NUMVERTICES;i++){
2194- already=false;
2195- for(j=0;j<i;j++){
2196- if (partition[i]==partition[j]){
2197- already=true;
2198- break;
2199- }
2200- }
2201- if(!already){
2202- partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
2203- partition_areas->SetValue(partition[i],area,ADD_VAL);
2204- };
2205- }
2206+void Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
2207+ TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
2208 }
2209 /*}}}*/
2210-IssmDouble Tria::IceVolume(void){/*{{{*/
2211-
2212- /*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
2213- IssmDouble base,surface,bed;
2214- IssmDouble xyz_list[NUMVERTICES][3];
2215-
2216- if(!IsIceInElement())return 0;
2217-
2218- /*First get back the area of the base*/
2219- base=this->GetArea();
2220-
2221- /*Now get the average height*/
2222- Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
2223- Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
2224- surface_input->GetInputAverage(&surface);
2225- base_input->GetInputAverage(&bed);
2226-
2227- /*Return: */
2228- int domaintype;
2229- parameters->FindParam(&domaintype,DomainTypeEnum);
2230- if(domaintype==Domain2DverticalEnum){
2231- return base;
2232- }
2233- else{
2234- return base*(surface-bed);
2235- }
2236+void Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
2237+ TriaRef::GetInputValue(pvalue,values,gauss,P1Enum);
2238 }
2239 /*}}}*/
2240-IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/
2241-
2242- /*The volume above floatation: H + rho_water/rho_ice * bathymetry */
2243- IssmDouble rho_ice,rho_water;
2244- IssmDouble base,surface,bed,bathymetry;
2245- IssmDouble xyz_list[NUMVERTICES][3];
2246-
2247- if(!IsIceInElement() || IsFloating())return 0;
2248-
2249- rho_ice=matpar->GetRhoIce();
2250- rho_water=matpar->GetRhoWater();
2251- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
2252-
2253- /*First calculate the area of the base (cross section triangle)
2254- * http://en.wikipedia.org/wiki/Triangle
2255- * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
2256- 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]));
2257-
2258- /*Now get the average height and bathymetry*/
2259- Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
2260- Input* base_input = inputs->GetInput(BaseEnum); _assert_(base_input);
2261- Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
2262- surface_input->GetInputAverage(&surface);
2263- base_input->GetInputAverage(&bed);
2264- bed_input->GetInputAverage(&bathymetry);
2265-
2266- /*Return: */
2267- return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
2268+int Tria::VelocityInterpolation(void){/*{{{*/
2269+ return TriaRef::VelocityInterpolation(this->element_type);
2270 }
2271 /*}}}*/
2272-IssmDouble Tria::MassFlux( IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
2273-
2274- int domaintype;
2275- IssmDouble mass_flux=0.;
2276- IssmDouble xyz_list[NUMVERTICES][3];
2277- IssmDouble normal[2];
2278- IssmDouble length,rho_ice;
2279- IssmDouble h1,h2;
2280- IssmDouble vx1,vx2,vy1,vy2;
2281- GaussTria* gauss_1=NULL;
2282- GaussTria* gauss_2=NULL;
2283-
2284- /*Get material parameters :*/
2285- rho_ice=matpar->GetRhoIce();
2286-
2287- /*First off, check that this segment belongs to this element: */
2288- if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id);
2289-
2290- /*Get xyz list: */
2291- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
2292-
2293- /*get area coordinates of 0 and 1 locations: */
2294- gauss_1=new GaussTria();
2295- gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
2296- gauss_2=new GaussTria();
2297- gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
2298-
2299- normal[0]=cos(atan2(x1-x2,y2-y1));
2300- normal[1]=sin(atan2(x1-x2,y2-y1));
2301-
2302- length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
2303-
2304- Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
2305- this->parameters->FindParam(&domaintype,DomainTypeEnum);
2306- Input* vx_input=NULL;
2307- Input* vy_input=NULL;
2308- if(domaintype==Domain2DhorizontalEnum){
2309- vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2310- vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2311- }
2312- else{
2313- vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
2314- vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
2315- }
2316-
2317- thickness_input->GetInputValue(&h1, gauss_1);
2318- thickness_input->GetInputValue(&h2, gauss_2);
2319- vx_input->GetInputValue(&vx1,gauss_1);
2320- vx_input->GetInputValue(&vx2,gauss_2);
2321- vy_input->GetInputValue(&vy1,gauss_1);
2322- vy_input->GetInputValue(&vy2,gauss_2);
2323-
2324- mass_flux= rho_ice*length*(
2325- (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
2326- (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
2327- );
2328-
2329- /*clean up and return:*/
2330- delete gauss_1;
2331- delete gauss_2;
2332- return mass_flux;
2333+int Tria::VertexConnectivity(int vertexindex){/*{{{*/
2334+ _assert_(this->vertices);
2335+ return this->vertices[vertexindex]->Connectivity();
2336 }
2337 /*}}}*/
2338-IssmDouble Tria::MassFlux( IssmDouble* segment){/*{{{*/
2339+void Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
2340
2341- int domaintype;
2342- IssmDouble mass_flux=0.;
2343- IssmDouble xyz_list[NUMVERTICES][3];
2344- IssmDouble normal[2];
2345- IssmDouble length,rho_ice;
2346- IssmDouble x1,y1,x2,y2,h1,h2;
2347- IssmDouble vx1,vx2,vy1,vy2;
2348- GaussTria* gauss_1=NULL;
2349- GaussTria* gauss_2=NULL;
2350+ int normal_orientation=0;
2351+ IssmDouble s1,s2;
2352+ IssmDouble levelset[NUMVERTICES];
2353
2354- /*Get material parameters :*/
2355- rho_ice=matpar->GetRhoIce();
2356+ /*Recover parameters and values*/
2357+ IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
2358+ GetInputListOnVertices(&levelset[0],levelsetenum);
2359
2360- /*First off, check that this segment belongs to this element: */
2361- if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
2362+ if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
2363+ /*Portion of the segments*/
2364+ s1=levelset[2]/(levelset[2]-levelset[1]);
2365+ s2=levelset[2]/(levelset[2]-levelset[0]);
2366
2367- /*Recover segment node locations: */
2368- x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
2369+ if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle depending on distribution of levelsetfunction
2370+ /*New point 1*/
2371+ xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
2372+ xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
2373+ xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
2374
2375- /*Get xyz list: */
2376- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
2377-
2378- /*get area coordinates of 0 and 1 locations: */
2379- gauss_1=new GaussTria();
2380- gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
2381- gauss_2=new GaussTria();
2382- gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
2383-
2384- normal[0]=cos(atan2(x1-x2,y2-y1));
2385- normal[1]=sin(atan2(x1-x2,y2-y1));
2386-
2387- length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
2388-
2389- Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
2390- this->parameters->FindParam(&domaintype,DomainTypeEnum);
2391- Input* vx_input=NULL;
2392- Input* vy_input=NULL;
2393- if(domaintype==Domain2DhorizontalEnum){
2394- vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
2395- vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
2396+ /*New point 0*/
2397+ xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
2398+ xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
2399+ xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
2400 }
2401- else{
2402- vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
2403- vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
2404- }
2405+ else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
2406+ /*Portion of the segments*/
2407+ s1=levelset[0]/(levelset[0]-levelset[2]);
2408+ s2=levelset[0]/(levelset[0]-levelset[1]);
2409
2410- thickness_input->GetInputValue(&h1, gauss_1);
2411- thickness_input->GetInputValue(&h2, gauss_2);
2412- vx_input->GetInputValue(&vx1,gauss_1);
2413- vx_input->GetInputValue(&vx2,gauss_2);
2414- vy_input->GetInputValue(&vy1,gauss_1);
2415- vy_input->GetInputValue(&vy2,gauss_2);
2416+ if(levelset[0]<0.) normal_orientation=1;
2417+ /*New point 1*/
2418+ xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
2419+ xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
2420+ xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
2421
2422- mass_flux= rho_ice*length*(
2423- (ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
2424- (ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
2425- );
2426-
2427- /*clean up and return:*/
2428- delete gauss_1;
2429- delete gauss_2;
2430- return mass_flux;
2431-}
2432-/*}}}*/
2433-void Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
2434-
2435- switch(response_enum){
2436- case MaterialsRheologyBbarEnum:
2437- *presponse=this->material->GetBbar();
2438- break;
2439-
2440- case VelEnum:{
2441-
2442- /*Get input:*/
2443- IssmDouble vel;
2444- Input* vel_input;
2445-
2446- vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
2447- vel_input->GetInputAverage(&vel);
2448-
2449- /*Assign output pointers:*/
2450- *presponse=vel;}
2451- break;
2452- default:
2453- _error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
2454+ /*New point 2*/
2455+ xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
2456+ xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
2457+ xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
2458 }
2459+ else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
2460+ /*Portion of the segments*/
2461+ s1=levelset[1]/(levelset[1]-levelset[0]);
2462+ s2=levelset[1]/(levelset[1]-levelset[2]);
2463
2464-}
2465-/*}}}*/
2466-IssmDouble Tria::TotalSmb(void){/*{{{*/
2467+ if(levelset[1]<0.) normal_orientation=1;
2468+ /*New point 0*/
2469+ xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
2470+ xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
2471+ xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
2472
2473- /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
2474- IssmDouble base,smb,rho_ice;
2475- IssmDouble Total_Smb=0;
2476- IssmDouble xyz_list[NUMVERTICES][3];
2477-
2478- /*Get material parameters :*/
2479- rho_ice=matpar->GetRhoIce();
2480-
2481- if(!IsIceInElement())return 0;
2482-
2483- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
2484-
2485- /*First calculate the area of the base (cross section triangle)
2486- * http://en.wikipedia.org/wiki/Triangle
2487- * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
2488- base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
2489-
2490- /*Now get the average SMB over the element*/
2491- Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
2492- smb_input->GetInputAverage(&smb); // average smb on element in m ice s-1
2493- Total_Smb=rho_ice*base*smb; // smb on element in kg s-1
2494-
2495- /*Return: */
2496- return Total_Smb;
2497-}
2498-/*}}}*/
2499-IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/
2500-
2501- /*Intermediaries*/
2502- IssmDouble weight;
2503- IssmDouble Jdet;
2504- IssmDouble Jelem = 0;
2505- IssmDouble xyz_list[NUMVERTICES][3];
2506- GaussTria *gauss = NULL;
2507-
2508- /*If on water, return 0: */
2509- if(!IsIceInElement())return 0;
2510-
2511- /*Retrieve all inputs we will be needing: */
2512- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
2513- Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input);
2514-
2515- /* Start looping on the number of gaussian points: */
2516- gauss=new GaussTria(2);
2517- for(int ig=gauss->begin();ig<gauss->end();ig++){
2518-
2519- gauss->GaussPoint(ig);
2520-
2521- /* Get Jacobian determinant: */
2522- GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
2523-
2524- /*Get parameters at gauss point*/
2525- weights_input->GetInputValue(&weight,gauss);
2526-
2527- /*compute misfit between model and observation */
2528- Jelem+=Jdet*weight*gauss->weight;
2529+ /*New point 2*/
2530+ xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
2531+ xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
2532+ xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
2533 }
2534+ else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
2535+ xyz_zero[3*0+0]=xyz_list[0*3+0];
2536+ xyz_zero[3*0+1]=xyz_list[0*3+1];
2537+ xyz_zero[3*0+2]=xyz_list[0*3+2];
2538
2539- /* clean up and Return: */
2540- delete gauss;
2541- return Jelem;
2542-}
2543-/*}}}*/
2544-IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/
2545-
2546- /*Intermediaries*/
2547- IssmDouble model,observation,weight;
2548- IssmDouble Jdet;
2549- IssmDouble Jelem = 0;
2550- IssmDouble xyz_list[NUMVERTICES][3];
2551- GaussTria *gauss = NULL;
2552-
2553- /*If on water, return 0: */
2554- if(!IsIceInElement())return 0;
2555-
2556- /*Retrieve all inputs we will be needing: */
2557- ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
2558- Input* model_input=inputs->GetInput(modelenum); _assert_(model_input);
2559- Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input);
2560- Input* weights_input =inputs->GetInput(weightsenum); _assert_(weights_input);
2561-
2562- /* Start looping on the number of gaussian points: */
2563- gauss=new GaussTria(2);
2564- for(int ig=gauss->begin();ig<gauss->end();ig++){
2565-
2566- gauss->GaussPoint(ig);
2567-
2568- /* Get Jacobian determinant: */
2569- GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
2570-
2571- /*Get parameters at gauss point*/
2572- model_input->GetInputValue(&model,gauss);
2573- observation_input->GetInputValue(&observation,gauss);
2574- weights_input->GetInputValue(&weight,gauss);
2575-
2576- /*compute misfit between model and observation */
2577- Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight;
2578+ /*New point 2*/
2579+ xyz_zero[3*1+0]=xyz_list[1*3+0];
2580+ xyz_zero[3*1+1]=xyz_list[1*3+1];
2581+ xyz_zero[3*1+2]=xyz_list[1*3+2];
2582 }
2583+ else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
2584+ xyz_zero[3*0+0]=xyz_list[2*3+0];
2585+ xyz_zero[3*0+1]=xyz_list[2*3+1];
2586+ xyz_zero[3*0+2]=xyz_list[2*3+2];
2587
2588- /* clean up and Return: */
2589- delete gauss;
2590- return Jelem;
2591-}
2592-/*}}}*/
2593-IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
2594-
2595-
2596- /*intermediary: */
2597- IssmDouble* values=NULL;
2598- Input* thickness_input=NULL;
2599- IssmDouble thickness;
2600- IssmDouble weight;
2601- IssmDouble Jdet;
2602- IssmDouble volume;
2603- IssmDouble rho_ice;
2604- IssmDouble* xyz_list=NULL;
2605- int point1;
2606- IssmDouble fraction1,fraction2;
2607- bool mainlynegative=true;
2608-
2609- /*Output:*/
2610- volume=0;
2611-
2612- /* Get node coordinates and dof list: */
2613- GetVerticesCoordinates(&xyz_list);
2614-
2615- /*Retrieve inputs required:*/
2616- thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
2617-
2618- /*Retrieve material parameters: */
2619- rho_ice=matpar->GetRhoIce();
2620-
2621- /*Retrieve values of the levelset defining the masscon: */
2622- values = xNew<IssmDouble>(NUMVERTICES);
2623- for(int i=0;i<NUMVERTICES;i++){
2624- values[i]=levelset[this->vertices[i]->Sid()];
2625+ /*New point 2*/
2626+ xyz_zero[3*1+0]=xyz_list[0*3+0];
2627+ xyz_zero[3*1+1]=xyz_list[0*3+1];
2628+ xyz_zero[3*1+2]=xyz_list[0*3+2];
2629 }
2630-
2631- /*Ok, use the level set values to figure out where we put our gaussian points:*/
2632- this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
2633- Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
2634+ else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
2635+ xyz_zero[3*0+0]=xyz_list[1*3+0];
2636+ xyz_zero[3*0+1]=xyz_list[1*3+1];
2637+ xyz_zero[3*0+2]=xyz_list[1*3+2];
2638
2639- volume=0;
2640-
2641- for(int ig=gauss->begin();ig<gauss->end();ig++){
2642- gauss->GaussPoint(ig);
2643-
2644- this->JacobianDeterminant(&Jdet,xyz_list,gauss);
2645- thickness_input->GetInputValue(&thickness, gauss);
2646-
2647- volume+=thickness*gauss->weight*Jdet;
2648+ /*New point 2*/
2649+ xyz_zero[3*1+0]=xyz_list[2*3+0];
2650+ xyz_zero[3*1+1]=xyz_list[2*3+1];
2651+ xyz_zero[3*1+2]=xyz_list[2*3+2];
2652 }
2653+ else _error_("Case not covered");
2654
2655- /* clean up and Return: */
2656- xDelete<IssmDouble>(xyz_list);
2657- xDelete<IssmDouble>(values);
2658- delete gauss;
2659- return rho_ice*volume;
2660+ /*Assign output pointer*/
2661+ *pxyz_zero= xyz_zero;
2662 }
2663 /*}}}*/
2664
2665@@ -2958,179 +3175,45 @@
2666 /*}}}*/
2667 #endif
2668
2669-void Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
2670+#ifdef _HAVE_DAKOTA_
2671+void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
2672
2673- /*Intermediary*/
2674- int num_controls;
2675- int* control_type=NULL;
2676- Input* input=NULL;
2677+ int i,t,row;
2678+ IssmDouble time;
2679+ TransientInput *transientinput = NULL;
2680+ IssmDouble values[3];
2681
2682- /*retrieve some parameters: */
2683- this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
2684- this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
2685+ /*Check that name is an element input*/
2686+ if (!IsInput(name)) return;
2687
2688- for(int i=0;i<num_controls;i++){
2689- input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
2690- if (input->ObjectEnum()!=ControlInputEnum){
2691- _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
2692- }
2693+ switch(type){
2694
2695- ((ControlInput*)input)->UpdateValue(scalar);
2696- ((ControlInput*)input)->Constrain();
2697- if (save_parameter) ((ControlInput*)input)->SaveValue();
2698+ case VertexEnum:
2699+ /*Create transient input: */
2700+ for(t=0;t<ncols;t++){ //ncols is the number of times
2701
2702- }
2703+ /*create input values: */
2704+ for(i=0;i<3;i++){
2705+ row=this->vertices[i]->Sid();
2706+ values[i]=matrix[ncols*row+t];
2707+ }
2708
2709- /*Clean up and return*/
2710- xDelete<int>(control_type);
2711-}
2712-/*}}}*/
2713-void Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
2714+ /*time:*/
2715+ time=matrix[(nrows-1)*ncols+t];
2716
2717- int vertexpidlist[NUMVERTICES];
2718- IssmDouble grad_list[NUMVERTICES];
2719- Input* grad_input=NULL;
2720+ if(t==0) transientinput=new TransientInput(name);
2721+ transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time);
2722+ transientinput->Configure(parameters);
2723+ }
2724+ this->inputs->AddInput(transientinput);
2725+ break;
2726
2727- Input* input=inputs->GetInput(enum_type);
2728- if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
2729- if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
2730-
2731- GradientIndexing(&vertexpidlist[0],control_index);
2732- for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
2733- grad_input=new TriaInput(GradientEnum,grad_list,P1Enum);
2734-
2735- ((ControlInput*)input)->SetGradient(grad_input);
2736-
2737-}/*}}}*/
2738-void Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
2739-
2740- Input* input=inputs->GetInput(control_enum);
2741- if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
2742- if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
2743-
2744- int sidlist[NUMVERTICES];
2745- int connectivity[NUMVERTICES];
2746- IssmPDouble values[NUMVERTICES];
2747- IssmPDouble gradients[NUMVERTICES];
2748- IssmDouble value,gradient;
2749-
2750- this->GetVerticesConnectivityList(&connectivity[0]);
2751- this->GetVerticesSidList(&sidlist[0]);
2752-
2753- GaussTria* gauss=new GaussTria();
2754- for (int iv=0;iv<NUMVERTICES;iv++){
2755- gauss->GaussVertex(iv);
2756-
2757- ((ControlInput*)input)->GetInputValue(&value,gauss);
2758- ((ControlInput*)input)->GetGradientValue(&gradient,gauss);
2759-
2760- values[iv] = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
2761- gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
2762+ default:
2763+ _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
2764 }
2765- delete gauss;
2766
2767- vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
2768- vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
2769-
2770-}/*}}}*/
2771-void Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
2772-
2773- int vertexidlist[NUMVERTICES];
2774- Input *input=NULL;
2775-
2776- /*Get out if this is not an element input*/
2777- if(!IsInput(control_enum)) return;
2778-
2779- /*Prepare index list*/
2780- GradientIndexing(&vertexidlist[0],control_index,onsid);
2781-
2782- /*Get input (either in element or material)*/
2783- input=(Input*)this->inputs->GetInput(control_enum); _assert_(input);
2784-
2785- /*Check that it is a ControlInput*/
2786- if (input->ObjectEnum()!=ControlInputEnum){
2787- _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
2788- }
2789-
2790- ((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
2791 }
2792 /*}}}*/
2793-void Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
2794-
2795- IssmDouble values[NUMVERTICES];
2796- int vertexpidlist[NUMVERTICES],control_init;
2797-
2798-
2799- /*Get Domain type*/
2800- int domaintype;
2801- parameters->FindParam(&domaintype,DomainTypeEnum);
2802-
2803- /*Specific case for depth averaged quantities*/
2804- control_init=control_enum;
2805- if(domaintype==Domain2DverticalEnum){
2806- if(control_enum==MaterialsRheologyBbarEnum){
2807- control_enum=MaterialsRheologyBEnum;
2808- if(!IsOnBase()) return;
2809- }
2810- if(control_enum==DamageDbarEnum){
2811- control_enum=DamageDEnum;
2812- if(!IsOnBase()) return;
2813- }
2814- }
2815-
2816- /*Get out if this is not an element input*/
2817- if(!IsInput(control_enum)) return;
2818-
2819- /*Prepare index list*/
2820- GradientIndexing(&vertexpidlist[0],control_index);
2821-
2822- /*Get values on vertices*/
2823- for(int i=0;i<NUMVERTICES;i++){
2824- values[i]=vector[vertexpidlist[i]];
2825- }
2826- Input* new_input = new TriaInput(control_enum,values,P1Enum);
2827- Input* input = (Input*)this->inputs->GetInput(control_enum); _assert_(input);
2828- if(input->ObjectEnum()!=ControlInputEnum){
2829- _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
2830- }
2831-
2832- ((ControlInput*)input)->SetInput(new_input);
2833-}
2834-/*}}}*/
2835-void Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
2836-
2837- int *doflist = NULL;
2838- IssmDouble value;
2839-
2840- /*Fetch number of nodes for this finite element*/
2841- int numnodes = this->NumberofNodes(this->element_type);
2842-
2843- /*Fetch dof list and allocate solution vector*/
2844- GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
2845- IssmDouble* values = xNew<IssmDouble>(numnodes);
2846-
2847- /*Get inputs*/
2848- Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
2849-
2850- /*Ok, we have the values, fill in the array: */
2851- GaussTria* gauss=new GaussTria();
2852- for(int i=0;i<numnodes;i++){
2853- gauss->GaussNode(this->element_type,i);
2854-
2855- enum_input->GetInputValue(&value,gauss);
2856- values[i]=value;
2857- }
2858-
2859- solution->SetValues(numnodes,doflist,values,INS_VAL);
2860-
2861- /*Free ressources:*/
2862- xDelete<int>(doflist);
2863- xDelete<IssmDouble>(values);
2864- delete gauss;
2865-}
2866-/*}}}*/
2867-
2868-#ifdef _HAVE_DAKOTA_
2869 void Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
2870
2871 int i,j;
2872@@ -3222,89 +3305,4 @@
2873
2874 }
2875 /*}}}*/
2876-void Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
2877-
2878- int i,t,row;
2879- IssmDouble time;
2880- TransientInput *transientinput = NULL;
2881- IssmDouble values[3];
2882-
2883- /*Check that name is an element input*/
2884- if (!IsInput(name)) return;
2885-
2886- switch(type){
2887-
2888- case VertexEnum:
2889- /*Create transient input: */
2890- for(t=0;t<ncols;t++){ //ncols is the number of times
2891-
2892- /*create input values: */
2893- for(i=0;i<3;i++){
2894- row=this->vertices[i]->Sid();
2895- values[i]=matrix[ncols*row+t];
2896- }
2897-
2898- /*time:*/
2899- time=matrix[(nrows-1)*ncols+t];
2900-
2901- if(t==0) transientinput=new TransientInput(name);
2902- transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time);
2903- transientinput->Configure(parameters);
2904- }
2905- this->inputs->AddInput(transientinput);
2906- break;
2907-
2908- default:
2909- _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
2910- }
2911-
2912-}
2913-/*}}}*/
2914 #endif
2915-
2916-void Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
2917-
2918- IssmDouble h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
2919- IssmDouble bed_hydro;
2920- IssmDouble rho_water,rho_ice,density;
2921-
2922- /*material parameters: */
2923- rho_water=matpar->GetRhoWater();
2924- rho_ice=matpar->GetRhoIce();
2925- density=rho_ice/rho_water;
2926- GetInputListOnVertices(&h[0],ThicknessEnum);
2927- GetInputListOnVertices(&r[0],BedEnum);
2928- GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
2929-
2930- /*go through vertices, and figure out which ones are grounded and want to unground: */
2931- for(int i=0;i<NUMVERTICES;i++){
2932- /*Find if grounded vertices want to start floating*/
2933- if (gl[i]>0.){
2934- bed_hydro=-density*h[i];
2935- if(bed_hydro>r[i]){
2936- /*Vertex that could potentially unground, flag it*/
2937- potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
2938- }
2939- }
2940- }
2941-}
2942-/*}}}*/
2943-int Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
2944-
2945- int i;
2946- int nflipped=0;
2947-
2948- /*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
2949- for(i=0;i<3;i++){
2950- if (reCast<bool>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
2951- vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
2952-
2953- /*If node was not on ice shelf, we flipped*/
2954- if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
2955- nflipped++;
2956- }
2957- }
2958- }
2959- return nflipped;
2960-}
2961-/*}}}*/
2962Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
2963===================================================================
2964--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18910)
2965+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 18911)
2966@@ -41,105 +41,101 @@
2967 Object *copy();
2968 /*}}}*/
2969 /*Update virtual functions resolution: {{{*/
2970- void InputUpdateFromVector(IssmDouble* vector, int name, int type);
2971 #ifdef _HAVE_DAKOTA_
2972+ void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
2973 void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
2974- void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
2975 #endif
2976 void InputUpdateFromIoModel(int index, IoModel* iomodel);
2977+ void InputUpdateFromVector(IssmDouble* vector, int name, int type);
2978 /*}}}*/
2979 /*Element virtual functions definitions: {{{*/
2980+ void AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
2981+ void CalvingRateLevermann();
2982 IssmDouble CharacteristicLength(void);
2983 void ComputeBasalStress(Vector<IssmDouble>* sigma_b);
2984+ void ComputeDeviatoricStressTensor();
2985 void ComputeSigmaNN();
2986 void ComputeStressTensor();
2987- void ComputeDeviatoricStressTensor();
2988- void StrainRateparallel();
2989- void StrainRateperpendicular();
2990 void ComputeSurfaceNormalVelocity();
2991- void StressIntensityFactor(void){_error_("not implemented yet");};
2992- void CalvingRateLevermann();
2993 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
2994- void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
2995- void ResetHooks();
2996+ void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
2997+ void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
2998 void Delta18oParameterization(void);
2999+ int EdgeOnBaseIndex();
3000+ void EdgeOnBaseIndices(int* pindex1,int* pindex);
3001+ int EdgeOnSurfaceIndex();
3002+ void EdgeOnSurfaceIndices(int* pindex1,int* pindex);
3003+ void ElementResponse(IssmDouble* presponse,int response_enum);
3004 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
3005+ int FiniteElement(void);
3006 void FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
3007- int FiniteElement(void);
3008- Element* GetUpperElement(void){_error_("not implemented yet");};
3009 Element* GetBasalElement(void){_error_("not implemented yet");};
3010 void GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues);
3011 void GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
3012 IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
3013+ void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
3014+ void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
3015 int GetNodeIndex(Node* node);
3016 int GetNumberOfNodes(void);
3017 int GetNumberOfNodes(int enum_type);
3018 int GetNumberOfVertices(void);
3019- bool IsOnBase();
3020- bool IsOnSurface();
3021- bool HasEdgeOnBase();
3022- bool HasEdgeOnSurface();
3023- void EdgeOnSurfaceIndices(int* pindex1,int* pindex);
3024- void EdgeOnBaseIndices(int* pindex1,int* pindex);
3025- int EdgeOnBaseIndex();
3026- int EdgeOnSurfaceIndex();
3027- bool IsNodeOnShelfFromFlags(IssmDouble* flags);
3028- int NumberofNodesVelocity(void);
3029- int NumberofNodesPressure(void);
3030 void GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
3031+ Element* GetUpperElement(void){_error_("not implemented yet");};
3032+ void GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
3033 void GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
3034 void GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
3035+ bool HasEdgeOnBase();
3036+ bool HasEdgeOnSurface();
3037+ IssmDouble IceVolume(void);
3038+ IssmDouble IceVolumeAboveFloatation(void);
3039+ void InputControlUpdate(IssmDouble scalar,bool save_parameter);
3040 void InputDepthAverageAtBase(int enum_type,int average_enum_type);
3041 void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
3042 void InputScale(int enum_type,IssmDouble scale_factor);
3043+ bool IsFaceOnBoundary(void);
3044+ bool IsIcefront(void);
3045+ bool IsNodeOnShelfFromFlags(IssmDouble* flags);
3046+ bool IsOnBase();
3047+ bool IsOnSurface();
3048+ bool IsZeroLevelset(int levelset_enum);
3049+ IssmDouble Masscon(IssmDouble* levelset);
3050+ IssmDouble MassFlux(IssmDouble* segment);
3051+ IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
3052 void MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
3053+ IssmDouble Misfit(int modelenum,int observationenum,int weightsenum);
3054+ IssmDouble MisfitArea(int weightsenum);
3055 int NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
3056+ int NumberofNodesPressure(void);
3057+ int NumberofNodesVelocity(void);
3058 void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
3059+ void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
3060+ int PressureInterpolation();
3061 void ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
3062 void ResetFSBasalBoundaryCondition(void);
3063+ void ResetHooks();
3064+ void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
3065+ void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
3066 Element* SpawnBasalElement(void);
3067 Element* SpawnTopElement(void);
3068- int VelocityInterpolation();
3069- int PressureInterpolation();
3070+ void StrainRateparallel();
3071+ void StrainRateperpendicular();
3072+ void StressIntensityFactor(void){_error_("not implemented yet");};
3073+ IssmDouble SurfaceArea(void);
3074 int TensorInterpolation();
3075- IssmDouble SurfaceArea(void);
3076+ IssmDouble TimeAdapt();
3077+ IssmDouble TotalSmb(void);
3078 void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
3079- IssmDouble TimeAdapt();
3080+ int UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
3081+ void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
3082 void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
3083- void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
3084+ int VelocityInterpolation();
3085 int VertexConnectivity(int vertexindex);
3086 void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
3087 void ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
3088- void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
3089- void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
3090- bool IsZeroLevelset(int levelset_enum);
3091- bool IsIcefront(void);
3092- bool IsFaceOnBoundary(void);
3093
3094- void AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
3095- IssmDouble IceVolume(void);
3096- IssmDouble IceVolumeAboveFloatation(void);
3097- IssmDouble TotalSmb(void);
3098- IssmDouble MassFlux(IssmDouble* segment);
3099- IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
3100- void ElementResponse(IssmDouble* presponse,int response_enum);
3101- IssmDouble Masscon(IssmDouble* levelset);
3102- IssmDouble Misfit(int modelenum,int observationenum,int weightsenum);
3103- IssmDouble MisfitArea(int weightsenum);
3104-
3105 #ifdef _HAVE_GIA_
3106 void GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
3107 #endif
3108-
3109- void GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
3110- void SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
3111- void ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
3112- void ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
3113- void InputControlUpdate(IssmDouble scalar,bool save_parameter);
3114-
3115- void PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
3116- int UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
3117-
3118 /*}}}*/
3119 /*Tria specific routines:{{{*/
3120 void AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
3121@@ -147,18 +143,15 @@
3122 IssmDouble GetArea(void);
3123 void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
3124 int GetElementType(void);
3125- void NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
3126- void NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
3127- void NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
3128 void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
3129 void GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
3130 Node* GetNode(int node_number);
3131 void InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
3132 void InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
3133 void JacobianDeterminant(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss);
3134+ void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
3135 void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
3136 void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss);
3137- void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
3138 void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
3139 IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
3140 Gauss* NewGauss(void);
3141@@ -170,23 +163,25 @@
3142 Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
3143 Gauss* NewGaussTop(int order);
3144 void NodalFunctions(IssmDouble* basis,Gauss* gauss);
3145- void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
3146- void NodalFunctionsP2(IssmDouble* basis,Gauss* gauss);
3147 void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
3148- void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
3149+ void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
3150 void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
3151- void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
3152- void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
3153 void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
3154+ void NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
3155+ void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
3156+ void NodalFunctionsP2(IssmDouble* basis,Gauss* gauss);
3157 void NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss);
3158+ void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
3159+ void NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
3160+ void NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
3161+ void NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
3162 void SetClone(int* minranks);
3163 void SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
3164 Seg* SpawnSeg(int index1,int index2);
3165 IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
3166+ void UpdateConstraintsExtrudeFromBase(void);
3167+ void UpdateConstraintsExtrudeFromTop(void);
3168 void ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
3169-
3170- void UpdateConstraintsExtrudeFromBase(void);
3171- void UpdateConstraintsExtrudeFromTop(void);
3172 /*}}}*/
3173
3174 };
Note: See TracBrowser for help on using the repository browser.