source: issm/oecreview/Archive/22819-23185/ISSM-22873-22874.diff

Last change on this file was 23186, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22819-23185
File size: 33.7 KB
RevLine 
[23186]1Index: ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22873)
4+++ ../trunk-jpl/src/c/analyses/ThermalAnalysis.cpp (revision 22874)
5@@ -11,7 +11,7 @@
6 int M,N;
7 int finiteelement;
8 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
9- _assert_(finiteelement==P1Enum);
10+ _assert_(finiteelement==P1Enum);
11
12 /*Output*/
13 IssmDouble *spcvector = NULL;
14@@ -93,8 +93,8 @@
15
16 int finiteelement;
17 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
18- _assert_(finiteelement==P1Enum);
19-
20+ _assert_(finiteelement==P1Enum);
21+
22 if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
23 ::CreateNodes(nodes,iomodel,ThermalAnalysisEnum,finiteelement);
24 iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
25@@ -112,7 +112,7 @@
26 /*Update elements: */
27 int finiteelement;
28 iomodel->FindConstant(&finiteelement,"md.thermal.fe");
29- _assert_(finiteelement==P1Enum);
30+ _assert_(finiteelement==P1Enum);
31 int counter=0;
32 for(int i=0;i<iomodel->numberofelements;i++){
33 if(iomodel->my_elements[i]){
34@@ -188,8 +188,10 @@
35 iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
36 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
37 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
38- if (FrictionCoupling==1){
39- iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
40+ if (FrictionCoupling==3){
41+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
42+ else if(FrictionCoupling==4){
43+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
44 }
45 break;
46 case 2:
47@@ -201,8 +203,10 @@
48 iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
49 iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
50 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
51- if (FrictionCoupling==1){
52- iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
53+ if (FrictionCoupling==3){
54+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
55+ else if(FrictionCoupling==4){
56+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
57 }
58 break;
59 case 4:
60@@ -230,8 +234,10 @@
61 iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
62 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
63 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
64- if (FrictionCoupling==1){
65- iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
66+ if (FrictionCoupling==3){
67+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
68+ else if(FrictionCoupling==4){
69+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
70 }
71 break;
72 case 9:
73@@ -452,7 +458,7 @@
74 K[1][0]=h/(2.*vel)*fabs(vy*vx); K[1][1]=h/(2.*vel)*fabs(vy*vy); K[1][2]=h/(2.*vel)*fabs(vy*vz);
75 K[2][0]=h/(2.*vel)*fabs(vz*vx); K[2][1]=h/(2.*vel)*fabs(vz*vy); K[2][2]=h/(2.*vel)*fabs(vz*vz);
76 for(int i=0;i<3;i++) for(int j=0;j<3;j++) K[i][j] = D_scalar*K[i][j];
77- GetBAdvecprime(Bprime,element,xyz_list,gauss);
78+ GetBAdvecprime(Bprime,element,xyz_list,gauss);
79
80 TripleMultiply(Bprime,3,numnodes,1,
81 &K[0][0],3,3,0,
82@@ -487,7 +493,7 @@
83 return Ke;
84 }/*}}}*/
85 ElementVector* ThermalAnalysis::CreatePVector(Element* element){/*{{{*/
86-
87+
88 /* Check if ice in element */
89 if(!element->IsIceInElement()) return NULL;
90
91@@ -700,9 +706,9 @@
92
93 }/*}}}*/
94 void ThermalAnalysis::GetBAdvec(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
95- /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
96+ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
97 * For node i, Bi' can be expressed in the actual coordinate system
98- * by:
99+ * by:
100 * Bi_advec =[ h ]
101 * [ h ]
102 * [ h ]
103@@ -729,9 +735,9 @@
104 xDelete<IssmDouble>(basis);
105 }/*}}}*/
106 void ThermalAnalysis::GetBAdvecprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
107- /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
108+ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
109 * For node i, Bi' can be expressed in the actual coordinate system
110- * by:
111+ * by:
112 * Biprime_advec=[ dh/dx ]
113 * [ dh/dy ]
114 * [ dh/dz ]
115@@ -758,9 +764,9 @@
116 xDelete<IssmDouble>(dbasis);
117 }/*}}}*/
118 void ThermalAnalysis::GetBConduct(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
119- /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
120+ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1.
121 * For node i, Bi' can be expressed in the actual coordinate system
122- * by:
123+ * by:
124 * Bi_conduct=[ dh/dx ]
125 * [ dh/dy ]
126 * [ dh/dz ]
127Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
128===================================================================
129--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 22873)
130+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 22874)
131@@ -610,7 +610,7 @@
132 }
133 break;
134 case L1L2ApproximationEnum: numdofs =2; break;
135- case HOApproximationEnum:
136+ case HOApproximationEnum:
137 switch(domaintype){
138 case Domain3DEnum: numdofs=2; break;
139 case Domain2DverticalEnum: numdofs=1; break;
140@@ -784,8 +784,10 @@
141 iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
142 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
143 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
144- if(FrictionCoupling==1){
145- iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
146+ if(FrictionCoupling==3){
147+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
148+ else if(FrictionCoupling==4){
149+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
150 }
151 break;
152 case 2:
153@@ -797,8 +799,10 @@
154 iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
155 iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
156 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
157- if(FrictionCoupling==1){
158- iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
159+ if(FrictionCoupling==3){
160+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
161+ else if(FrictionCoupling==4){
162+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
163 }
164 break;
165 case 4:
166@@ -826,8 +830,11 @@
167 iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
168 iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
169 iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
170- if(FrictionCoupling==1){
171- iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
172+ if(FrictionCoupling==3){
173+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
174+ else if(FrictionCoupling==4){
175+ iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",EffectivePressureEnum);
176+
177 }
178 break;
179 case 9:
180@@ -933,9 +940,9 @@
181 else if(newton>0)
182 solutionsequence_newton(femmodel);
183 else
184- solutionsequence_nonlinear(femmodel,conserve_loads);
185+ solutionsequence_nonlinear(femmodel,conserve_loads);
186 }
187- else if(!isFS && (isSSA || isHO || isL1L2)){
188+ else if(!isFS && (isSSA || isHO || isL1L2)){
189 if(VerboseSolution()) _printf0_(" computing velocities\n");
190
191 femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
192@@ -942,7 +949,7 @@
193 if(newton>0)
194 solutionsequence_newton(femmodel);
195 else
196- solutionsequence_nonlinear(femmodel,conserve_loads);
197+ solutionsequence_nonlinear(femmodel,conserve_loads);
198
199 if(domaintype==Domain2DverticalEnum && isSSA){
200 femmodel->parameters->SetParam(VxEnum,InputToExtrudeEnum);
201@@ -965,7 +972,7 @@
202 int approximation;
203 element->GetInputValue(&approximation,ApproximationEnum);
204 switch(approximation){
205- case FSApproximationEnum:
206+ case FSApproximationEnum:
207 return CreateDVectorFS(element);
208 default:
209 return NULL; //no need for doftypes outside of FS approximation
210@@ -978,11 +985,11 @@
211 int approximation;
212 element->GetInputValue(&approximation,ApproximationEnum);
213 switch(approximation){
214- case SSAApproximationEnum:
215+ case SSAApproximationEnum:
216 return CreateJacobianMatrixSSA(element);
217- case HOApproximationEnum:
218+ case HOApproximationEnum:
219 return CreateJacobianMatrixHO(element);
220- case FSApproximationEnum:
221+ case FSApproximationEnum:
222 return CreateJacobianMatrixFS(element);
223 case NoneApproximationEnum:
224 return NULL;
225@@ -996,19 +1003,19 @@
226 switch(approximation){
227 case SIAApproximationEnum:
228 return NULL;
229- case SSAApproximationEnum:
230+ case SSAApproximationEnum:
231 return CreateKMatrixSSA(element);
232- case L1L2ApproximationEnum:
233+ case L1L2ApproximationEnum:
234 return CreateKMatrixL1L2(element);
235- case HOApproximationEnum:
236+ case HOApproximationEnum:
237 return CreateKMatrixHO(element);
238- case FSApproximationEnum:
239+ case FSApproximationEnum:
240 return CreateKMatrixFS(element);
241- case SSAHOApproximationEnum:
242+ case SSAHOApproximationEnum:
243 return CreateKMatrixSSAHO(element);
244- case HOFSApproximationEnum:
245+ case HOFSApproximationEnum:
246 return CreateKMatrixHOFS(element);
247- case SSAFSApproximationEnum:
248+ case SSAFSApproximationEnum:
249 return CreateKMatrixSSAFS(element);
250 case NoneApproximationEnum:
251 return NULL;
252@@ -1023,19 +1030,19 @@
253 switch(approximation){
254 case SIAApproximationEnum:
255 return NULL;
256- case SSAApproximationEnum:
257+ case SSAApproximationEnum:
258 return CreatePVectorSSA(element);
259- case L1L2ApproximationEnum:
260+ case L1L2ApproximationEnum:
261 return CreatePVectorL1L2(element);
262- case HOApproximationEnum:
263+ case HOApproximationEnum:
264 return CreatePVectorHO(element);
265- case FSApproximationEnum:
266+ case FSApproximationEnum:
267 return CreatePVectorFS(element);
268- case SSAHOApproximationEnum:
269+ case SSAHOApproximationEnum:
270 return CreatePVectorSSAHO(element);
271- case HOFSApproximationEnum:
272+ case HOFSApproximationEnum:
273 return CreatePVectorHOFS(element);
274- case SSAFSApproximationEnum:
275+ case SSAFSApproximationEnum:
276 return CreatePVectorSSAFS(element);
277 case NoneApproximationEnum:
278 return NULL;
279@@ -1125,12 +1132,12 @@
280 case FSApproximationEnum: case NoneApproximationEnum:
281 InputUpdateFromSolutionFS(solution,element);
282 return;
283- case SIAApproximationEnum:
284+ case SIAApproximationEnum:
285 return;
286- case SSAApproximationEnum:
287+ case SSAApproximationEnum:
288 InputUpdateFromSolutionSSA(solution,element);
289 return;
290- case HOApproximationEnum:
291+ case HOApproximationEnum:
292 InputUpdateFromSolutionHO(solution,element);
293 return;
294 case L1L2ApproximationEnum:
295@@ -1580,7 +1587,7 @@
296
297 /*Retrieve all inputs and parameters*/
298 element->GetVerticesCoordinates(&xyz_list);
299- Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
300+ Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
301 Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input);
302 IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
303
304@@ -1688,13 +1695,13 @@
305 return pe;
306 }/*}}}*/
307 void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
308- /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
309+ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
310 * For node i, Bi can be expressed in the actual coordinate system
311- * by:
312+ * by:
313 * 2D 1D
314 * Bi=[ dN/dx 0 ] Bi=[ dN/dx ]
315- * [ 0 dN/dy ]
316- * [ 1/2*dN/dy 1/2*dN/dx ]
317+ * [ 0 dN/dy ]
318+ * [ 1/2*dN/dy 1/2*dN/dx ]
319 * where N is the finiteelement function for node i.
320 *
321 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
322@@ -1728,9 +1735,9 @@
323 xDelete<IssmDouble>(dbasis);
324 }/*}}}*/
325 void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
326- /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2.
327+ /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2.
328 * For node i, Bi can be expressed in the actual coordinate system
329- * by:
330+ * by:
331 * 2D 1D
332 * Bi=[ N 0 ] Bi = N
333 * [ 0 N ]
334@@ -1765,9 +1772,9 @@
335 xDelete<IssmDouble>(basis);
336 }/*}}}*/
337 void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
338- /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
339+ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
340 * For node i, Bi' can be expressed in the actual coordinate system
341- * by:
342+ * by:
343 * 2D 1D
344 * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ]
345 * [ dN/dx 2*dN/dy ]
346@@ -2049,7 +2056,7 @@
347
348 /*Retrieve all inputs and parameters*/
349 element->GetVerticesCoordinates(&xyz_list);
350- Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
351+ Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
352 Input* surface_input =element->GetInput(SurfaceEnum); _assert_(surface_input);
353 IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
354
355@@ -2676,16 +2683,16 @@
356 return pe;
357 }/*}}}*/
358 void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
359- /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
360+ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
361 * For node i, Bi can be expressed in the actual coordinate system
362- * by:
363+ * by:
364 * 3D 2D
365 *
366 * Bi=[ dh/dx 0 ] Bi=[ dh/dx]
367 * [ 0 dh/dy ] [ dh/dy]
368- * [ 1/2*dh/dy 1/2*dh/dx ]
369- * [ 1/2*dh/dz 0 ]
370- * [ 0 1/2*dh/dz ]
371+ * [ 1/2*dh/dy 1/2*dh/dx ]
372+ * [ 1/2*dh/dz 0 ]
373+ * [ 0 1/2*dh/dz ]
374 * where h is the interpolation function for node i.
375 *
376 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
377@@ -2725,9 +2732,9 @@
378 xDelete<IssmDouble>(dbasis);
379 }/*}}}*/
380 void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
381- /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2.
382+ /*Compute B matrix. B=[B1 B2 B3] where Bi is square and of size 2.
383 * For node i, Bi can be expressed in the actual coordinate system
384- * by:
385+ * by:
386 * 3D 2D
387 * Bi=[ N 0 ] Bi=N
388 * [ 0 N ]
389@@ -2762,13 +2769,13 @@
390 xDelete<IssmDouble>(basis);
391 }/*}}}*/
392 void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
393- /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
394+ /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
395 * For node i, Bi' can be expressed in the actual coordinate system
396- * by:
397+ * by:
398 * 3D 2D
399 * Bi_prime=[ 2*dN/dx dN/dy ] Bi_prime=[ 2*dN/dx ]
400 * [ dN/dx 2*dN/dy ] [ dN/dy ]
401- * [ dN/dy dN/dx ]
402+ * [ dN/dy dN/dx ]
403 * where hNis the finiteelement function for node i.
404 *
405 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
406@@ -3083,7 +3090,7 @@
407 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
408 if(dim==2) slope2=slope[0]*slope[0];
409 else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
410- scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt;
411+ scalar = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt;
412 for(i=0;i<vnumnodes;i++){
413 for(j=0;j<vnumnodes;j++){
414 Ke->values[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j];
415@@ -3285,7 +3292,7 @@
416 BtBUzawa,pnumdof,numdof,0,
417 &Ke->values[0],1);
418
419- if(element->IsOnBase() && 0){
420+ if(element->IsOnBase() && 0){
421 element->FindParam(&rl,AugmentedLagrangianRlambdaEnum);
422 element->GetVerticesCoordinatesBase(&xyz_list_base);
423 element->NormalBase(&normal[0],xyz_list_base);
424@@ -4084,7 +4091,7 @@
425 for(int ig=gauss->begin();ig<gauss->end();ig++){
426 gauss->GaussPoint(ig);
427 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
428-
429+
430 pressure_input->GetInputValue(&pressure, gauss);
431 element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss);
432
433@@ -4095,7 +4102,7 @@
434 }
435 }
436
437- if(element->IsOnBase() && 0){
438+ if(element->IsOnBase() && 0){
439 IssmDouble sigmann;
440 IssmDouble* vbasis = xNew<IssmDouble>(numnodes);
441
442@@ -4304,7 +4311,7 @@
443 return pe;
444 }/*}}}*/
445 void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
446- /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
447+ /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
448 * For node i, Bvi can be expressed in the actual coordinate system
449 * by: Bvi=[ dphi/dx 0 ]
450 * [ 0 dphi/dy ]
451@@ -4416,9 +4423,9 @@
452 xDelete<IssmDouble>(pbasis);
453 }/*}}}*/
454 void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
455- /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
456+ /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
457 * For node i, Li can be expressed in the actual coordinate system
458- * by in 3d
459+ * by in 3d
460 * Li=[ h 0 0 0 ]
461 * [ 0 h 0 0 ]
462 * in 2d:
463@@ -4476,9 +4483,9 @@
464 xDelete<IssmDouble>(vbasis);
465 }/*}}}*/
466 void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
467- /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
468+ /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
469 * For node i, Bi' can be expressed in the actual coordinate system
470- * by:
471+ * by:
472 * Bvi' = [ dphi/dx 0 ]
473 * [ 0 dphi/dy ]
474 * [ dphi/dy dphi/dx ]
475@@ -4589,9 +4596,9 @@
476 xDelete<IssmDouble>(pbasis);
477 }/*}}}*/
478 void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
479- /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2.
480+ /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2.
481 * For node i, Bi' can be expressed in the actual coordinate system
482- * by:
483+ * by:
484 * Bvi' = [ dphi/dx dphi/dy ]
485 *
486 * In 3d
487@@ -4625,9 +4632,9 @@
488 xDelete<IssmDouble>(vdbasis);
489 }/*}}}*/
490 void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
491- /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
492+ /* Compute B' matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
493 * For node i, Bi' can be expressed in the actual coordinate system
494- * by:
495+ * by:
496 * Bvi' = [ dphi/dx 0 ]
497 * [ 0 dphi/dy ]
498 * [ dphi/dy dphi/dx ]
499@@ -4688,7 +4695,7 @@
500 xDelete<IssmDouble>(vdbasis);
501 }/*}}}*/
502 void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
503- /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi.
504+ /*Compute B matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi.
505 */
506
507 /*Fetch number of nodes for this finite element*/
508@@ -4710,7 +4717,7 @@
509 xDelete<IssmDouble>(basis);
510 }/*}}}*/
511 void StressbalanceAnalysis::GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
512- /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
513+ /*Compute B matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3.
514 * For node i, Bvi can be expressed in the actual coordinate system
515 * by: Bvi=[ dphi/dx 0 ]
516 * [ 0 dphi/dy ]
517@@ -4775,7 +4782,7 @@
518 }/*}}}*/
519 void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
520 /*Compute C matrix. C=[Cp1 Cp2 ...] where:
521- * Cpi=[phi phi].
522+ * Cpi=[phi phi].
523 */
524
525 /*Fetch number of nodes for this finite element*/
526@@ -4795,7 +4802,7 @@
527 xDelete<IssmDouble>(basis);
528 }/*}}}*/
529 void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
530- /* Compute C' matrix. C'=[C1' C2' ...]
531+ /* Compute C' matrix. C'=[C1' C2' ...]
532 * Ci' = [ phi 0 ]
533 * [ 0 phi ]
534 *
535@@ -4929,15 +4936,15 @@
536 /*Allocate new inputs*/
537 int tnumnodes = element->GetNumberOfVertices(); //Tensors, P1 DG
538 IssmDouble* epsxx = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxx = xNew<IssmDouble>(tnumnodes);
539- IssmDouble* epsyy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapyy = xNew<IssmDouble>(tnumnodes);
540- IssmDouble* epsxy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxy = xNew<IssmDouble>(tnumnodes);
541- IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL;
542- IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL;
543- IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL;
544+ IssmDouble* epsyy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapyy = xNew<IssmDouble>(tnumnodes);
545+ IssmDouble* epsxy = xNew<IssmDouble>(tnumnodes); IssmDouble* sigmapxy = xNew<IssmDouble>(tnumnodes);
546+ IssmDouble* epszz = NULL; IssmDouble* sigmapzz = NULL;
547+ IssmDouble* epsxz = NULL; IssmDouble* sigmapxz = NULL;
548+ IssmDouble* epsyz = NULL; IssmDouble* sigmapyz = NULL;
549 if(dim==3){
550 epszz = xNew<IssmDouble>(tnumnodes); sigmapzz = xNew<IssmDouble>(tnumnodes);
551- epsxz = xNew<IssmDouble>(tnumnodes); sigmapxz = xNew<IssmDouble>(tnumnodes);
552- epsyz = xNew<IssmDouble>(tnumnodes); sigmapyz = xNew<IssmDouble>(tnumnodes);
553+ epsxz = xNew<IssmDouble>(tnumnodes); sigmapxz = xNew<IssmDouble>(tnumnodes);
554+ epsyz = xNew<IssmDouble>(tnumnodes); sigmapyz = xNew<IssmDouble>(tnumnodes);
555 }
556
557 /*Get d and tau*/
558@@ -5171,9 +5178,9 @@
559 }
560 epsxx = dvx[0];
561 epsyy = dvy[1];
562- epsxy = 0.5*(dvx[1] + dvy[0]);
563+ epsxy = 0.5*(dvx[1] + dvy[0]);
564 if(dim==3){
565- epszz = dvz[2];
566+ epszz = dvz[2];
567 epsxz = 0.5*(dvx[2] + dvz[0]);
568 epsyz = 0.5*(dvy[2] + dvz[1]);
569 }
570@@ -5182,7 +5189,7 @@
571 IssmDouble coef1,coef2,coef3;
572 B_input->GetInputValue(&B,gauss);
573 n_input->GetInputValue(&n,gauss);
574- coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2) )^(n-1)/n )
575+ coef1 = B*pow(1./sqrt(2.),(1.-n)/n); //2 eta_0 = 2 * B/(2* (1/sqrt(2) )^(n-1)/n )
576 coef2 = r;
577 if(dim==2){
578 coef3 = sqrt(
579@@ -5206,7 +5213,7 @@
580 dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + 2.*epsxy_old*epsxy_old );
581 }
582 else{
583- dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old
584+ dnorm = sqrt( epsxx_old*epsxx_old + epsyy_old*epsyy_old + epszz_old*epszz_old
585 +2.*(epsxy_old*epsxy_old + epsxz_old*epsxz_old + epsyz_old*epsyz_old));
586 }
587 /*Initial guess cannot be 0 otherwise log(0) - inf*/
588@@ -5359,9 +5366,9 @@
589 }
590 epsxx = dvx[0];
591 epsyy = dvy[1];
592- epsxy = 0.5*(dvx[1] + dvy[0]);
593+ epsxy = 0.5*(dvx[1] + dvy[0]);
594 if(dim==3){
595- epszz = dvz[2];
596+ epszz = dvz[2];
597 epsxz = 0.5*(dvx[2] + dvz[0]);
598 epsyz = 0.5*(dvy[2] + dvz[1]);
599 }
600@@ -5613,7 +5620,7 @@
601
602 for(i=0;i<numdof2dm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftot+j+numdofm]+=Ke_drag[i][j];
603 for(i=0;i<numdof2d;i++) for(j=0;j<numdof2dm;j++) Ke->values[(i+numdofm)*numdoftot+j]+=Ke_drag2[i][j];
604-
605+
606 /*Transform Coordinate System*/
607 element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
608
609@@ -5646,8 +5653,8 @@
610 IssmDouble Bprime2[3][numdofs];
611 IssmDouble D[4][4]={0.0}; // material matrix, simple scalar matrix.
612 IssmDouble D2[3][3]={0.0}; // material matrix, simple scalar matrix.
613- IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix
614- IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix
615+ IssmDouble Ke_gg[numdofs][numdofm]={0.0}; //local element stiffness matrix
616+ IssmDouble Ke_gg2[numdofm][numdofs]={0.0}; //local element stiffness matrix
617 IssmDouble *xyz_list = NULL;
618
619 /*Find penta on bed as FS must be coupled to the dofs on the bed: */
620@@ -5718,7 +5725,7 @@
621 &Bprime2[0][0],3,numdofs,0,
622 &Ke_gg2[0][0],1);
623
624- }
625+ }
626 for(i=0;i<numdofs;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i][j];
627 for(i=0;i<numdofm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg2[i][j];
628
629@@ -5802,7 +5809,7 @@
630 this->GetBHOFriction(L,element,3,xyz_list_tria,gauss);
631
632 DL_scalar=alpha2*gauss->weight*Jdet2d;
633- for (i=0;i<2;i++) DL[i][i]=DL_scalar;
634+ for (i=0;i<2;i++) DL[i][i]=DL_scalar;
635
636 /* Do the triple producte tL*D*L: */
637 TripleMultiply( L,2,numdof,1,
638@@ -5883,7 +5890,7 @@
639
640 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
641 this->GetBSSAHO(B, element,xyz_list, gauss);
642- this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria);
643+ this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria);
644 element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
645
646 D_scalar=2*viscosity*gauss->weight*Jdet;
647@@ -5893,7 +5900,7 @@
648 &D[0][0],3,3,0,
649 Bprime,3,numdofm,0,
650 Ke_gg,1);
651- }
652+ }
653 for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
654 for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
655
656@@ -5902,7 +5909,7 @@
657
658 /*Clean-up and return*/
659 basaltria->DeleteMaterials(); delete basaltria;
660-
661+
662 delete gauss;
663 delete gauss_tria;
664 xDelete<IssmDouble>(B);
665@@ -5940,9 +5947,9 @@
666 int indices[3]={18,19,20};
667 Ke1->StaticCondensation(3,&indices[0]);
668 int init = element->FiniteElement();
669- element->SetTemporaryElementType(P1Enum);
670+ element->SetTemporaryElementType(P1Enum);
671 ElementMatrix* Ke2=CreateKMatrixSSA3d(element);
672- element->SetTemporaryElementType(init);
673+ element->SetTemporaryElementType(init);
674 ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element);
675 ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
676
677@@ -5986,7 +5993,7 @@
678 /*Initialize Element matrix and return if necessary*/
679 if(element->IsFloating() || !element->IsOnBase()) return NULL;
680
681- /*Build a tria element using the 3 nodes of the base of the penta. Then use
682+ /*Build a tria element using the 3 nodes of the base of the penta. Then use
683 * the tria functionality to build a friction stiffness matrix on these 3
684 * nodes: */
685 Element* basalelement = element->SpawnBasalElement();
686@@ -6222,7 +6229,7 @@
687 element->JacobianDeterminant(&Jdet, xyz_list,gauss);
688 element->NodalFunctionsP1(&basis[0],gauss);
689 element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
690-
691+
692 element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
693 vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
694
695@@ -6468,9 +6475,9 @@
696 return pe;
697 }/*}}}*/
698 void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
699- /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.
700+ /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2.
701 * For node i, Bprimei can be expressed in the actual coordinate system
702- * by:
703+ * by:
704 * Bprimei=[ 2*dh/dx dh/dy 0 0 ]
705 * [ dh/dx 2*dh/dy 0 0 ]
706 * [ dh/dy dh/dx 0 0 ]
707@@ -6517,9 +6524,9 @@
708 }
709 }/*}}}*/
710 void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
711- /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
712+ /*Compute Bprime matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2.
713 * For node i, Bprimei can be expressed in the actual coordinate system
714- * by:
715+ * by:
716 * Bprimei=[ dN/dx 0 ]
717 * [ 0 dN/dy ]
718 * [ dN/dy dN/dx ]
719@@ -6552,9 +6559,9 @@
720 xDelete<IssmDouble>(dbasis);
721 }/*}}}*/
722 void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
723- /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
724+ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2.
725 * For node i, Bi can be expressed in the actual coordinate system
726- * by:
727+ * by:
728 * Bi=[ dh/dx 0 0 0 ]
729 * [ 0 dh/dy 0 0 ]
730 * [ 1/2*dh/dy 1/2*dh/dx 0 0 ]
731@@ -6610,9 +6617,9 @@
732 }
733 }/*}}}*/
734 void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
735- /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
736+ /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
737 * For node i, Bi can be expressed in the actual coordinate system
738- * by:
739+ * by:
740 * Bi=[ dN/dx 0 ]
741 * [ 0 dN/dy ]
742 * [ 1/2*dN/dy 1/2*dN/dx ]
743@@ -6642,9 +6649,9 @@
744 xDelete<IssmDouble>(dbasis);
745 }/*}}}*/
746 void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
747- /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2.
748+ /*Compute B matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2.
749 * For node i, Bi can be expressed in the actual coordinate system
750- * by:
751+ * by:
752 * Bi=[ dh/dx 0 ]
753 * [ 0 dh/dy ]
754 * [ 1/2*dh/dy 1/2*dh/dx ]
755@@ -6673,9 +6680,9 @@
756 xDelete<IssmDouble>(dbasis);
757 }/*}}}*/
758 void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
759- /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
760+ /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
761 * For node i, Lpi can be expressed in the actual coordinate system
762- * by:
763+ * by:
764 * Lpi=[ h 0 ]
765 * [ 0 h ]
766 * [ h 0 ]
767@@ -6707,9 +6714,9 @@
768 }
769 }/*}}}*/
770 void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
771- /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
772+ /* Compute Lprime matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
773 * For node i, Lpi can be expressed in the actual coordinate system
774- * by:
775+ * by:
776 * Lpi=[ h 0 0 0]
777 * [ 0 h 0 0]
778 * [ 0 0 h 0]
779@@ -6812,9 +6819,9 @@
780 }
781 }/*}}}*/
782 void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
783- /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
784+ /* Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
785 * For node i, Li can be expressed in the actual coordinate system
786- * by:
787+ * by:
788 * Li=[ h 0 0 ]
789 * [ 0 h 0 ]
790 * [ 0 0 h ]
791@@ -6852,9 +6859,9 @@
792 }/*}}}*/
793 void StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
794 /*
795- * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
796+ * Compute L matrix. L=[L1 L2 L3] where Li is square and of size numdof.
797 * For node i, Li can be expressed in the actual coordinate system
798- * by:
799+ * by:
800 * Li=[ h 0 ]
801 * [ 0 h ]
802 * [ h 0 ]
803@@ -7157,7 +7164,7 @@
804 element->GetInputListOnNodes(&vz[0],VzEnum,0.);
805 for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
806
807- /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
808+ /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
809 *so the pressure is just the pressure at the bedrock: */
810 rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
811 g = element->GetMaterialParameter(ConstantsGEnum);
Note: See TracBrowser for help on using the repository browser.