source: issm/oecreview/Archive/16554-17801/ISSM-17390-17391.diff@ 17802

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

Added archives
File size: 15.2 KB
RevLine 
[17802]1Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h
2===================================================================
3--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 17390)
4+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 17391)
5@@ -25,8 +25,8 @@
6 ElementMatrix* CreateJacobianMatrix(Element* element);
7 ElementMatrix* CreateKMatrix(Element* element);
8 ElementVector* CreatePVector(Element* element);
9- void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
10- void GetBprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
11+ void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
12+ void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
13 void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
14 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
15 void UpdateConstraints(FemModel* femmodel);
16Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
17===================================================================
18--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17390)
19+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17391)
20@@ -116,38 +116,41 @@
21 }/*}}}*/
22 ElementMatrix* DamageEvolutionAnalysis::CreateKMatrix(Element* element){/*{{{*/
23
24+ /* Check if ice in element */
25+ if(!element->IsIceInElement()) return NULL;
26+
27 /*Intermediaries*/
28- int meshtype;
29- Element* basalelement;
30+ Element* basalelement;
31+ int meshtype,dim;
32+ int stabilization;
33+ IssmDouble Jdet,dt,D_scalar,h;
34+ IssmDouble vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
35+ IssmDouble *xyz_list = NULL;
36
37- /*Get basal element*/
38+ /*Get problem dimension and basal element*/
39 element->FindParam(&meshtype,MeshTypeEnum);
40 switch(meshtype){
41 case Mesh2DhorizontalEnum:
42 basalelement = element;
43+ dim = 2;
44 break;
45 case Mesh3DEnum:
46 if(!element->IsOnBed()) return NULL;
47 basalelement = element->SpawnBasalElement();
48+ dim = 2;
49 break;
50 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
51 }
52
53- /*Intermediaries */
54- int stabilization;
55- IssmDouble Jdet,dt,D_scalar;
56- IssmDouble vel,vx,vy,dvxdx,dvydy,dvx[2],dvy[2];
57- IssmDouble *xyz_list = NULL;
58-
59 /*Fetch number of nodes and dof for this finite element*/
60 int numnodes = basalelement->GetNumberOfNodes();
61
62 /*Initialize Element vector*/
63- ElementMatrix* Ke = basalelement->NewElementMatrix(NoneApproximationEnum);
64+ ElementMatrix* Ke = basalelement->NewElementMatrix();
65 IssmDouble* basis = xNew<IssmDouble>(numnodes);
66- IssmDouble* B = xNew<IssmDouble>(2*numnodes);
67- IssmDouble* Bprime = xNew<IssmDouble>(2*numnodes);
68- IssmDouble D[2][2]={0.};
69+ IssmDouble* B = xNew<IssmDouble>(dim*numnodes);
70+ IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes);
71+ IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
72
73 /*Retrieve all inputs and parameters*/
74 basalelement->GetVerticesCoordinates(&xyz_list);
75@@ -160,10 +163,15 @@
76 vyaverage_input=basalelement->GetInput(VyEnum); _assert_(vyaverage_input);
77 }
78 else{
79- vxaverage_input=basalelement->GetInput(VxAverageEnum); _assert_(vxaverage_input);
80- vyaverage_input=basalelement->GetInput(VyAverageEnum); _assert_(vyaverage_input);
81+ if(dim==1){
82+ vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
83+ }
84+ if(dim==2){
85+ vxaverage_input=element->GetInput(VxAverageEnum); _assert_(vxaverage_input);
86+ vyaverage_input=element->GetInput(VyAverageEnum); _assert_(vyaverage_input);
87+ }
88 }
89- IssmDouble h=basalelement->CharacteristicLength();
90+ h=basalelement->CharacteristicLength();
91
92 /* Start looping on the number of gaussian points: */
93 Gauss* gauss=basalelement->NewGauss(2);
94@@ -172,64 +180,75 @@
95
96 basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
97 basalelement->NodalFunctions(basis,gauss);
98- D_scalar=gauss->weight*Jdet;
99-
100+
101 vxaverage_input->GetInputValue(&vx,gauss);
102- vyaverage_input->GetInputValue(&vy,gauss);
103 vxaverage_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
104- vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
105+ if(dim==2){
106+ vyaverage_input->GetInputValue(&vy,gauss);
107+ vyaverage_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
108+ }
109
110+ D_scalar=gauss->weight*Jdet;
111 TripleMultiply(basis,1,numnodes,1,
112 &D_scalar,1,1,0,
113 basis,1,numnodes,0,
114 &Ke->values[0],1);
115
116- GetB(B,element,xyz_list,gauss);
117- GetBprime(Bprime,element,xyz_list,gauss);
118+ GetB(B,basalelement,dim,xyz_list,gauss);
119+ GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
120
121 dvxdx=dvx[0];
122- dvydy=dvy[1];
123+ if(dim==2) dvydy=dvy[1];
124 D_scalar=dt*gauss->weight*Jdet;
125
126- D[0][0]=D_scalar*dvxdx;
127- D[1][1]=D_scalar*dvydy;
128- TripleMultiply(B,2,numnodes,1,
129- &D[0][0],2,2,0,
130- B,2,numnodes,0,
131+ D[0*dim+0]=D_scalar*dvxdx;
132+ if(dim==2) D[1*dim+1]=D_scalar*dvydy;
133+
134+ TripleMultiply(B,dim,numnodes,1,
135+ D,dim,dim,0,
136+ B,dim,numnodes,0,
137 &Ke->values[0],1);
138
139- D[0][0]=D_scalar*vx;
140- D[1][1]=D_scalar*vy;
141- TripleMultiply(B,2,numnodes,1,
142- &D[0][0],2,2,0,
143- Bprime,2,numnodes,0,
144+ D[0*dim+0]=D_scalar*vx;
145+ if(dim==2) D[1*dim+1]=D_scalar*vy;
146+
147+ TripleMultiply(B,dim,numnodes,1,
148+ D,dim,dim,0,
149+ Bprime,dim,numnodes,0,
150 &Ke->values[0],1);
151
152 if(stabilization==2){
153- /*Streamline upwinding*/
154- vel=sqrt(vx*vx+vy*vy)+1.e-8;
155- D[0][0]=h/(2*vel)*vx*vx;
156- D[1][0]=h/(2*vel)*vy*vx;
157- D[0][1]=h/(2*vel)*vx*vy;
158- D[1][1]=h/(2*vel)*vy*vy;
159+ if(dim==1){
160+ vel=fabs(vx)+1.e-8;
161+ D[0]=h/(2*vel)*vx*vx;
162+ }
163+ else{
164+ /*Streamline upwinding*/
165+ vel=sqrt(vx*vx+vy*vy)+1.e-8;
166+ D[0*dim+0]=h/(2*vel)*vx*vx;
167+ D[1*dim+0]=h/(2*vel)*vy*vx;
168+ D[0*dim+1]=h/(2*vel)*vx*vy;
169+ D[1*dim+1]=h/(2*vel)*vy*vy;
170+ }
171 }
172 else if(stabilization==1){
173- /*SSA*/
174 vxaverage_input->GetInputAverage(&vx);
175- vyaverage_input->GetInputAverage(&vy);
176- D[0][0]=h/2.0*fabs(vx);
177- D[0][1]=0.;
178- D[1][0]=0.;
179- D[1][1]=h/2.0*fabs(vy);
180+ if(dim==2) vyaverage_input->GetInputAverage(&vy);
181+ D[0*dim+0]=h/2.0*fabs(vx);
182+ if(dim==2) D[1*dim+1]=h/2.0*fabs(vy);
183 }
184 if(stabilization==1 || stabilization==2){
185- D[0][0]=D_scalar*D[0][0];
186- D[1][0]=D_scalar*D[1][0];
187- D[0][1]=D_scalar*D[0][1];
188- D[1][1]=D_scalar*D[1][1];
189- TripleMultiply(Bprime,2,numnodes,1,
190- &D[0][0],2,2,0,
191- Bprime,2,numnodes,0,
192+ if(dim==1) D[0]=D_scalar*D[0];
193+ else{
194+ D[0*dim+0]=D_scalar*D[0*dim+0];
195+ D[1*dim+0]=D_scalar*D[1*dim+0];
196+ D[0*dim+1]=D_scalar*D[0*dim+1];
197+ D[1*dim+1]=D_scalar*D[1*dim+1];
198+ }
199+
200+ TripleMultiply(Bprime,dim,numnodes,1,
201+ D,dim,dim,0,
202+ Bprime,dim,numnodes,0,
203 &Ke->values[0],1);
204 }
205
206@@ -240,15 +259,22 @@
207 xDelete<IssmDouble>(basis);
208 xDelete<IssmDouble>(B);
209 xDelete<IssmDouble>(Bprime);
210+ xDelete<IssmDouble>(D);
211 delete gauss;
212 if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
213 return Ke;
214 }/*}}}*/
215 ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
216
217+ /* Check if ice in element */
218+ if(!element->IsIceInElement()) return NULL;
219+
220 /*Intermediaries*/
221 int meshtype;
222 Element* basalelement;
223+ IssmDouble Jdet,dt;
224+ IssmDouble f,damage;
225+ IssmDouble* xyz_list = NULL;
226
227 /*Get basal element*/
228 element->FindParam(&meshtype,MeshTypeEnum);
229@@ -263,37 +289,35 @@
230 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
231 }
232
233- /*Intermediaries */
234- IssmDouble Jdet,dt;
235- IssmDouble f,damage;
236- IssmDouble* xyz_list = NULL;
237-
238 /*Fetch number of nodes and dof for this finite element*/
239- int numnodes = element->GetNumberOfNodes();
240+ int numnodes = basalelement->GetNumberOfNodes();
241
242 /*Initialize Element vector and other vectors*/
243- ElementVector* pe = element->NewElementVector();
244+ ElementVector* pe = basalelement->NewElementVector();
245 IssmDouble* basis = xNew<IssmDouble>(numnodes);
246
247 /*Retrieve all inputs and parameters*/
248- element->GetVerticesCoordinates(&xyz_list);
249- element->FindParam(&dt,TimesteppingTimeStepEnum);
250- this->CreateDamageFInput(element);
251- Input* damaged_input = element->GetInput(DamageDbarEnum); _assert_(damaged_input);
252- Input* damagef_input = element->GetInput(DamageFEnum); _assert_(damagef_input);
253+ basalelement->GetVerticesCoordinates(&xyz_list);
254+ basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
255+ this->CreateDamageFInput(basalelement);
256+ Input* damaged_input = basalelement->GetInput(DamageDEnum); _assert_(damaged_input);
257+ Input* damagef_input = basalelement->GetInput(DamageFEnum); _assert_(damagef_input);
258
259+
260 /* Start looping on the number of gaussian points: */
261- Gauss* gauss=element->NewGauss(2);
262+ Gauss* gauss=basalelement->NewGauss(2);
263 for(int ig=gauss->begin();ig<gauss->end();ig++){
264 gauss->GaussPoint(ig);
265
266- element->JacobianDeterminant(&Jdet,xyz_list,gauss);
267- element->NodalFunctions(basis,gauss);
268+ basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
269+ basalelement->NodalFunctions(basis,gauss);
270
271 damaged_input->GetInputValue(&damage,gauss);
272 damagef_input->GetInputValue(&f,gauss);
273
274- for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
275+ for(int i=0;i<numnodes;i++){
276+ pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
277+ }
278 }
279
280 /*Clean up and return*/
281@@ -303,7 +327,7 @@
282 delete gauss;
283 return pe;
284 }/*}}}*/
285-void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
286+void DamageEvolutionAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
287 /*Compute B matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2.
288 * For node i, Bi can be expressed in the actual coordinate system
289 * by:
290@@ -323,14 +347,15 @@
291
292 /*Build B: */
293 for(int i=0;i<numnodes;i++){
294- B[numnodes*0+i] = basis[i];
295- B[numnodes*1+i] = basis[i];
296+ for(int j=0;j<dim;j++){
297+ B[numnodes*j+i] = basis[i];
298+ }
299 }
300
301 /*Clean-up*/
302 xDelete<IssmDouble>(basis);
303 }/*}}}*/
304-void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
305+void DamageEvolutionAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
306 /*Compute B' matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2.
307 * For node i, Bi' can be expressed in the actual coordinate system
308 * by:
309@@ -345,13 +370,14 @@
310 int numnodes = element->GetNumberOfNodes();
311
312 /*Get nodal functions derivatives*/
313- IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
314+ IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
315 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
316
317 /*Build B': */
318 for(int i=0;i<numnodes;i++){
319- Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
320- Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
321+ for(int j=0;j<dim;j++){
322+ Bprime[numnodes*j+i] = dbasis[j*numnodes+i];
323+ }
324 }
325
326 /*Clean-up*/
327@@ -363,35 +389,46 @@
328 }/*}}}*/
329 void DamageEvolutionAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
330
331+ int meshtype;
332 IssmDouble max_damage;
333 int *doflist = NULL;
334+ Element* basalelement=NULL;
335
336+ element->FindParam(&meshtype,MeshTypeEnum);
337+ if(meshtype!=Mesh2DhorizontalEnum){
338+ if(!element->IsOnBed()) return;
339+ basalelement=element->SpawnBasalElement();
340+ }
341+ else{
342+ basalelement = element;
343+ }
344 /*Fetch number of nodes and dof for this finite element*/
345- int numnodes = element->GetNumberOfNodes();
346+ int numnodes = basalelement->GetNumberOfNodes();
347
348 /*Fetch dof list and allocate solution vector*/
349- element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
350- IssmDouble* values = xNew<IssmDouble>(numnodes);
351+ basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
352+ IssmDouble* newdamage = xNew<IssmDouble>(numnodes);
353
354 /*Get user-supplied max_damage: */
355- element->FindParam(&max_damage,DamageMaxDamageEnum);
356+ basalelement->FindParam(&max_damage,DamageMaxDamageEnum);
357
358 /*Use the dof list to index into the solution vector: */
359 for(int i=0;i<numnodes;i++){
360- values[i]=solution[doflist[i]];
361+ newdamage[i]=solution[doflist[i]];
362 /*Check solution*/
363- if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
364+ if(xIsNan<IssmDouble>(newdamage[i])) _error_("NaN found in solution vector");
365 /*Enforce D < max_damage and D > 0 */
366- if(values[i]>max_damage) values[i]=max_damage;
367- else if(values[i]<0.) values[i]=0.;
368+ if(newdamage[i]>max_damage) newdamage[i]=max_damage;
369+ else if(newdamage[i]<0.) newdamage[i]=0.;
370 }
371
372 /*Get all inputs and parameters*/
373- element->AddInput(DamageDbarEnum,values,P1Enum);
374+ basalelement->AddBasalInput(DamageDEnum,newdamage,P1Enum);
375
376 /*Free ressources:*/
377- xDelete<IssmDouble>(values);
378+ xDelete<IssmDouble>(newdamage);
379 xDelete<int>(doflist);
380+ if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
381 }/*}}}*/
382 void DamageEvolutionAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
383 /*Default, do nothing*/
384@@ -426,7 +463,7 @@
385 Input* sigma_xx_input = element->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input);
386 Input* sigma_xy_input = element->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input);
387 Input* sigma_yy_input = element->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input);
388- Input* damage_input = element->GetInput(DamageDbarEnum); _assert_(damage_input);
389+ Input* damage_input = element->GetInput(DamageDEnum); _assert_(damage_input);
390
391 /*retrieve the desired type of equivalent stress*/
392 element->FindParam(&equivstress,DamageEquivStressEnum);
393@@ -445,13 +482,12 @@
394 s_xy=sigma_xy/(1.-damage);
395 s_yy=sigma_yy/(1.-damage);
396
397- if(equivstress==1){ /* von Mises */
398- J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy);
399- Chi=sqrt(3.0)*J2s;
400+ if(equivstress==0){ /* von Mises */
401+ Chi=sqrt(s_xx*s_xx - s_xx*s_yy + s_yy*s_yy + 3*s_xy*s_xy);
402 }
403 Psi=Chi-stress_threshold;
404 PosPsi=max(Psi,0.);
405- NegPsi=max(-Psi,0.);
406+ NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
407
408 f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3);
409 }
410Index: ../trunk-jpl/src/c/cores/damage_core.cpp
411===================================================================
412--- ../trunk-jpl/src/c/cores/damage_core.cpp (revision 17390)
413+++ ../trunk-jpl/src/c/cores/damage_core.cpp (revision 17391)
414@@ -18,7 +18,6 @@
415 int numoutputs = 0;
416 char **requested_outputs = NULL;
417
418- if(VerboseSolution()) _printf0_(" computing damage\n");
419
420 //first recover parameters common to all solutions
421 femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
422@@ -32,6 +31,7 @@
423 ResetConstraintsx(femmodel);
424 }
425
426+ if(VerboseSolution()) _printf0_(" computing damage\n");
427 femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
428 solutionsequence_linear(femmodel);
429
430Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp
431===================================================================
432--- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17390)
433+++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17391)
434@@ -712,7 +712,7 @@
435 }
436
437 //First recover damage using the element: */
438- element->GetInputValue(&damage,node,DamageDbarEnum);
439+ element->GetInputValue(&damage,node,DamageDEnum);
440
441 //Recover our data:
442 parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);
Note: See TracBrowser for help on using the repository browser.