source: issm/oecreview/Archive/17984-18295/ISSM-18283-18284.diff@ 18296

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

Added 17984-18295
File size: 13.0 KB
RevLine 
[18296]1Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
2===================================================================
3--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18283)
4+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 18284)
5@@ -78,6 +78,8 @@
6 ElementVector* CreatePVectorFSViscousXTH(Element* element);
7 ElementVector* CreatePVectorFSShelf(Element* element);
8 ElementVector* CreatePVectorFSFront(Element* element);
9+ ElementVector* CreatePVectorFSFriction(Element* element);
10+ ElementVector* CreatePVectorFSStress(Element* element);
11 void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
12 void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
13 void GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
14Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
15===================================================================
16--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18283)
17+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 18284)
18@@ -244,6 +244,7 @@
19 iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
20 if(fe_FS==LATaylorHoodEnum){
21 iomodel->FetchDataToInput(elements,PressureEnum);
22+ InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
23 }
24
25 /*Friction law variables*/
26@@ -2895,13 +2896,8 @@
27
28 /*Fetch number of nodes and dof for this finite element*/
29 int vnumnodes = element->NumberofNodesVelocity();
30- int pnumnodes;
31- if(dim==2) pnumnodes=3;
32- else pnumnodes=6;
33- int lnumnodes;
34- if(dim==2) lnumnodes=2;
35- else lnumnodes=3;
36- //int pnumnodes = element->NumberofNodes(P1Enum);
37+ int pnumnodes = element->GetNumberOfNodes(P1Enum);
38+ int lnumnodes = element->GetNumberOfNodes(P2Enum);
39 int numdof = vnumnodes*dim;
40 int pnumdof = pnumnodes;
41 int lnumdof = lnumnodes;
42@@ -2919,16 +2915,13 @@
43 IssmDouble* BU = xNew<IssmDouble>(pnumdof);
44 IssmDouble* BprimeU = xNew<IssmDouble>(numdof);
45 IssmDouble* D = xNewZeroInit<IssmDouble>(epssize*epssize);
46- IssmDouble* CtCUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
47- IssmDouble* C = xNew<IssmDouble>(lnumdof);
48- IssmDouble* Cprime = xNew<IssmDouble>(numdof);
49
50 /*Retrieve all inputs and parameters*/
51 element->GetVerticesCoordinates(&xyz_list);
52- Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
53- Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
54- Input* vz_input;
55- if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
56+ Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input);
57+ Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input);
58+ Input* vz_input = NULL;
59+ if(dim==3){vz_input = element->GetInput(VzEnum); _assert_(vz_input);}
60
61 /* Start looping on the number of gaussian points: */
62 Gauss* gauss = element->NewGauss(5);
63@@ -2956,7 +2949,6 @@
64 &DU,1,1,0,
65 BprimeU,1,numdof,0,
66 BtBUzawa,1);
67-
68 }
69
70 if(element->IsOnBase()){
71@@ -2964,30 +2956,37 @@
72 element->GetVerticesCoordinatesBase(&xyz_list_base);
73 element->NormalBase(&normal[0],xyz_list_base);
74
75- int lsize;
76- IssmDouble* Dlambda = xNewZeroInit<IssmDouble>(dim*dim);
77- IssmDouble* C = xNewZeroInit<IssmDouble>(dim*lnumdof);
78- IssmDouble* Cprime = xNewZeroInit<IssmDouble>(dim*numdof);
79+ IssmDouble* Dlambda = xNewZeroInit<IssmDouble>(dim*dim);
80+ IssmDouble* C = xNewZeroInit<IssmDouble>(dim*lnumdof);
81+ IssmDouble* Cprime = xNewZeroInit<IssmDouble>(dim*numdof);
82+ IssmDouble* CtCUzawa = xNewZeroInit<IssmDouble>(numdof*lnumdof);
83
84 delete gauss;
85- Gauss* gauss=element->NewGaussBase(5);
86+ gauss = element->NewGaussBase(5);
87 for(int ig=gauss->begin();ig<gauss->end();ig++){
88 gauss->GaussPoint(ig);
89
90 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
91 this->GetCFS(C,element,dim,xyz_list,gauss);
92 this->GetCFSprime(Cprime,element,dim,xyz_list,gauss);
93- for(i=0;i<dim;i++) Dlambda[i*epssize+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
94+ for(i=0;i<dim;i++) Dlambda[i*dim+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
95 TripleMultiply(C,dim,lnumdof,1,
96 Dlambda,dim,dim,0,
97 Cprime,dim,numdof,0,
98 CtCUzawa,1);
99 }
100+
101+ /*The sigma naugmentation should not be transformed*/
102+ MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
103+ CtCUzawa,lnumdof,numdof,0,
104+ &Ke->values[0],1);
105+
106 /*Delete base part*/
107- xDelete<IssmDouble>(xyz_list_base);
108 xDelete<IssmDouble>(Dlambda);
109 xDelete<IssmDouble>(C);
110 xDelete<IssmDouble>(Cprime);
111+ xDelete<IssmDouble>(CtCUzawa);
112+ xDelete<IssmDouble>(xyz_list_base);
113 }
114
115 /*Transform Coordinate System*/
116@@ -2998,11 +2997,6 @@
117 BtBUzawa,pnumdof,numdof,0,
118 &Ke->values[0],1);
119
120- /*The sigma naugmentation should not be transformed*/
121- MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
122- CtCUzawa,lnumdof,numdof,0,
123- &Ke->values[0],1);
124-
125 /*Clean up and return*/
126 delete gauss;
127 xDelete<IssmDouble>(xyz_list);
128@@ -3012,9 +3006,6 @@
129 xDelete<IssmDouble>(BprimeU);
130 xDelete<IssmDouble>(BU);
131 xDelete<IssmDouble>(BtBUzawa);
132- xDelete<IssmDouble>(Cprime);
133- xDelete<IssmDouble>(C);
134- xDelete<IssmDouble>(CtCUzawa);
135 xDelete<int>(cs_list);
136 return Ke;
137 }/*}}}*/
138@@ -3297,6 +3288,21 @@
139 }/*}}}*/
140 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
141
142+ ElementVector* pe = NULL;
143+
144+ ElementVector* pe1=CreatePVectorFSViscous(element);
145+ ElementVector* pe2=CreatePVectorFSFriction(element);
146+ ElementVector* pe3=CreatePVectorFSStress(element);
147+ pe =new ElementVector(pe1,pe2,pe3);
148+ delete pe1;
149+ delete pe2;
150+ delete pe3;
151+
152+ /*clean-up and return*/
153+ return pe;
154+}/*}}}*/
155+ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
156+
157 int i,dim,fe_FS;
158 IssmDouble x_coord,y_coord,z_coord;
159 IssmDouble Jdet,forcex,forcey,forcez;
160@@ -3370,6 +3376,116 @@
161 }
162 return pe;
163 }/*}}}*/
164+ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/
165+
166+ if(!element->IsOnBase()) return NULL;
167+
168+ /*Intermediaries*/
169+ int dim;
170+ IssmDouble alpha2,Jdet;
171+ IssmDouble bed_normal[3];
172+ IssmDouble *xyz_list_base = NULL;
173+ Gauss* gauss = NULL;
174+
175+ /*Get problem dimension*/
176+ element->FindParam(&dim,DomainDimensionEnum);
177+
178+ /*Fetch number of nodes and dof for this finite element*/
179+ int vnumnodes = element->NumberofNodesVelocity();
180+
181+ /*Initialize Element matrix and vectors*/
182+ ElementVector* pe = element->NewElementVector(FSvelocityEnum);
183+ IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
184+
185+ /*Retrieve all inputs and parameters*/
186+ element->GetVerticesCoordinatesBase(&xyz_list_base);
187+ Input* alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input);
188+
189+ /* Start looping on the number of gaussian points: */
190+ gauss=element->NewGaussBase(3);
191+ for(int ig=gauss->begin();ig<gauss->end();ig++){
192+ gauss->GaussPoint(ig);
193+
194+ alpha2_input->GetInputValue(&alpha2, gauss);
195+ element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
196+ element->NodalFunctionsVelocity(vbasis,gauss);
197+ element->NormalBase(&bed_normal[0],xyz_list_base);
198+
199+ for(int i=0;i<vnumnodes;i++){
200+ pe->values[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1];
201+ pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0];
202+ if(dim==3){
203+ pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
204+ }
205+ }
206+
207+ }
208+
209+ /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
210+
211+ /*Clean up and return*/
212+ delete gauss;
213+ xDelete<IssmDouble>(xyz_list_base);
214+ xDelete<IssmDouble>(vbasis);
215+ return pe;
216+}/*}}}*/
217+ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/
218+
219+ if(!element->IsOnBase()) return NULL;
220+
221+ /*Intermediaries*/
222+ int dim;
223+ IssmDouble sigmann,sigmant,Jdet,bedslope,beta;
224+ IssmDouble *xyz_list_base = NULL;
225+ Gauss* gauss = NULL;
226+
227+ /*Get problem dimension*/
228+ element->FindParam(&dim,DomainDimensionEnum);
229+
230+ /*Fetch number of nodes and dof for this finite element*/
231+ int vnumnodes = element->NumberofNodesVelocity();
232+
233+ /*Initialize Element matrix and vectors*/
234+ ElementVector* pe = element->NewElementVector(FSvelocityEnum);
235+ IssmDouble* vbasis = xNew<IssmDouble>(vnumnodes);
236+
237+ /*Retrieve all inputs and parameters*/
238+ element->GetVerticesCoordinatesBase(&xyz_list_base);
239+ Input* sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input);
240+ Input* sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input);
241+ Input* bedslope_input=element->GetInput(BedSlopeXEnum); _assert_(bedslope_input);
242+
243+ /* Start looping on the number of gaussian points: */
244+ gauss=element->NewGaussBase(3);
245+ for(int ig=gauss->begin();ig<gauss->end();ig++){
246+ gauss->GaussPoint(ig);
247+
248+ sigmann_input->GetInputValue(&sigmann, gauss);
249+ sigmant_input->GetInputValue(&sigmant, gauss);
250+ bedslope_input->GetInputValue(&bedslope, gauss);
251+ element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
252+ element->NodalFunctionsVelocity(vbasis,gauss);
253+
254+ beta=sqrt(1+bedslope*bedslope);
255+ for(int i=0;i<vnumnodes;i++){
256+ pe->values[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i];
257+ pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i];
258+ if(dim==3){
259+ //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
260+ _error_("3d not supported yet");
261+ }
262+ }
263+
264+ }
265+
266+ /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
267+
268+ /*Clean up and return*/
269+ delete gauss;
270+ xDelete<IssmDouble>(xyz_list_base);
271+ xDelete<IssmDouble>(vbasis);
272+ return pe;
273+}/*}}}*/
274 #else
275 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
276
277@@ -3501,7 +3617,6 @@
278 /*clean-up and return*/
279 return pe;
280 }/*}}}*/
281-#endif
282 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
283
284 int i,dim;
285@@ -3571,6 +3686,7 @@
286 xDelete<IssmDouble>(xyz_list);
287 return pe;
288 }/*}}}*/
289+#endif
290 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/
291
292 int i,tausize,dim;
293@@ -3771,7 +3887,7 @@
294 element->FindParam(&r,AugmentedLagrangianREnum);
295 element->GetVerticesCoordinates(&xyz_list);
296
297- /*Get d and tau*/
298+ /*Get pressure and sigmann*/
299 Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
300 Input* sigmann_input =element->GetInput(SigmaNNEnum); _assert_(sigmann_input);
301
302@@ -3786,14 +3902,11 @@
303 for(i=0;i<numnodes;i++){
304 pe->values[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
305 pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
306- if(dim==3){
307- pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
308- }
309+ if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
310 }
311 }
312
313 if(element->IsOnBase()){
314-
315 IssmDouble sigmann;
316 IssmDouble* vbasis = xNew<IssmDouble>(numnodes);
317
318@@ -3801,7 +3914,7 @@
319 element->NormalBase(&bed_normal[0],xyz_list_base);
320
321 delete gauss;
322- Gauss* gauss=element->NewGaussBase(5);
323+ gauss=element->NewGaussBase(5);
324 for(int ig=gauss->begin();ig<gauss->end();ig++){
325 gauss->GaussPoint(ig);
326
327@@ -3810,11 +3923,9 @@
328 sigmann_input->GetInputValue(&sigmann, gauss);
329
330 for(i=0;i<numnodes;i++){
331- pe->values[i*dim+0] += - sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
332- pe->values[i*dim+1] += - sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
333- if(dim==3){
334- pe->values[i*dim+2]+= - sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
335- }
336+ pe->values[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
337+ pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
338+ if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
339 }
340 }
341 xDelete<IssmDouble>(xyz_list_base);
342@@ -3822,7 +3933,7 @@
343 }
344
345 /*Transform coordinate system*/
346- //element->TransformLoadVectorCoord(pe,cs_list); Do not transform pressure augmentation
347+ //element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation
348
349 /*Clean up and return*/
350 delete gauss;
351@@ -4470,19 +4581,16 @@
352 */
353
354 /*Fetch number of nodes for this finite element*/
355- int lnumnodes;
356- if(dim==2) lnumnodes=3;
357- else lnumnodes=6;
358- //int pnumnodes = element->NumberofNodes(P1Enum);
359+ int lnumnodes = element->GetNumberOfNodes(P2Enum);
360
361 /*Get nodal functions derivatives*/
362 IssmDouble* basis =xNew<IssmDouble>(lnumnodes);
363- element->NodalFunctions(basis,gauss);
364+ element->NodalFunctionsP2(basis,gauss);
365
366 /*Build B: */
367 for(int i=0;i<lnumnodes;i++){
368- C[i*lnumnodes+0] = basis[i];
369- C[i*lnumnodes+1] = basis[i];
370+ C[lnumnodes*0+i] = basis[i];
371+ C[lnumnodes*1+i] = basis[i];
372 }
373
374 /*Clean up*/
375@@ -4501,7 +4609,7 @@
376 */
377
378 /*Fetch number of nodes for this finite element*/
379- int vnumnodes = element->GetNumberOfNodes();
380+ int vnumnodes = element->NumberofNodesVelocity();
381 int vnumdof = vnumnodes*dim;
382
383 IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
Note: See TracBrowser for help on using the repository browser.