source: issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

Last change on this file was 28013, checked in by Mathieu Morlighem, 16 months ago

merged trunk-jpl and trunk for revision 28011
File size: 19.0 KB
Line 
1/*!\file SurfaceMassBalancex
2 * \brief: calculates SMB
3 */
4
5#include <config.h>
6#include "./SurfaceMassBalancex.h"
7#include "../../shared/shared.h"
8#include "../../toolkits/toolkits.h"
9#include "../modules.h"
10#include "../../classes/Inputs/TransientInput.h"
11#include "../../shared/Random/random.h"
12
13void SmbForcingx(FemModel* femmodel){/*{{{*/
14
15 // void SmbForcingx(smb,ni){
16 // INPUT parameters: ni: working size of arrays
17 // OUTPUT: mass-balance (m/yr ice): agd(NA)
18
19}/*}}}*/
20void SmbGradientsx(FemModel* femmodel){/*{{{*/
21
22 // void SurfaceMassBalancex(hd,agd,ni){
23 // INPUT parameters: ni: working size of arrays
24 // INPUT: surface elevation (m): hd(NA)
25 // OUTPUT: mass-balance (m/yr ice): agd(NA)
26 int v;
27 IssmDouble rho_water; // density of fresh water
28 IssmDouble rho_ice; // density of ice
29 IssmDouble yts; // conversion factor year to second
30
31 /*Loop over all the elements of this partition*/
32 for(Object* & object : femmodel->elements->objects){
33 Element* element=xDynamicCast<Element*>(object);
34
35 /*Allocate all arrays*/
36 int numvertices = element->GetNumberOfVertices();
37 IssmDouble* Href = xNew<IssmDouble>(numvertices); // reference elevation from which deviations are used to calculate the SMB adjustment
38 IssmDouble* Smbref = xNew<IssmDouble>(numvertices); // reference SMB to which deviations are added
39 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // Hs-SMB relation parameter
40 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // Hs-SMB relation paremeter
41 IssmDouble* s = xNew<IssmDouble>(numvertices); // surface elevation (m)
42 IssmDouble* smb = xNew<IssmDouble>(numvertices);
43
44 /*Recover SmbGradients*/
45 element->GetInputListOnVertices(Href,SmbHrefEnum);
46 element->GetInputListOnVertices(Smbref,SmbSmbrefEnum);
47 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
48 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
49
50 /*Recover surface elevation at vertices: */
51 element->GetInputListOnVertices(s,SurfaceEnum);
52
53 /*Get material parameters :*/
54 rho_ice=element->FindParam(MaterialsRhoIceEnum);
55 rho_water=element->FindParam(MaterialsRhoFreshwaterEnum);
56
57 /* Get constants */
58 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
59
60 // loop over all vertices
61 for(v=0;v<numvertices;v++){
62 if(Smbref[v]>0){
63 smb[v]=Smbref[v]+b_pos[v]*(s[v]-Href[v]);
64 }
65 else{
66 smb[v]=Smbref[v]+b_neg[v]*(s[v]-Href[v]);
67 }
68
69 } //end of the loop over the vertices
70
71 /*Add input to element and Free memory*/
72 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
73 xDelete<IssmDouble>(Href);
74 xDelete<IssmDouble>(Smbref);
75 xDelete<IssmDouble>(b_pos);
76 xDelete<IssmDouble>(b_neg);
77 xDelete<IssmDouble>(s);
78 xDelete<IssmDouble>(smb);
79 }
80
81}/*}}}*/
82void SmbGradientsElax(FemModel* femmodel){/*{{{*/
83
84 // void SurfaceMassBalancex(hd,agd,ni){
85 // INPUT parameters: ni: working size of arrays
86 // INPUT: surface elevation (m): hd(NA)
87 // OUTPUT: surface mass-balance (m/yr ice): agd(NA)
88 int v;
89
90 /*Loop over all the elements of this partition*/
91 for(Object* & object : femmodel->elements->objects){
92 Element* element=xDynamicCast<Element*>(object);
93
94 /*Allocate all arrays*/
95 int numvertices = element->GetNumberOfVertices();
96 IssmDouble* ela = xNew<IssmDouble>(numvertices); // Equilibrium Line Altitude (m a.s.l) to which deviations are used to calculate the SMB
97 IssmDouble* b_pos = xNew<IssmDouble>(numvertices); // SMB gradient above ELA (m ice eq. per m elevation change)
98 IssmDouble* b_neg = xNew<IssmDouble>(numvertices); // SMB gradient below ELA (m ice eq. per m elevation change)
99 IssmDouble* b_max = xNew<IssmDouble>(numvertices); // Upper cap on SMB rate (m/y ice eq.)
100 IssmDouble* b_min = xNew<IssmDouble>(numvertices); // Lower cap on SMB rate (m/y ice eq.)
101 IssmDouble* s = xNew<IssmDouble>(numvertices); // Surface elevation (m a.s.l.)
102 IssmDouble* smb = xNew<IssmDouble>(numvertices); // SMB (m/y ice eq.)
103
104 /*Recover ELA, SMB gradients, and caps*/
105 element->GetInputListOnVertices(ela,SmbElaEnum);
106 element->GetInputListOnVertices(b_pos,SmbBPosEnum);
107 element->GetInputListOnVertices(b_neg,SmbBNegEnum);
108 element->GetInputListOnVertices(b_max,SmbBMaxEnum);
109 element->GetInputListOnVertices(b_min,SmbBMinEnum);
110
111 /*Recover surface elevation at vertices: */
112 element->GetInputListOnVertices(s,SurfaceEnum);
113
114 /*Loop over all vertices, calculate SMB*/
115 for(v=0;v<numvertices;v++){
116 // if surface is above the ELA
117 if(s[v]>ela[v]){
118 smb[v]=b_pos[v]*(s[v]-ela[v]);
119 }
120 // if surface is below or equal to the ELA
121 else{
122 smb[v]=b_neg[v]*(s[v]-ela[v]);
123 }
124
125 // if SMB is larger than upper cap, set SMB to upper cap
126 if(smb[v]>b_max[v]){
127 smb[v]=b_max[v];
128 }
129 // if SMB is smaller than lower cap, set SMB to lower cap
130 if(smb[v]<b_min[v]){
131 smb[v]=b_min[v];
132 }
133 } //end of the loop over the vertices
134
135 /*Add input to element and Free memory*/
136 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
137 xDelete<IssmDouble>(ela);
138 xDelete<IssmDouble>(b_pos);
139 xDelete<IssmDouble>(b_neg);
140 xDelete<IssmDouble>(b_max);
141 xDelete<IssmDouble>(b_min);
142 xDelete<IssmDouble>(s);
143 xDelete<IssmDouble>(smb);
144
145 }
146
147}/*}}}*/
148void Smbarmax(FemModel* femmodel){/*{{{*/
149
150 /*Get time parameters*/
151 IssmDouble time,dt,starttime,tstep_arma;
152 femmodel->parameters->FindParam(&time,TimeEnum);
153 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
154 femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
155 femmodel->parameters->FindParam(&tstep_arma,SmbARMATimestepEnum);
156
157 /*Determine if this is a time step for the ARMA model*/
158 bool isstepforarma = false;
159
160 #ifndef _HAVE_AD_
161 if((fmod(time,tstep_arma)<fmod((time-dt),tstep_arma)) || (time<=starttime+dt) || tstep_arma==dt) isstepforarma = true;
162 #else
163 _error_("not implemented yet");
164 #endif
165
166 /*Load parameters*/
167 bool isstochastic;
168 bool issmbstochastic = false;
169 int M,N,arorder,maorder,numbasins,numparams,numbreaks,numelevbins,my_rank;
170 femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
171 femmodel->parameters->FindParam(&numparams,SmbNumParamsEnum);
172 femmodel->parameters->FindParam(&numbreaks,SmbNumBreaksEnum);
173 femmodel->parameters->FindParam(&arorder,SmbARMAarOrderEnum);
174 femmodel->parameters->FindParam(&maorder,SmbARMAmaOrderEnum);
175 femmodel->parameters->FindParam(&numelevbins,SmbNumElevationBinsEnum);
176 IssmDouble* datebreaks = NULL;
177 IssmDouble* arlagcoefs = NULL;
178 IssmDouble* malagcoefs = NULL;
179 IssmDouble* polyparams = NULL;
180 IssmDouble* lapserates = NULL;
181 IssmDouble* elevbins = NULL;
182 IssmDouble* refelevation = NULL;
183
184 femmodel->parameters->FindParam(&datebreaks,&M,&N,SmbARMAdatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1));
185 femmodel->parameters->FindParam(&polyparams,&M,&N,SmbARMApolyparamsEnum); _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams);
186 femmodel->parameters->FindParam(&arlagcoefs,&M,&N,SmbARMAarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder);
187 femmodel->parameters->FindParam(&malagcoefs,&M,&N,SmbARMAmalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder);
188 femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum); _assert_(M==numbasins); _assert_(N==numelevbins*12);
189 femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum); _assert_(M==numbasins); _assert_(N==(numelevbins-1)*12);
190 femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum); _assert_(M==numbasins);
191
192 femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
193 if(isstochastic){
194 int numstochasticfields;
195 int* stochasticfields;
196 femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
197 femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
198 for(int i=0;i<numstochasticfields;i++){
199 if(stochasticfields[i]==SMBarmaEnum) issmbstochastic = true;
200 }
201 xDelete<int>(stochasticfields);
202 }
203
204 /*Loop over each element to compute SMB at vertices*/
205 for(Object* &object:femmodel->elements->objects){
206 Element* element = xDynamicCast<Element*>(object);
207 /*Compute ARMA*/
208 element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,issmbstochastic,SMBarmaEnum);
209 /*Compute lapse rate adjustment*/
210 element->LapseRateBasinSMB(numelevbins,lapserates,elevbins,refelevation);
211 }
212
213 /*Cleanup*/
214 xDelete<IssmDouble>(arlagcoefs);
215 xDelete<IssmDouble>(malagcoefs);
216 xDelete<IssmDouble>(polyparams);
217 xDelete<IssmDouble>(datebreaks);
218 xDelete<IssmDouble>(lapserates);
219 xDelete<IssmDouble>(elevbins);
220 xDelete<IssmDouble>(refelevation);
221}/*}}}*/
222void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/
223
224 for(Object* & object : femmodel->elements->objects){
225 Element* element=xDynamicCast<Element*>(object);
226 element->Delta18oParameterization();
227 }
228
229}/*}}}*/
230void MungsmtpParameterizationx(FemModel* femmodel){/*{{{*/
231
232 for(Object* & object : femmodel->elements->objects){
233 Element* element=xDynamicCast<Element*>(object);
234 element->MungsmtpParameterization();
235 }
236
237}/*}}}*/
238void Delta18opdParameterizationx(FemModel* femmodel){/*{{{*/
239
240 for(Object* & object : femmodel->elements->objects){
241 Element* element=xDynamicCast<Element*>(object);
242 element->Delta18opdParameterization();
243 }
244
245}/*}}}*/
246void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/
247
248 // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){
249 // note "v" prefix means 12 monthly means, ie time dimension
250 // INPUT parameters: ni: working size of arrays
251 // INPUT: surface elevation (m): hd(NA)
252 // monthly mean surface sealevel temperature (degrees C): vTempsea(NA
253 // ,NTIME)
254 // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME)
255 // OUTPUT: mass-balance (m/yr ice): agd(NA)
256 // mean annual surface temperature (degrees C): Tsurf(NA)
257
258 int it, jj, itm;
259 IssmDouble DT = 0.02, sigfac, snormfac;
260 IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day
261 IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm
262 IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day
263 IssmDouble siglimc, siglim0, siglim0c;
264 IssmDouble tstep, tsint, tint, tstepc;
265 int NPDMAX = 1504, NPDCMAX = 1454;
266 //IssmDouble pdds[NPDMAX]={0};
267 //IssmDouble pds[NPDCMAX]={0};
268 IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow
269 IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
270 IssmDouble tstar; // monthly mean surface temp
271
272 bool ismungsm;
273 bool issetpddfac;
274
275 IssmDouble *pdds = NULL;
276 IssmDouble *pds = NULL;
277 Element *element = NULL;
278
279 pdds=xNew<IssmDouble>(NPDMAX+1);
280 pds=xNew<IssmDouble>(NPDCMAX+1);
281
282 // Get ismungsm parameter
283 femmodel->parameters->FindParam(&ismungsm,SmbIsmungsmEnum);
284
285 // Get issetpddfac parameter
286 femmodel->parameters->FindParam(&issetpddfac,SmbIssetpddfacEnum);
287
288 /* initialize PDD (creation of a lookup table)*/
289 tstep = 0.1;
290 tsint = tstep*0.5;
291 sigfac = -1.0/(2.0*pow(signorm,2));
292 snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
293 siglim = 2.5*signorm;
294 siglimc = 2.5*signormc;
295 siglim0 = siglim/DT + 0.5;
296 siglim0c = siglimc/DT + 0.5;
297 PDup = siglimc+PDCUT;
298
299 itm = reCast<int,IssmDouble>((2*siglim/DT + 1.5));
300
301 if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp");
302 for(it = 0; it < itm; it++){
303 // tstar = REAL(it)*DT-siglim;
304 tstar = it*DT-siglim;
305 tint = tsint;
306 pddt = 0.;
307 for ( jj = 0; jj < 600; jj++){
308 if (tint > (tstar+siglim)){break;}
309 pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep;
310 tint = tint+tstep;
311 }
312 pdds[it] = pddt*snormfac;
313 }
314 pdds[itm+1] = siglim + DT;
315
316 //*********compute PD(T) : snow/precip fraction. precipitation falling as snow
317 tstepc = 0.1;
318 tsint = PDCUT-tstepc*0.5;
319 signormc = signorm - 0.5;
320 sigfac = -1.0/(2.0*pow(signormc,2));
321 snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0)));
322 siglimc = 2.5*signormc ;
323 itm = reCast<int,IssmDouble>((PDCUT+2.*siglimc)/DT + 1.5);
324 if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com");
325 for(it = 0; it < itm; it++ ){
326 tstar = it*DT-siglimc;
327 // tstar = REAL(it)*DT-siglimc;
328 tint = tsint; // start against upper bound
329 pd = 0.;
330 for (jj = 0; jj < 600; jj++){
331 if (tint<(tstar-siglimc)) {break;}
332 pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc;
333 tint = tint-tstepc;
334 }
335 pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction
336 }
337 pds[itm+1] = 0.;
338 // *******END initialize PDD
339
340 for(Object* & object : femmodel->elements->objects){
341 element=xDynamicCast<Element*>(object);
342 element->PositiveDegreeDay(pdds,pds,signorm,ismungsm,issetpddfac);
343 }
344 /*free ressouces: */
345 xDelete<IssmDouble>(pdds);
346 xDelete<IssmDouble>(pds);
347}/*}}}*/
348void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/
349
350 bool isfirnwarming;
351 femmodel->parameters->FindParam(&isfirnwarming,SmbIsfirnwarmingEnum);
352
353 for(Object* & object : femmodel->elements->objects){
354 Element* element=xDynamicCast<Element*>(object);
355 element->PositiveDegreeDaySicopolis(isfirnwarming);
356 }
357
358}/*}}}*/
359void SmbHenningx(FemModel* femmodel){/*{{{*/
360
361 /*Intermediaries*/
362 IssmDouble z_critical = 1675.;
363 IssmDouble dz = 0;
364 IssmDouble a = -15.86;
365 IssmDouble b = 0.00969;
366 IssmDouble c = -0.235;
367 IssmDouble f = 1.;
368 IssmDouble g = -0.0011;
369 IssmDouble h = -1.54e-5;
370 IssmDouble smb,smbref,anomaly,yts,z;
371
372 /* Get constants */
373 femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
374 /*iomodel->FindConstant(&yts,"md.constants.yts");*/
375 /*this->parameters->FindParam(&yts,ConstantsYtsEnum);*/
376 /*Mathieu original*/
377 /*IssmDouble smb,smbref,z;*/
378
379 /*Loop over all the elements of this partition*/
380 for(Object* & object : femmodel->elements->objects){
381 Element* element=xDynamicCast<Element*>(object);
382
383 /*Get reference SMB (uncorrected) and allocate all arrays*/
384 int numvertices = element->GetNumberOfVertices();
385 IssmDouble* surfacelist = xNew<IssmDouble>(numvertices);
386 IssmDouble* smblistref = xNew<IssmDouble>(numvertices);
387 IssmDouble* smblist = xNew<IssmDouble>(numvertices);
388 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
389 element->GetInputListOnVertices(smblistref,SmbSmbrefEnum);
390
391 /*Loop over all vertices of element and correct SMB as a function of altitude z*/
392 for(int v=0;v<numvertices;v++){
393
394 /*Get vertex elevation, anoma smb*/
395 z = surfacelist[v];
396 anomaly = smblistref[v];
397
398 /* Henning edited acc. to Riannes equations*/
399 /* Set SMB maximum elevation, if dz = 0 -> z_critical = 1675 */
400 z_critical = z_critical + dz;
401
402 /* Calculate smb acc. to the surface elevation z */
403 if(z<z_critical){
404 smb = a + b*z + c;
405 }
406 else{
407 smb = (a + b*z)*(f + g*(z-z_critical) + h*(z-z_critical)*(z-z_critical)) + c;
408 }
409
410 /* Compute smb including anomaly,
411 correct for number of seconds in a year [s/yr]*/
412 smb = smb/yts + anomaly;
413
414 /*Update array accordingly*/
415 smblist[v] = smb;
416
417 }
418
419 /*Add input to element and Free memory*/
420 element->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
421 xDelete<IssmDouble>(surfacelist);
422 xDelete<IssmDouble>(smblistref);
423 xDelete<IssmDouble>(smblist);
424 }
425
426}/*}}}*/
427void SmbComponentsx(FemModel* femmodel){/*{{{*/
428
429 // void SmbComponentsx(acc,evap,runoff,ni){
430 // INPUT parameters: ni: working size of arrays
431 // INPUT: surface accumulation (m/yr water equivalent): acc
432 // surface evaporation (m/yr water equivalent): evap
433 // surface runoff (m/yr water equivalent): runoff
434 // OUTPUT: mass-balance (m/yr ice): agd(NA)
435
436 /*Loop over all the elements of this partition*/
437 for(Object* & object : femmodel->elements->objects){
438 Element* element=xDynamicCast<Element*>(object);
439
440 /*Allocate all arrays*/
441 int numvertices = element->GetNumberOfVertices();
442 IssmDouble* acc = xNew<IssmDouble>(numvertices);
443 IssmDouble* evap = xNew<IssmDouble>(numvertices);
444 IssmDouble* runoff = xNew<IssmDouble>(numvertices);
445 IssmDouble* smb = xNew<IssmDouble>(numvertices);
446
447 /*Recover Smb Components*/
448 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
449 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
450 element->GetInputListOnVertices(runoff,SmbRunoffEnum);
451
452 // loop over all vertices
453 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-runoff[v];
454
455 /*Add input to element and Free memory*/
456 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
457 xDelete<IssmDouble>(acc);
458 xDelete<IssmDouble>(evap);
459 xDelete<IssmDouble>(runoff);
460 xDelete<IssmDouble>(smb);
461 }
462
463}/*}}}*/
464void SmbMeltComponentsx(FemModel* femmodel){/*{{{*/
465
466 // void SmbMeltComponentsx(acc,evap,melt,refreeze,ni){
467 // INPUT parameters: ni: working size of arrays
468 // INPUT: surface accumulation (m/yr water equivalent): acc
469 // surface evaporation (m/yr water equivalent): evap
470 // surface melt (m/yr water equivalent): melt
471 // refreeze of surface melt (m/yr water equivalent): refreeze
472 // OUTPUT: mass-balance (m/yr ice): agd(NA)
473
474 /*Loop over all the elements of this partition*/
475 for(Object* & object : femmodel->elements->objects){
476 Element* element=xDynamicCast<Element*>(object);
477
478 /*Allocate all arrays*/
479 int numvertices = element->GetNumberOfVertices();
480 IssmDouble* acc = xNew<IssmDouble>(numvertices);
481 IssmDouble* evap = xNew<IssmDouble>(numvertices);
482 IssmDouble* melt = xNew<IssmDouble>(numvertices);
483 IssmDouble* refreeze = xNew<IssmDouble>(numvertices);
484 IssmDouble* smb = xNew<IssmDouble>(numvertices);
485
486 /*Recover Smb Components*/
487 element->GetInputListOnVertices(acc,SmbAccumulationEnum);
488 element->GetInputListOnVertices(evap,SmbEvaporationEnum);
489 element->GetInputListOnVertices(melt,SmbMeltEnum);
490 element->GetInputListOnVertices(refreeze,SmbRefreezeEnum);
491
492 // loop over all vertices
493 for(int v=0;v<numvertices;v++) smb[v]=acc[v]-evap[v]-melt[v]+refreeze[v];
494
495 /*Add input to element and Free memory*/
496 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
497 xDelete<IssmDouble>(acc);
498 xDelete<IssmDouble>(evap);
499 xDelete<IssmDouble>(melt);
500 xDelete<IssmDouble>(refreeze);
501 xDelete<IssmDouble>(smb);
502 }
503
504}/*}}}*/
505void SmbDebrisEvattx(FemModel* femmodel){/*{{{*/
506 for(Object* & object : femmodel->elements->objects){
507 Element* element=xDynamicCast<Element*>(object);
508 element->SmbDebrisEvatt();
509 }
510}/*}}}*/
511void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
512
513 for(Object* & object : femmodel->elements->objects){
514 Element* element=xDynamicCast<Element*>(object);
515 element->SmbGradCompParameterization();
516 }
517
518}/*}}}*/
519#ifdef _HAVE_SEMIC_
520void SmbSemicx(FemModel* femmodel,int ismethod){/*{{{*/
521
522 for(Object* & object : femmodel->elements->objects){
523 Element* element=xDynamicCast<Element*>(object);
524 if (ismethod == 1) element->SmbSemicTransient(); // Inwoo's version.
525 else element->SmbSemic(); // original SmbSEMIC
526 }
527
528}/*}}}*/
529#else
530void SmbSemicx(FemModel* femmodel){_error_("SEMIC not installed");}
531#endif //_HAVE_SEMIC_
Note: See TracBrowser for help on using the repository browser.