source: issm/oecreview/Archive/16133-16554/ISSM-16349-16350.diff@ 16556

Last change on this file since 16556 was 16556, checked in by Mathieu Morlighem, 11 years ago

NEW: added Archive/16133-16554
File size: 36.7 KB
RevLine 
[16556]1Index: ../trunk-jpl/src/c/classes/Materials/Matice.cpp
2===================================================================
3--- ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16349)
4+++ ../trunk-jpl/src/c/classes/Materials/Matice.cpp (revision 16350)
5@@ -291,6 +291,72 @@
6 *pviscosity=viscosity;
7 }
8 /*}}}*/
9+/*FUNCTION Matice::GetViscosity2dvertical {{{*/
10+void Matice::GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* epsilon){
11+ /*From a string tensor and a material object, return viscosity, using Glen's flow law.
12+ (1-D) B
13+ viscosity= --------------------------------------
14+ 2[ exx^2+eyy^2+ 2exy^2]^[(n-1)/2n]
15+
16+ where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity
17+ vector, and n the flow law exponent.
18+
19+ If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we
20+ return 10^14, initial viscosity.
21+ */
22+
23+ /*output: */
24+ IssmDouble viscosity;
25+
26+ /*input strain rate: */
27+ IssmDouble exx,eyy,exy;
28+
29+ /*Intermediary: */
30+ IssmDouble A,e;
31+ IssmDouble B,D,n;
32+
33+ /*Get B and n*/
34+ B=GetB();
35+ n=GetN();
36+ D=GetD();
37+
38+ if (n==1){
39+ /*Viscous behaviour! viscosity=B: */
40+ viscosity=(1-D)*B/2;
41+ }
42+ else{
43+ if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
44+ viscosity=0.5*pow(10.,14);
45+ }
46+ else{
47+ /*Retrive strain rate components: */
48+ exx=epsilon[0];
49+ eyy=epsilon[1];
50+ exy=epsilon[2];
51+
52+ /*Build viscosity: viscosity=B/(2*A^e) */
53+ A=exx*exx+eyy*eyy+2.*exy*exy;
54+ if(A==0.){
55+ /*Maxiviscositym viscosity for 0 shear areas: */
56+ viscosity=2.5*2.e+17;
57+ }
58+ else{
59+ e=(n-1.)/(2.*n);
60+ viscosity=(1.-D)*B/(2.*pow(A,e));
61+ }
62+ }
63+ }
64+
65+ /*Checks in debugging mode*/
66+ if(viscosity<=0) _error_("Negative viscosity");
67+ _assert_(B>0);
68+ _assert_(n>0);
69+ _assert_(D>=0 && D<1);
70+
71+ /*Return: */
72+ *pviscosity=viscosity;
73+}
74+/*}}}*/
75 /*FUNCTION Matice::GetViscosity3d {{{*/
76 void Matice::GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* epsilon){
77
78@@ -792,10 +858,10 @@
79 if(iomodel->meshtype==Mesh2DhorizontalEnum){
80
81 /*Intermediaries*/
82- const int num_vertices = 3; //Tria has 3 vertices
83- IssmDouble nodeinputs[num_vertices];
84- IssmDouble cmmininputs[num_vertices];
85- IssmDouble cmmaxinputs[num_vertices];
86+ const int num_vertices = 3; //Tria has 3 vertices
87+ IssmDouble nodeinputs[num_vertices];
88+ IssmDouble cmmininputs[num_vertices];
89+ IssmDouble cmmaxinputs[num_vertices];
90
91 /*Get B*/
92 if (iomodel->Data(MaterialsRheologyBEnum)) {
93@@ -866,7 +932,7 @@
94 /*Get D:*/
95 if (iomodel->Data(DamageDEnum)) {
96 for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1];
97- this->inputs->AddInput(new TriaInput(DamageDbarEnum,nodeinputs,P1Enum));
98+ this->inputs->AddInput(new TriaInput(DamageDEnum,nodeinputs,P1Enum));
99 }
100 }
101 #ifdef _HAVE_3D_
102Index: ../trunk-jpl/src/c/classes/Materials/Material.h
103===================================================================
104--- ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16349)
105+++ ../trunk-jpl/src/c/classes/Materials/Material.h (revision 16350)
106@@ -26,6 +26,7 @@
107 virtual void Configure(Elements* elements)=0;
108 virtual void GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum)=0;
109 virtual void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
110+ virtual void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
111 virtual void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0;
112 virtual void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0;
113 virtual void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon)=0;
114Index: ../trunk-jpl/src/c/classes/Materials/Matice.h
115===================================================================
116--- ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16349)
117+++ ../trunk-jpl/src/c/classes/Materials/Matice.h (revision 16350)
118@@ -53,6 +53,7 @@
119 void GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum);
120 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
121 void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
122+ void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon);
123 void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
124 void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon);
125 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
126Index: ../trunk-jpl/src/c/classes/Materials/Matpar.h
127===================================================================
128--- ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16349)
129+++ ../trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16350)
130@@ -79,10 +79,11 @@
131 void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");};
132 /*}}}*/
133 /*Material virtual functions resolution: {{{*/
134- void InputDuplicate(int original_enum,int new_enum);
135- void Configure(Elements* elements);
136- void GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
137+ void InputDuplicate(int original_enum,int new_enum);
138+ void Configure(Elements* elements);
139+ void GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
140 void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
141+ void GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
142 void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");};
143 void GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");};
144 void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon){_error_("not supported");};
145Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
146===================================================================
147--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16349)
148+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16350)
149@@ -19,7 +19,9 @@
150 /*}}}*/
151
152 /*Element macros*/
153-#define NUMVERTICES 3
154+#define NUMVERTICES 3
155+#define NUMVERTICES1D 2
156+
157 /*Constructors/destructor/copy*/
158 /*FUNCTION Tria::Tria(){{{*/
159 Tria::Tria(){
160@@ -218,7 +220,18 @@
161 switch(analysis_type){
162 #ifdef _HAVE_STRESSBALANCE_
163 case StressbalanceAnalysisEnum:
164- return CreateKMatrixStressbalanceSSA();
165+ int approximation;
166+ inputs->GetInputValue(&approximation,ApproximationEnum);
167+ switch(approximation){
168+ case SSAApproximationEnum:
169+ return CreateKMatrixStressbalanceSSA();
170+ case FSApproximationEnum:
171+ return CreateKMatrixStressbalanceFS();
172+ case NoneApproximationEnum:
173+ return NULL;
174+ default:
175+ _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
176+ }
177 break;
178 case StressbalanceSIAAnalysisEnum:
179 return CreateKMatrixStressbalanceSIA();
180@@ -372,7 +385,18 @@
181 switch(analysis_type){
182 #ifdef _HAVE_STRESSBALANCE_
183 case StressbalanceAnalysisEnum:
184- return CreatePVectorStressbalanceSSA();
185+ int approximation;
186+ inputs->GetInputValue(&approximation,ApproximationEnum);
187+ switch(approximation){
188+ case SSAApproximationEnum:
189+ return CreatePVectorStressbalanceSSA();
190+ case FSApproximationEnum:
191+ return CreatePVectorStressbalanceFS();
192+ case NoneApproximationEnum:
193+ return NULL;
194+ default:
195+ _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
196+ }
197 break;
198 case StressbalanceSIAAnalysisEnum:
199 return CreatePVectorStressbalanceSIA();
200@@ -1678,7 +1702,17 @@
201 switch(analysis_type){
202 #ifdef _HAVE_STRESSBALANCE_
203 case StressbalanceAnalysisEnum:
204- InputUpdateFromSolutionStressbalanceHoriz(solution);
205+ int approximation;
206+ inputs->GetInputValue(&approximation,ApproximationEnum);
207+ if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
208+ InputUpdateFromSolutionStressbalanceFS(solution);
209+ }
210+ else if (approximation==SSAApproximationEnum || approximation==SIAApproximationEnum){
211+ InputUpdateFromSolutionStressbalanceHoriz(solution);
212+ }
213+ else{
214+ _error_("approximation not supported yet");
215+ }
216 break;
217 case StressbalanceSIAAnalysisEnum:
218 InputUpdateFromSolutionStressbalanceHoriz(solution);
219@@ -2055,6 +2089,50 @@
220 return onbed;
221 }
222 /*}}}*/
223+/*FUNCTION Tria::HasEdgeOnBed {{{*/
224+bool Tria::HasEdgeOnBed(){
225+
226+ IssmDouble values[NUMVERTICES];
227+ IssmDouble sum;
228+
229+ /*Retrieve all inputs and parameters*/
230+ GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
231+ sum = values[0]+values[1]+values[2];
232+
233+ _assert_(sum==0. || sum==1. || sum==2.);
234+
235+ if(sum>1.){
236+ return true;
237+ }
238+ else{
239+ return false;
240+ }
241+}
242+/*}}}*/
243+/*FUNCTION Tria::EdgeOnBedIndices{{{*/
244+void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){
245+
246+ bool found=false;
247+ IssmDouble values[NUMVERTICES];
248+
249+ /*Retrieve all inputs and parameters*/
250+ GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
251+
252+ for(int i=0;i<NUMVERTICES;i++){
253+ if(values[i]==1.){
254+ if(found){
255+ *pindex2 = i;
256+ return;
257+ }
258+ else{
259+ *pindex1 = i;
260+ }
261+ }
262+ }
263+
264+ _error_("Could not find 2 vertices on bed");
265+}
266+/*}}}*/
267 /*FUNCTION Tria::IsFloating {{{*/
268 bool Tria::IsFloating(){
269
270@@ -3036,6 +3114,112 @@
271 #endif
272
273 #ifdef _HAVE_STRESSBALANCE_
274+/*FUNCTION Tria::CreateKMatrixStressbalanceFS{{{*/
275+ElementMatrix* Tria::CreateKMatrixStressbalanceFS(void){
276+
277+ ElementMatrix* Ke1 = NULL;
278+ ElementMatrix* Ke2 = NULL;
279+ ElementMatrix* Ke = NULL;
280+
281+ /*compute all stiffness matrices for this element*/
282+ Ke1=CreateKMatrixStressbalanceFSViscous();
283+ Ke2=CreateKMatrixStressbalanceFSFriction();
284+ Ke =new ElementMatrix(Ke1,Ke2);
285+
286+ /*clean-up and return*/
287+ delete Ke1;
288+ delete Ke2;
289+ return Ke;
290+
291+}
292+/*}}}*/
293+/*FUNCTION Tria::CreateKMatrixStressbalanceFSViscous {{{*/
294+ElementMatrix* Tria::CreateKMatrixStressbalanceFSViscous(void){
295+
296+ /*Intermediaries */
297+ int i,approximation;
298+ IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
299+ IssmDouble xyz_list[NUMVERTICES][3];
300+ IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
301+ GaussTria *gauss=NULL;
302+
303+ /*Fetch number of nodes and dof for this finite element*/
304+ int vnumnodes = this->NumberofNodesVelocity();
305+ int pnumnodes = this->NumberofNodesPressure();
306+ int numdof = vnumnodes*NDOF2 + pnumnodes*NDOF1;
307+
308+ /*Prepare coordinate system list*/
309+ int* cs_list = xNew<int>(vnumnodes+pnumnodes);
310+ for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
311+ for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
312+
313+ /*Initialize Element matrix and vectors*/
314+ ElementMatrix* Ke = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
315+ IssmDouble* B = xNew<IssmDouble>(5*numdof);
316+ IssmDouble* Bprime = xNew<IssmDouble>(5*numdof);
317+ IssmDouble* D = xNewZeroInit<IssmDouble>(5*5);
318+
319+ /*Retrieve all inputs and parameters*/
320+ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
321+ parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
322+ Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
323+ Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
324+
325+ /* Start looping on the number of gaussian points: */
326+ gauss=new GaussTria(5);
327+ for(int ig=gauss->begin();ig<gauss->end();ig++){
328+
329+ gauss->GaussPoint(ig);
330+
331+ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
332+ GetBFS(B,&xyz_list[0][0],gauss);
333+ GetBprimeFS(Bprime,&xyz_list[0][0],gauss);
334+
335+ this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
336+ material->GetViscosity2dvertical(&viscosity,&epsilon[0]);
337+
338+ D_scalar=gauss->weight*Jdet;
339+ for(i=0;i<3;i++) D[i*5+i] = +D_scalar*2.*viscosity;
340+ for(i=3;i<5;i++) D[i*5+i] = -D_scalar*FSreconditioning;
341+
342+ TripleMultiply(B,5,numdof,1,
343+ D,5,5,0,
344+ Bprime,5,numdof,0,
345+ Ke->values,1);
346+ }
347+
348+ /*Transform Coordinate System*/
349+ TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
350+
351+ /*Clean up and return*/
352+ delete gauss;
353+ xDelete<int>(cs_list);
354+ xDelete<IssmDouble>(B);
355+ xDelete<IssmDouble>(Bprime);
356+ xDelete<IssmDouble>(D);
357+ return Ke;
358+}
359+/*}}}*/
360+/*FUNCTION Tria::CreateKMatrixStressbalanceFSFriction{{{*/
361+ElementMatrix* Tria::CreateKMatrixStressbalanceFSFriction(void){
362+
363+ /*Intermediaries */
364+ int i,j;
365+ int analysis_type,approximation;
366+ IssmDouble alpha2,Jdet2d;
367+ IssmDouble FSreconditioning,viscosity;
368+ IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
369+ IssmDouble xyz_list[NUMVERTICES][3];
370+ IssmDouble xyz_list_seg[NUMVERTICES1D][3];
371+ Friction *friction = NULL;
372+ GaussTria *gauss = NULL;
373+
374+ /*Initialize Element matrix and return if necessary*/
375+ if(IsFloating() || !HasEdgeOnBed()) return NULL;
376+
377+ _error_("STOP");
378+}
379+/*}}}*/
380 /*FUNCTION Tria::CreateKMatrixStressbalanceSSA {{{*/
381 ElementMatrix* Tria::CreateKMatrixStressbalanceSSA(void){
382
383@@ -3225,6 +3409,183 @@
384 return Ke;
385 }
386 /*}}}*/
387+/*FUNCTION Tria::CreatePVectorStressbalanceFS {{{*/
388+ElementVector* Tria::CreatePVectorStressbalanceFS(void){
389+
390+ ElementVector* pe1;
391+ ElementVector* pe2;
392+ ElementVector* pe3;
393+ ElementVector* pe;
394+
395+ /*compute all stiffness matrices for this element*/
396+ pe1=CreatePVectorStressbalanceFSViscous();
397+ pe2=CreatePVectorStressbalanceFSShelf();
398+ pe3=CreatePVectorStressbalanceFSFront();
399+ pe =new ElementVector(pe1,pe2,pe3);
400+
401+ /*clean-up and return*/
402+ delete pe1;
403+ delete pe2;
404+ delete pe3;
405+ return pe;
406+}
407+/*}}}*/
408+/*FUNCTION Tria::CreatePVectorStressbalanceFSFront{{{*/
409+ElementVector* Tria::CreatePVectorStressbalanceFSFront(void){
410+
411+ /*Intermediaries */
412+ int i;
413+ IssmDouble ls[NUMVERTICES];
414+ IssmDouble xyz_list[NUMVERTICES][3];
415+ bool isfront;
416+
417+ /*Retrieve all inputs and parameters*/
418+ GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
419+
420+ /*If the level set is awlays <=0, there is no ice front here*/
421+ isfront = false;
422+ if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
423+ if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
424+ isfront = true;
425+ }
426+ }
427+
428+ /*If no front, return NULL*/
429+ if(!isfront) return NULL;
430+
431+ _error_("STOP");
432+}
433+/*}}}*/
434+/*FUNCTION Tria::CreatePVectorStressbalanceFSViscous {{{*/
435+ElementVector* Tria::CreatePVectorStressbalanceFSViscous(void){
436+
437+ /*Intermediaries*/
438+ int i;
439+ int approximation;
440+ IssmDouble Jdet,gravity,rho_ice;
441+ IssmDouble forcex,forcey;
442+ IssmDouble xyz_list[NUMVERTICES][3];
443+ GaussTria *gauss=NULL;
444+
445+ /*Fetch number of nodes and dof for this finite element*/
446+ int vnumnodes = this->NumberofNodesVelocity();
447+ int pnumnodes = this->NumberofNodesPressure();
448+
449+ /*Prepare coordinate system list*/
450+ int* cs_list = xNew<int>(vnumnodes+pnumnodes);
451+ for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
452+ for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
453+
454+ /*Initialize Element matrix and vectors*/
455+ ElementVector* pe = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
456+ IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
457+
458+ /*Retrieve all inputs and parameters*/
459+ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
460+ Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input);
461+ Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input);
462+ rho_ice = matpar->GetRhoIce();
463+ gravity = matpar->GetG();
464+
465+ /* Start looping on the number of gaussian points: */
466+ gauss=new GaussTria(5);
467+ for(int ig=gauss->begin();ig<gauss->end();ig++){
468+
469+ gauss->GaussPoint(ig);
470+
471+ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
472+ GetNodalFunctionsVelocity(vbasis, gauss);
473+
474+ loadingforcex_input->GetInputValue(&forcex,gauss);
475+ loadingforcey_input->GetInputValue(&forcey,gauss);
476+
477+ for(i=0;i<vnumnodes;i++){
478+ pe->values[i*NDOF2+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
479+ pe->values[i*NDOF2+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
480+ pe->values[i*NDOF2+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
481+ }
482+ }
483+
484+ /*Transform coordinate system*/
485+ TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
486+
487+ /*Clean up and return*/
488+ delete gauss;
489+ xDelete<int>(cs_list);
490+ xDelete<IssmDouble>(vbasis);
491+ return pe;
492+}
493+/*}}}*/
494+/*FUNCTION Tria::CreatePVectorStressbalanceFSShelf{{{*/
495+ElementVector* Tria::CreatePVectorStressbalanceFSShelf(void){
496+
497+ /*Intermediaries*/
498+ int i,j;
499+ int indices[2];
500+ IssmDouble gravity,rho_water,bed,water_pressure;
501+ IssmDouble normal_vel,vx,vy,dt;
502+ IssmDouble xyz_list_seg[NUMVERTICES1D][3];
503+ IssmDouble xyz_list[NUMVERTICES][3];
504+ IssmDouble normal[2];
505+ IssmDouble Jdet;
506+
507+ /*Initialize Element vector and return if necessary*/
508+ if(!HasEdgeOnBed() || !IsFloating()) return NULL;
509+
510+ /*Fetch number of nodes and dof for this finite element*/
511+ int vnumnodes = this->NumberofNodesVelocity();
512+ int pnumnodes = this->NumberofNodesPressure();
513+
514+ /*Prepare coordinate system list*/
515+ int* cs_list = xNew<int>(vnumnodes+pnumnodes);
516+ for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
517+ for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
518+
519+ /*Initialize Element matrix and vectors*/
520+ ElementVector* pe = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
521+ IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
522+
523+ /*Retrieve all inputs and parameters*/
524+ rho_water=matpar->GetRhoWater();
525+ gravity=matpar->GetG();
526+ GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
527+ Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input);
528+
529+ /*Get vertex indices that lie on bed*/
530+ this->EdgeOnBedIndices(&indices[0],&indices[1]);
531+
532+ for(i=0;i<NUMVERTICES1D;i++) for(j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
533+
534+ /* Start looping on the number of gauss 1d (nodes on the bedrock) */
535+ GaussTria* gauss=new GaussTria(indices[0],indices[1],2);
536+ for(int ig=gauss->begin();ig<gauss->end();ig++){
537+
538+ gauss->GaussPoint(ig);
539+
540+ GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
541+ GetNodalFunctionsVelocity(vbasis, gauss);
542+
543+ GetSegmentNormal(&normal[0],xyz_list_seg);
544+ _assert_(normal[1]<0.);
545+ bed_input->GetInputValue(&bed, gauss);
546+ water_pressure=gravity*rho_water*bed;
547+
548+ for(i=0;i<vnumnodes;i++){
549+ pe->values[2*i+0]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[0];
550+ pe->values[2*i+1]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[1];
551+ }
552+ }
553+
554+ /*Transform coordinate system*/
555+ TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
556+
557+ /*Clean up and return*/
558+ delete gauss;
559+ xDelete<int>(cs_list);
560+ xDelete<IssmDouble>(vbasis);
561+ return pe;
562+}
563+/*}}}*/
564 /*FUNCTION Tria::CreatePVectorStressbalanceSSA {{{*/
565 ElementVector* Tria::CreatePVectorStressbalanceSSA(){
566
567@@ -3695,6 +4056,79 @@
568
569 }
570 /*}}}*/
571+/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceFS {{{*/
572+void Tria::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){
573+
574+ int i;
575+ int* vdoflist=NULL;
576+ int* pdoflist=NULL;
577+ IssmDouble FSreconditioning;
578+
579+ /*Fetch number of nodes and dof for this finite element*/
580+ int vnumnodes = this->NumberofNodesVelocity();
581+ int pnumnodes = this->NumberofNodesPressure();
582+ int vnumdof = vnumnodes*NDOF2;
583+ int pnumdof = pnumnodes*NDOF1;
584+
585+ /*Initialize values*/
586+ IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
587+ IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
588+ IssmDouble* vx = xNew<IssmDouble>(vnumnodes);
589+ IssmDouble* vy = xNew<IssmDouble>(vnumnodes);
590+ IssmDouble* vel = xNew<IssmDouble>(vnumnodes);
591+ IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
592+
593+ /*Get dof list: */
594+ GetDofListVelocity(&vdoflist,GsetEnum);
595+ GetDofListPressure(&pdoflist,GsetEnum);
596+
597+ /*Use the dof list to index into the solution vector: */
598+ for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
599+ for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
600+
601+ /*Transform solution in Cartesian Space*/
602+ TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYEnum);
603+
604+ /*Ok, we have vx and vy in values, fill in all arrays: */
605+ for(i=0;i<vnumnodes;i++){
606+ vx[i] = vvalues[i*NDOF2+0];
607+ vy[i] = vvalues[i*NDOF2+1];
608+ if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
609+ if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
610+ }
611+ for(i=0;i<pnumnodes;i++){
612+ pressure[i] = pvalues[i];
613+ if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
614+ }
615+
616+ /*Recondition pressure and compute vel: */
617+ this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
618+ for(i=0;i<pnumnodes;i++) pressure[i]=pressure[i]*FSreconditioning;
619+ for(i=0;i<vnumnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
620+
621+ /*Now, we have to move the previous inputs to old
622+ * status, otherwise, we'll wipe them off: */
623+ this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
624+ this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
625+ this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
626+
627+ /*Add vx and vy as inputs to the tria element: */
628+ this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
629+ this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
630+ this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
631+ this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum));
632+
633+ /*Free ressources:*/
634+ xDelete<IssmDouble>(pressure);
635+ xDelete<IssmDouble>(vel);
636+ xDelete<IssmDouble>(vy);
637+ xDelete<IssmDouble>(vx);
638+ xDelete<IssmDouble>(vvalues);
639+ xDelete<IssmDouble>(pvalues);
640+ xDelete<int>(vdoflist);
641+ xDelete<int>(pdoflist);
642+}
643+/*}}}*/
644 /*FUNCTION Tria::InputUpdateFromSolutionStressbalanceSIA {{{*/
645 void Tria::InputUpdateFromSolutionStressbalanceSIA(IssmDouble* solution){
646
647Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.h
648===================================================================
649--- ../trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 16349)
650+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.h (revision 16350)
651@@ -22,8 +22,10 @@
652 void SetElementType(int type,int type_counter);
653
654 /*Numerics*/
655+ void GetBFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss);
656 void GetBSSA(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
657 void GetBSSAFS(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss);
658+ void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss);
659 void GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
660 void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
661 void GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
662@@ -35,10 +37,14 @@
663 void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss);
664 void GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussTria* gauss);
665 void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
666+ void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTria* gauss);
667+ void GetNodalFunctionsPressure(IssmDouble* basis, GaussTria* gauss);
668 void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2);
669 void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
670 void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
671 void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss);
672+ void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
673+ void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
674 void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
675 void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
676 void GetInputDerivativeValue(IssmDouble* pp, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss);
677Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
678===================================================================
679--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 16349)
680+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 16350)
681@@ -333,9 +333,9 @@
682
683 /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
684 * For node i, Bvi can be expressed in the actual coordinate system
685- * by: Bvi=[ dh/dx 0 0 ]
686+ * by: Bvi=[ dh/dx 0 0 ]
687 * [ 0 dh/dy 0 ]
688- * [ 0 0 dh/dy ]
689+ * [ 0 0 dh/dz ]
690 * [ 1/2*dh/dy 1/2*dh/dx 0 ]
691 * [ 1/2*dh/dz 0 1/2*dh/dx ]
692 * [ 0 1/2*dh/dz 1/2*dh/dy ]
693Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
694===================================================================
695--- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16349)
696+++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 16350)
697@@ -83,6 +83,8 @@
698 int GetNumberOfNodes(void);
699 int Sid();
700 bool IsOnBed();
701+ bool HasEdgeOnBed();
702+ void EdgeOnBedIndices(int* pindex1,int* pindex);
703 bool IsFloating();
704 bool IsNodeOnShelfFromFlags(IssmDouble* flags);
705 bool NoIceInElement();
706@@ -232,15 +234,24 @@
707 ElementMatrix* CreateKMatrixStressbalanceSSAViscous(void);
708 ElementMatrix* CreateKMatrixStressbalanceSSAFriction(void);
709 ElementMatrix* CreateKMatrixStressbalanceSIA(void);
710+ ElementMatrix* CreateKMatrixStressbalanceFS(void);
711+ ElementMatrix* CreateKMatrixStressbalanceFSViscous(void);
712+ ElementMatrix* CreateKMatrixStressbalanceFSFriction(void);
713 ElementVector* CreatePVectorStressbalanceSSA(void);
714 ElementVector* CreatePVectorStressbalanceSSADrivingStress(void);
715 ElementVector* CreatePVectorStressbalanceSSAFront(void);
716 ElementVector* CreatePVectorStressbalanceSIA(void);
717+ ElementVector* CreatePVectorStressbalanceFS(void);
718+ ElementVector* CreatePVectorStressbalanceFSFront(void);
719+ ElementVector* CreatePVectorStressbalanceFSViscous(void);
720+ void PVectorGLSstabilization(ElementVector* pe);
721+ ElementVector* CreatePVectorStressbalanceFSShelf(void);
722 ElementMatrix* CreateJacobianStressbalanceSSA(void);
723 void GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
724 void GetSolutionFromInputsStressbalanceHoriz(Vector<IssmDouble>* solution);
725 void GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
726 void InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
727+ void InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
728 void InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
729 #endif
730
731Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
732===================================================================
733--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16349)
734+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16350)
735@@ -10657,7 +10657,7 @@
736 IssmDouble* vy = xNew<IssmDouble>(vnumnodes);
737 IssmDouble* vz = xNew<IssmDouble>(vnumnodes);
738 IssmDouble* vel = xNew<IssmDouble>(vnumnodes);
739- IssmDouble* pressure = xNew<IssmDouble>(vnumnodes);
740+ IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
741
742 /*Get dof list: */
743 GetDofListVelocity(&vdoflist,GsetEnum);
744Index: ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp
745===================================================================
746--- ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 16349)
747+++ ../trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 16350)
748@@ -199,6 +199,116 @@
749 xDelete<IssmDouble>(basis);
750 }
751 /*}}}*/
752+/*FUNCTION TriaRef::GetBFS {{{*/
753+void TriaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
754+ /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
755+ * For node i, Bvi can be expressed in the actual coordinate system
756+ * by: Bvi=[ dphi/dx 0 ]
757+ * [ 0 dphi/dy ]
758+ * [ 1/2*dphi/dy 1/2*dphi/dx]
759+ * [ 0 0 ]
760+ * [ dphi/dx dphi/dy ]
761+ *
762+ * by: Bpi=[ 0 ]
763+ * [ 0 ]
764+ * [ 0 ]
765+ * [ phi_p ]
766+ * [ 0 ]
767+ * where phi is the interpolation function for node i.
768+ * Same thing for Bb except the last column that does not exist.
769+ */
770+
771+ /*Fetch number of nodes for this finite element*/
772+ int pnumnodes = this->NumberofNodesPressure();
773+ int vnumnodes = this->NumberofNodesVelocity();
774+
775+ /*Get nodal functions derivatives*/
776+ IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);
777+ IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
778+ GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
779+ GetNodalFunctionsPressure(pbasis,gauss);
780+
781+ /*Build B: */
782+ for(int i=0;i<vnumnodes;i++){
783+ B[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
784+ B[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
785+ B[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
786+ B[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
787+ B[(2*vnumnodes+pnumnodes)*2+2*i+0] = .5*vdbasis[1*vnumnodes+i];
788+ B[(2*vnumnodes+pnumnodes)*2+2*i+1] = .5*vdbasis[0*vnumnodes+i];
789+ B[(2*vnumnodes+pnumnodes)*3+2*i+0] = 0.;
790+ B[(2*vnumnodes+pnumnodes)*3+2*i+1] = 0.;
791+ B[(2*vnumnodes+pnumnodes)*4+2*i+0] = vdbasis[0*vnumnodes+i];
792+ B[(2*vnumnodes+pnumnodes)*4+2*i+1] = vdbasis[1*vnumnodes+i];
793+ }
794+ for(int i=0;i<pnumnodes;i++){
795+ B[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
796+ B[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
797+ B[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
798+ B[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = pbasis[i];
799+ B[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = 0.;
800+ }
801+
802+ /*Clean up*/
803+ xDelete<IssmDouble>(vdbasis);
804+ xDelete<IssmDouble>(pbasis);
805+}
806+/*}}}*/
807+/*FUNCTION TriaRef::GetBprimeFS {{{*/
808+void TriaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss){
809+ /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
810+ * For node i, Bi' can be expressed in the actual coordinate system
811+ * by:
812+ * Bvi' = [ dphi/dx 0 ]
813+ * [ 0 dphi/dy ]
814+ * [ dphi/dy dphi/dx ]
815+ * [ dphi/dx dphi/dy ]
816+ * [ 0 0 ]
817+ *
818+ * by: Bpi=[ 0 ]
819+ * [ 0 ]
820+ * [ 0 ]
821+ * [ 0 ]
822+ * [ phi ]
823+ * where phi is the interpolation function for node i.
824+ */
825+
826+ /*Fetch number of nodes for this finite element*/
827+ int pnumnodes = this->NumberofNodesPressure();
828+ int vnumnodes = this->NumberofNodesVelocity();
829+
830+ /*Get nodal functions derivatives*/
831+ IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);
832+ IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
833+ GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
834+ GetNodalFunctionsPressure(pbasis,gauss);
835+
836+ /*Build B_prime: */
837+ for(int i=0;i<vnumnodes;i++){
838+ B_prime[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
839+ B_prime[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
840+ B_prime[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
841+ B_prime[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
842+ B_prime[(2*vnumnodes+pnumnodes)*2+2*i+0] = vdbasis[1*vnumnodes+i];
843+ B_prime[(2*vnumnodes+pnumnodes)*2+2*i+1] = vdbasis[0*vnumnodes+i];
844+ B_prime[(2*vnumnodes+pnumnodes)*3+2*i+0] = vdbasis[0*vnumnodes+i];
845+ B_prime[(2*vnumnodes+pnumnodes)*3+2*i+1] = vdbasis[1*vnumnodes+i];
846+ B_prime[(2*vnumnodes+pnumnodes)*4+2*i+0] = 0.;
847+ B_prime[(2*vnumnodes+pnumnodes)*4+2*i+1] = 0.;
848+ }
849+ for(int i=0;i<pnumnodes;i++){
850+ B_prime[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
851+ B_prime[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
852+ B_prime[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
853+ B_prime[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = 0.;
854+ B_prime[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = pbasis[i];
855+ }
856+
857+ /*Clean up*/
858+ xDelete<IssmDouble>(vdbasis);
859+ xDelete<IssmDouble>(pbasis);
860+}
861+/*}}}*/
862 /*FUNCTION TriaRef::GetBMasstransport{{{*/
863 void TriaRef::GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
864 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
865@@ -457,6 +567,76 @@
866 }
867 }
868 /*}}}*/
869+/*FUNCTION TriaRef::GetNodalFunctionsVelocity{{{*/
870+void TriaRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussTria* gauss){
871+ /*This routine returns the values of the nodal functions at the gaussian point.*/
872+
873+ switch(this->element_type){
874+ case P1P1Enum:
875+ this->element_type = P1Enum;
876+ this->GetNodalFunctions(basis,gauss);
877+ this->element_type = P1P1Enum;
878+ return;
879+ case P1P1GLSEnum:
880+ this->element_type = P1Enum;
881+ this->GetNodalFunctions(basis,gauss);
882+ this->element_type = P1P1GLSEnum;
883+ return;
884+ case MINIcondensedEnum:
885+ this->element_type = P1bubbleEnum;
886+ this->GetNodalFunctions(basis,gauss);
887+ this->element_type = MINIcondensedEnum;
888+ return;
889+ case MINIEnum:
890+ this->element_type = P1bubbleEnum;
891+ this->GetNodalFunctions(basis,gauss);
892+ this->element_type = MINIEnum;
893+ return;
894+ case TaylorHoodEnum:
895+ this->element_type = P2Enum;
896+ this->GetNodalFunctions(basis,gauss);
897+ this->element_type = TaylorHoodEnum;
898+ return;
899+ default:
900+ _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
901+ }
902+}
903+/*}}}*/
904+/*FUNCTION TriaRef::GetNodalFunctionsPressure{{{*/
905+void TriaRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussTria* gauss){
906+ /*This routine returns the values of the nodal functions at the gaussian point.*/
907+
908+ switch(this->element_type){
909+ case P1P1Enum:
910+ this->element_type = P1Enum;
911+ this->GetNodalFunctions(basis,gauss);
912+ this->element_type = P1P1Enum;
913+ return;
914+ case P1P1GLSEnum:
915+ this->element_type = P1Enum;
916+ this->GetNodalFunctions(basis,gauss);
917+ this->element_type = P1P1GLSEnum;
918+ return;
919+ case MINIcondensedEnum:
920+ this->element_type = P1Enum;
921+ this->GetNodalFunctions(basis,gauss);
922+ this->element_type = MINIcondensedEnum;
923+ return;
924+ case MINIEnum:
925+ this->element_type = P1Enum;
926+ this->GetNodalFunctions(basis,gauss);
927+ this->element_type = MINIEnum;
928+ return;
929+ case TaylorHoodEnum:
930+ this->element_type = P1Enum;
931+ this->GetNodalFunctions(basis,gauss);
932+ this->element_type = TaylorHoodEnum;
933+ return;
934+ default:
935+ _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
936+ }
937+}
938+/*}}}*/
939 /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/
940 void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){
941 /*This routine returns the values of the nodal functions at the gaussian point.*/
942@@ -528,6 +708,72 @@
943
944 }
945 /*}}}*/
946+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesPressure{{{*/
947+void TriaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){
948+ switch(this->element_type){
949+ case P1P1Enum:
950+ this->element_type = P1Enum;
951+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
952+ this->element_type = P1P1Enum;
953+ return;
954+ case P1P1GLSEnum:
955+ this->element_type = P1Enum;
956+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
957+ this->element_type = P1P1GLSEnum;
958+ return;
959+ case MINIcondensedEnum:
960+ this->element_type = P1Enum;
961+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
962+ this->element_type = MINIcondensedEnum;
963+ return;
964+ case MINIEnum:
965+ this->element_type = P1Enum;
966+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
967+ this->element_type = MINIEnum;
968+ return;
969+ case TaylorHoodEnum:
970+ this->element_type = P1Enum;
971+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
972+ this->element_type = TaylorHoodEnum;
973+ return;
974+ default:
975+ _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
976+ }
977+}
978+/*}}}*/
979+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesVelocity{{{*/
980+void TriaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){
981+ switch(this->element_type){
982+ case P1P1Enum:
983+ this->element_type = P1Enum;
984+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
985+ this->element_type = P1P1Enum;
986+ return;
987+ case P1P1GLSEnum:
988+ this->element_type = P1Enum;
989+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
990+ this->element_type = P1P1GLSEnum;
991+ return;
992+ case MINIcondensedEnum:
993+ this->element_type = P1bubbleEnum;
994+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
995+ this->element_type = MINIcondensedEnum;
996+ return;
997+ case MINIEnum:
998+ this->element_type = P1bubbleEnum;
999+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
1000+ this->element_type = MINIEnum;
1001+ return;
1002+ case TaylorHoodEnum:
1003+ this->element_type = P2Enum;
1004+ this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
1005+ this->element_type = TaylorHoodEnum;
1006+ return;
1007+ default:
1008+ _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
1009+ }
1010+}
1011+/*}}}*/
1012 /*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/
1013 void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss){
1014 /*This routine returns the values of the nodal functions derivatives (with respect to the
Note: See TracBrowser for help on using the repository browser.