proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
ElastoPlastic.h
Go to the documentation of this file.
1#ifndef ELASTOPLASTIC_H
2#define ELASTOPLASTIC_H
3#include <cmath>
4#include <iostream>
5#include <iomanip>
6#include <vector>
7//extern "C"
8//{
9#include PROTEUS_LAPACK_H
10#include PROTEUS_BLAS_H
11//}
12#include "proteus_blas.h"
13#include "CompKernel.h"
14#include "ModelFactory.h"
15#include "mohrCoulomb.h"
16#include "mohrCoulomb2.h"
17#include "vonMises.h"
18#include "tresca.h"
19#include "druckerPrager.h"
20#include "mcdp.h"
21#include "../mprans/ArgumentsDict.h"
22#include "xtensor-python/pyarray.hpp"
23#include "pybind11/pybind11.h"
24
25namespace py = pybind11;
26
27namespace proteus
28{
29 class ElastoPlastic_base
30 {
31 public:
32 virtual ~ElastoPlastic_base(){}
33 virtual void calculateResidual(arguments_dict& args)=0;
34 virtual void calculateJacobian(arguments_dict& args)=0;
35 };
36
37 template<class CompKernelType,
38 int nSpace,
39 int nQuadraturePoints_element,
40 int nDOF_mesh_trial_element,
41 int nDOF_trial_element,
42 int nDOF_test_element,
43 int nQuadraturePoints_elementBoundary>
44 class ElastoPlastic : public ElastoPlastic_base
45 {
46 public:
47 CompKernelType ck;
48
49 const int nDOF_test_X_trial_element;
50
51 EIndex<nSpace> ei;
52
53 const int X,Y,Z,
54 XX,XY,XZ,
55 YX,YY,YZ,
56 ZX,ZY,ZZ,
57 sXX,sXY,sXZ,
58 sYX,sYY,sYZ,
59 sZX,sZY,sZZ,
60 nSymTen,
61 XHX,XHY,
62 YHX,YHY,
63 ZHX,ZHY,
64 HXHX,HXHY,
65 HYHX,HYHY;
66
67 ElastoPlastic():
68 ck(),
69 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
70 X(ei.X),
71 Y(ei.Y),
72 Z(ei.Z),
73 XX(ei.XX),XY(ei.XY),XZ(ei.XZ),
74 YX(ei.YX),YY(ei.YY),YZ(ei.YZ),
75 ZX(ei.ZX),ZY(ei.ZY),ZZ(ei.ZZ),
76 sXX(ei.sXX),sXY(ei.sXY),sXZ(ei.sXZ),
77 sYX(ei.sYX),sYY(ei.sYY),sYZ(ei.sYZ),
78 sZX(ei.sZX),sZY(ei.sZY),sZZ(ei.sZZ),
79 nSymTen(ei.nSymTen),
80 XHX(ei.XHX),XHY(ei.XHY),
81 YHX(ei.YHX),YHY(ei.YHY),
82 ZHX(ei.ZHX),ZHY(ei.ZHY),
83 HXHX(ei.HXHX),HXHY(ei.HXHY),
84 HYHX(ei.HYHX),HYHY(ei.HYHY)
85 {}
86
87 inline void calculateStrain(double* D, double* strain)
88 {
89 //Voigt notation from Belytschko, Liu, Moran
90 strain[sXX] = D[XX];//du/dx
91 strain[sYY] = D[YY];//dv/dy
92 strain[sZZ] = D[ZZ];//dw/dz
93 strain[sYZ] = D[YZ]+D[ZY];//(dv/dz + dz/dy)
94 strain[sXZ] = D[XZ]+D[ZX];//(du/dz + dw/dx)
95 strain[sXY] = D[XY]+D[YX];//(du/dy + dv/dx)
96 }
97
98 inline void elasticStress(const double* C, const double* strain, double* stress)
99 {
100 stress[sXX] = C[0]*strain[sXX] + C[1]*(strain[sYY] +strain[sZZ]);
101 stress[sYY] = C[0]*strain[sYY] + C[1]*(strain[sXX] +strain[sZZ]);
102 stress[sZZ] = C[0]*strain[sZZ] + C[1]*(strain[sXX] +strain[sYY]);
103 stress[sYZ] = C[sYZ+sYZ*nSymTen]*strain[sYZ];
104 stress[sXZ] = C[sYZ+sYZ*nSymTen]*strain[sXZ];
105 stress[sXY] = C[sYZ+sYZ*nSymTen]*strain[sXY];
106 }
107
108 inline void elasticStrain(const double* Cinv, const double* stress, double* strain)
109 {
110 strain[sXX] = Cinv[0]*stress[sXX] + Cinv[1]*(stress[sYY] +stress[sZZ]);
111 strain[sYY] = Cinv[0]*stress[sYY] + Cinv[1]*(stress[sXX] +stress[sZZ]);
112 strain[sZZ] = Cinv[0]*stress[sZZ] + Cinv[1]*(stress[sXX] +stress[sYY]);
113 strain[sYZ] = Cinv[sYZ+sYZ*nSymTen]*stress[sYZ];
114 strain[sXZ] = Cinv[sYZ+sYZ*nSymTen]*stress[sXZ];
115 strain[sXY] = Cinv[sYZ+sYZ*nSymTen]*stress[sXY];
116 }
117
118 inline void elasticModuli(const double& E,const double& nu,double* C,double* Cinv)
119 {
120 for (int i=0;i<nSymTen;i++)
121 for (int j=0;j<nSymTen;j++)
122 {
123 C[i+nSymTen*j] = 0.0;
124 Cinv[i+nSymTen*j] = 0.0;
125 }
126 C[sXX+sXX*nSymTen] = (E/(1.0+nu))*(1.0 + (nu/(1.0-2.0*nu)));
127 C[sXX+sYY*nSymTen] = (E/(1.0+nu))*(nu/(1.0-2.0*nu));
128 C[sXX+sZZ*nSymTen] = (E/(1.0+nu))*(nu/(1.0-2.0*nu));
129 C[sYY+sXX*nSymTen] = (E/(1.0+nu))*(nu/(1.0-2.0*nu));
130 C[sYY+sYY*nSymTen] = (E/(1.0+nu))*(1.0 + (nu/(1.0-2.0*nu)));
131 C[sYY+sZZ*nSymTen] = (E/(1.0+nu))*(nu/(1.0-2.0*nu));
132 C[sZZ+sXX*nSymTen] = (E/(1.0+nu))*(nu/(1.0-2.0*nu));
133 C[sZZ+sYY*nSymTen] = (E/(1.0+nu))*(nu/(1.0-2.0*nu));
134 C[sZZ+sZZ*nSymTen] = (E/(1.0+nu))*(1.0 + (nu/(1.0-2.0*nu)));
135 C[sYZ+sYZ*nSymTen] = (E/(1.0+nu))*0.5;
136 C[sXZ+sXZ*nSymTen] = (E/(1.0+nu))*0.5;
137 C[sXY+sXY*nSymTen] = (E/(1.0+nu))*0.5;
138
139 Cinv[sXX+sXX*nSymTen] = 1.0/E;
140 Cinv[sXX+sYY*nSymTen] = -nu/E;
141 Cinv[sXX+sZZ*nSymTen] = -nu/E;
142 Cinv[sYY+sXX*nSymTen] = -nu/E;
143 Cinv[sYY+sYY*nSymTen] = 1.0/E;
144 Cinv[sYY+sZZ*nSymTen] = -nu/E;
145 Cinv[sZZ+sXX*nSymTen] = -nu/E;
146 Cinv[sZZ+sYY*nSymTen] = -nu/E;
147 Cinv[sZZ+sZZ*nSymTen] = 1.0/E;
148 Cinv[sYZ+sYZ*nSymTen] = ((1.0+nu)/E)*2.0;
149 Cinv[sXZ+sXZ*nSymTen] = ((1.0+nu)/E)*2.0;
150 Cinv[sXY+sXY*nSymTen] = ((1.0+nu)/E)*2.0;
151 //cek debug
152 double I[nSymTen*nSymTen];
153 for (int i=0;i<nSymTen;i++)
154 {
155 for (int j=0;j<nSymTen;j++)
156 {
157 I[i+j*nSymTen] = 0.0;
158 for(int k=0;k<nSymTen;k++)
159 I[i+j*nSymTen] += C[i+k*nSymTen]*Cinv[k+j*nSymTen];
160 if(i==j)
161 if(fabs(1.0-I[i+j*nSymTen]) > 1.0e-12)
162 std::cout<<"I["<<i<<"]["<<j<<"]="<<I[i+j*nSymTen]<<std::endl;
163 if(i!=j)
164 if(fabs(I[i+j*nSymTen]) > 1.0e-12)
165 std::cout<<"I["<<i<<"]["<<j<<"]="<<I[i+j*nSymTen]<<std::endl;
166 }
167 }
168 }
169
170 inline void evaluateConstitutiveModel(const double* materialProperties, const double* stress,
171 double& f, double* df, double* r, double* dr, double& stress_3)
172 {
173 double g;
174 /* std::cout<<materialProperties[5]<<"phi \n" */
175 /* <<materialProperties[7]<<"psi \n" */
176 /* <<materialProperties[6]<<"c \n"<<std::endl; */
177 /* mohrCoulomb(materialProperties[5],//phi */
178 /* materialProperties[7],//psi */
179 /* materialProperties[6],//c */
180 /* stress, */
181 /* f, */
182 /* df, */
183 /* g, */
184 /* r, */
185 /* dr); */
186 /* mohrCoulomb2(materialProperties[5],//phi */
187 /* materialProperties[7],//psi */
188 /* materialProperties[6],//c */
189 /* stress, */
190 /* f, */
191 /* df, */
192 /* g, */
193 /* r, */
194 /* dr); */
195 /* if (! (std::isfinite(f) && */
196 /* std::isfinite(df[sXX]) && std::isfinite(df[sYY]) && std::isfinite(df[sZZ]) && std::isfinite(df[sYZ]) && std::isfinite(df[sXZ]) && std::isfinite(df[sXY]) && */
197 /* std::isfinite(r[sXX]) && std::isfinite(r[sYY]) && std::isfinite(r[sZZ]) && std::isfinite(r[sYZ]) && std::isfinite(r[sXZ]) && std::isfinite(r[sXY]))) */
198 /* { */
199 /* std::cout<<"f "<<f<<std::endl */
200 /* <<"stress "<<stress[sXX]<<'\t'<<stress[sYY]<<'\t'<<stress[sZZ]<<'\t'<<stress[sYZ]<<'\t'<<stress[sXZ]<<'\t'<<stress[sXY]<<std::endl */
201 /* <<"df "<<df[sXX]<<'\t'<<df[sYY]<<'\t'<<df[sZZ]<<'\t'<<df[sYZ]<<'\t'<<df[sXZ]<<'\t'<<df[sXY]<<std::endl */
202 /* <<"r "<<r[sXX]<<'\t'<<r[sYY]<<'\t'<<r[sZZ]<<'\t'<<r[sYZ]<<'\t'<<r[sXZ]<<'\t'<<r[sXY]<<std::endl; */
203 /* } */
204 /* vonMises(materialProperties[5],//phi */
205 /* materialProperties[7],//psi */
206 /* materialProperties[6],//c */
207 /* stress, */
208 /* f, */
209 /* df, */
210 /* g, */
211 /* r, */
212 /* dr); */
213 /* tresca(materialProperties[5],//phi */
214 /* materialProperties[7],//psi */
215 /* materialProperties[6],//c */
216 /* stress, */
217 /* f, */
218 /* df, */
219 /* g, */
220 /* r, */
221 /* dr); */
222 /* druckerPrager(materialProperties[5],//phi */
223 /* materialProperties[7],//psi */
224 /* materialProperties[6],//c */
225 /* stress, */
226 /* f, */
227 /* df, */
228 /* g, */
229 /* r, */
230 /* dr); */
231 mcdp(materialProperties[5],//phi
232 materialProperties[7],//psi
233 materialProperties[6],//c
234 stress,
235 f,
236 df,
237 g,
238 r,
239 dr,
240 stress_3);
241 }
242
243 inline void differenceJacobian(const double* materialProperties,const double* stress,const double& f,double *df, const double* r, double* dr)
244 {
245 double f_delta,
246 df_delta[nSymTen],
247 stress_plus_delta[nSymTen],
248 r_delta[nSymTen],
249 dr_delta[nSymTen*nSymTen],stress_3;
250 for (int i=0;i<nSymTen;i++)
251 {
252 stress_plus_delta[i] = stress[i];
253 }
254 for (int j=0;j<nSymTen;j++)
255 {
256 double Delta_stress = 1.0e-8*stress[j] + 1.0e-8;
257 stress_plus_delta[j] += Delta_stress;
258 evaluateConstitutiveModel(materialProperties,stress_plus_delta,f_delta,df_delta,r_delta,dr_delta,stress_3);
259 df[j] = (f_delta - f)/Delta_stress;
260 for(int k=0;k<nSymTen;k++)
261 {
262 dr[k+nSymTen*j] = (r_delta[k] - r[k])/Delta_stress;
263 }
264 stress_plus_delta[j] = stress[j];
265 }
266 }
267 inline void evaluateCoefficients(int usePicard,
268 const double pore_pressure,
269 const double* materialProperties,
270 const double* strain0,
271 const double* strain_last,
272 const double* Delta_strain,
273 const double* plasticStrain_last,
274 double* plasticStrain,
275 double* stress,
276 double* dstress)
277 {
278 int its=0, maxIts=25;
279 double E=materialProperties[0];
280 const double nu=materialProperties[1],
281 phi=materialProperties[5];
282 double f,
283 df[nSymTen],
284 r[nSymTen],
285 dr[nSymTen*nSymTen],
286 r0[nSymTen],
287 C[nSymTen*nSymTen],
288 Cinv[nSymTen*nSymTen],
289 A[nSymTen*nSymTen],
290 B[(nSymTen+1)*(nSymTen+1)],
291 a[nSymTen],
292 Delta_lambda,
293 Delta_lambda_old,
294 Delta_lambda_incr,
295 Delta_stress[nSymTen],
296 Delta_stress_old[nSymTen],
297 stress_old[nSymTen],
298 Delta_stress_incr[nSymTen],
299 minusDelta_plasticStrain[nSymTen],
300 effectiveStrain[nSymTen],
301 aNorm,
302 stressNorm,
303 f_atol=1.0e-6,
304 a_atol=1.0e-6,
305 dfA[nSymTen],
306 Aa[nSymTen],
307 Ar[nSymTen],
308 dfAa,
309 dfAr,
310 WORK[nSymTen],
311 R[nSymTen+1],
312 WORKR[nSymTen+1],
313 stress_3_init,stress_3;
314 PROTEUS_LAPACK_INTEGER N=nSymTen,NR=nSymTen+1,
315 INFO=0,
316 NRHS=1,
317 IPIV[nSymTen],
318 JPIV[nSymTen],
319 LWORK=nSymTen,
320 INFOR=0,
321 IPIVR[nSymTen+1],
322 JPIVR[nSymTen+1],
323 LWORKR=nSymTen+1;
324 double SCALE,SCALER;
325 char TRANS='N';
326 bool useSemiImplicit=false;
327 bool useDifferenceJacobian=false;
328 bool useFullPivoting=false;
329 bool useDumbNewton=false;
330 bool useContinuumTangentModulus=false;
331 bool useLineSearch=false;
332 bool useInnerPicard=false;
333 bool useOuterPicard=bool(usePicard);
334 bool predictorPhase = false;
335 int predictorPhaseIts=5;
336 std::vector<double> fhist,ahist;
337 //get modulus from initial strain
338 elasticModuli(E,nu,C,Cinv);
339 elasticStress(C,strain0,stress);
340 evaluateConstitutiveModel(materialProperties,stress,f,df,r,dr,stress_3_init);//cek hack, this is just to get stress_3_init
341 const double n_e = 0.6,
342 P_a = materialProperties[12],
343 K_i = 500.0;
344 if (P_a > 0.0)//cek hack, if not a soil then set P_a=0.0
345 {
346 E =K_i*P_a*pow(fmax(stress_3_init,P_a)/P_a,n_e);
347 //std::cout<<"E "<<E<<'\t'<<K_i<<'\t'<<P_a<<'\t'<<n_e<<'\t'<<stress_3_init<<std::endl;
348 elasticModuli(E,nu,C,Cinv);
349 }
350 //
351 //fully implicit Backward Euler integration
352 //
353 //input is strain_last, Delta_strain, and plasticStrain_last
354 //output is stress,dstress= dstress/dDelta_strain, and plasticStrain
355 //
356 //get stress at last step and evaluate r0 there for semi-implicit scheme because it lies on the yield surface
357 //
358 for (int i=0;i<nSymTen;i++)
359 {
360 effectiveStrain[i] = strain_last[i] - plasticStrain_last[i];
361 }
362 elasticStress(C,effectiveStrain,stress);
363 evaluateConstitutiveModel(materialProperties,stress,f,df,r,dr,stress_3);
364 for (int i=0;i<nSymTen;i++)
365 {
366 r0[i]=r[i];
367 }
368 //
369 //get elastic predictor stress
370 //
371 for (int i=0;i<nSymTen;i++)
372 {
373 Delta_lambda = 0.0;
374 Delta_stress[i] = 0.0;//tricky, this is not (stress - stress_last) but rather (stress - stress_elasticPredictor) which is 0 initially
375 minusDelta_plasticStrain[i] = 0.0;
376 plasticStrain[i] = plasticStrain_last[i];
377 effectiveStrain[i] = strain_last[i] + Delta_strain[i] - plasticStrain_last[i];
378 a[i] = 0.0;
379 }
380 elasticStress(C,effectiveStrain,stress);
381 evaluateConstitutiveModel(materialProperties,stress,f,df,r,dr,stress_3);
382 if (useDifferenceJacobian)
383 {
384 differenceJacobian(materialProperties,stress,f,df,r,dr);//recalculate df and dr
385 }
386 if (useSemiImplicit)//overwrite r and dr
387 {
388 for (int i=0;i<nSymTen;i++)
389 {
390 r[i]=r0[i];
391 for(int j=0;j<nSymTen;j++)
392 {
393 dr[i+j*nSymTen] = 0.0;
394 }
395 }
396 }
397 aNorm=0.0;//a = minusDelta_plasticStrain + Delta_lambda*r (==0 initially)
398 double f_last=f;
399 //main Newton iteration loop
400 while ((f >= f_atol || aNorm >= a_atol) && its < maxIts) //since aNorm=0 this will be skipped if f < f_atol i.e. within the yield surface up to the tolerance)
401 {
402 //fhist.push_back(f);
403 //ahist.push_back(aNorm);
404 //std::cout<<"plastic it "<<its<<" aNorm "<<aNorm<<" f "<<f<<std::endl;
405 //
406 //calculate Newton increment lambda_incr using block gaussian elimination
407 //
408 for(int i=0;i<nSymTen;i++)
409 {
410 for(int j=0;j<nSymTen;j++)
411 {
412 //
413 //set to A^{-1} and factor below
414 //
415 A[i+j*nSymTen] = Cinv[i+j*nSymTen];
416 if (!useInnerPicard)
417 A[i+j*nSymTen] += Delta_lambda*dr[i+j*nSymTen];
418 }
419 //dfA[i] = df[i];
420 Aa[i] = a[i];
421 Ar[i] = r[i];
422 //load the whole Jacobian if true to do brute force factorization instead of block factor
423 if (useDumbNewton)
424 {
425 for(int j=0;j<nSymTen;j++)
426 {
427 B[i+j*(nSymTen+1)] = A[i+j*nSymTen];
428 }
429 B[i+nSymTen*(nSymTen+1)] = r[i];
430 B[nSymTen + i*(nSymTen+1)] = df[i];
431 R[i] = -a[i];
432 }
433 }
434 if (useDumbNewton)
435 {
436 B[nSymTen + nSymTen*(nSymTen+1)] = 0.0;
437 R[nSymTen] = -f;
438 }
439 TRANS='N';
440 if(useFullPivoting)
441 {
442 dgetc2_(&N,A,&N,IPIV,JPIV,&INFO);
443 dgesc2_(&N,A,&N,Aa,IPIV,JPIV,&SCALE);
444 dgesc2_(&N,A,&N,Ar,IPIV,JPIV,&SCALE);
445 }
446 else
447 {
448 dgetrf_(&N,&N,A,&N,IPIV,&INFO);
449 dgetrs_(&TRANS,&N,&NRHS,A,&N,IPIV,Aa,&N,&INFO);
450 dgetrs_(&TRANS,&N,&NRHS,A,&N,IPIV,Ar,&N,&INFO);
451 }
452 //
453 //calculate increment in Delta_lamba
454 //
455 dfAa = 0.0;
456 dfAr = 0.0;
457 for(int i=0;i<nSymTen;i++)
458 {
459 dfAa += df[i]*Aa[i];
460 dfAr += df[i]*Ar[i];
461 }
462 if (fabs(dfAr) < 1.0e-8)
463 dfAr = copysign(1.0e-8,dfAr);
464 Delta_lambda_incr = (f-dfAa)/dfAr;
465 //calculate increment in Delta_stress
466 for(int i=0;i<nSymTen;i++)
467 {
468 Delta_stress_incr[i] = - Aa[i] - Delta_lambda_incr*Ar[i];
469 }
470 //overwrite increments if using brute force Newton
471 if(useDumbNewton)
472 {
473 dgetc2_(&NR,B,&NR,IPIVR,JPIVR,&INFOR);
474 dgesc2_(&NR,B,&NR,R,IPIVR,JPIVR,&SCALER);
475 for(int i=0;i<nSymTen;i++)
476 {
477 Delta_stress_incr[i] = R[i];
478 }
479 Delta_lambda_incr = R[nSymTen];
480 }
481 //increment unknowns
482 Delta_lambda_old = Delta_lambda;
483 Delta_lambda += Delta_lambda_incr;
484 for(int i=0;i<nSymTen;i++)
485 {
486 Delta_stress_old[i] = Delta_stress[i];
487 Delta_stress[i] += Delta_stress_incr[i];
488 stress_old[i] = stress[i];
489 stress[i] += Delta_stress_incr[i];
490 }
491 evaluateConstitutiveModel(materialProperties,stress,f,df,r,dr,stress_3);
492 if (useDifferenceJacobian)
493 {
494 differenceJacobian(materialProperties,stress,f,df,r,dr);
495 }
496 //reset r and r0 to do sem-implicit BackwardEuler integration if true
497 if (useSemiImplicit)
498 {
499 for(int i=0;i<nSymTen;i++)
500 {
501 r[i]=r0[i];
502 for(int j=0;j<nSymTen;j++)
503 {
504 dr[i + j*nSymTen]=0.0;
505 }
506 }
507 }
508 //try linesearch if true and f has increased
509 int linesearches=0,max_linesearches=100;
510 while(useLineSearch && (f > 0.99*f_last) && (linesearches < max_linesearches))
511 {
512 std::cout<<"+";
513 double ds = pow(0.5,linesearches+1);
514 Delta_lambda = Delta_lambda_old + ds*Delta_lambda_incr;
515 for(int i=0;i<nSymTen;i++)
516 {
517 Delta_stress[i] = Delta_stress_old[i] + ds*Delta_stress_incr[i];
518 stress[i] = stress_old[i] + ds*Delta_stress_incr[i];
519 }
520 evaluateConstitutiveModel(materialProperties,stress,f,df,r,dr,stress_3);
521 if (useDifferenceJacobian)
522 {
523 differenceJacobian(materialProperties,stress,f,df,r,dr);
524 }
525 //reset r and r0 to do sem-implicit BackwardEuler integration if true
526 if (useSemiImplicit || predictorPhase)
527 {
528 for(int i=0;i<nSymTen;i++)
529 {
530 r[i]=r0[i];
531 for(int j=0;j<nSymTen;j++)
532 {
533 dr[i + j*nSymTen]=0.0;
534 }
535 }
536 }
537 linesearches +=1;
538 }
539 if (linesearches > 0)
540 {
541 if (f > 0.99*f_last)
542 {
543 std::cout<<"f"<<std::endl;
544 its = maxIts-1;//force terminate
545 }
546 else
547 std::cout<<"s"<<std::endl;
548
549 }
550 //calculate residuals and norm(f is already calculated)
551 f_last = f;
552 aNorm=0.0;
553 stressNorm=0.0;
554 elasticStrain(Cinv,Delta_stress,minusDelta_plasticStrain);
555 for (int i=0;i<nSymTen;i++)
556 {
557 a[i] = minusDelta_plasticStrain[i] + Delta_lambda*r[i];
558 aNorm += a[i]*a[i];
559 stressNorm+=stress[i]*stress[i];
560 }
561 aNorm = sqrt(aNorm);
562 stressNorm = sqrt(stressNorm);
563 its++;
564 /* if (its > predictorPhaseIts) */
565 /* predictorPhase = false; */
566 /* // */
567 /* //if we haven't converged then try switchign to semi-implicit integration */
568 /* // */
569 /* bool resetSolution = false; */
570 /* if (its == maxIts && useSemiImplicit==true && useDifferenceJacobian==true */
571 /* && useFullPivoting==true && useDumbNewton == true && useLineSearch == false) */
572 /* { */
573 /* its = 0; */
574 /* useLineSearch=true; */
575 /* resetSolution = true; */
576 /* std::cout<<"linesearch failed"<<std::endl; */
577 /* } */
578 /* if (its == maxIts && useSemiImplicit==true && useDifferenceJacobian==true && useFullPivoting==true && useDumbNewton == false) */
579 /* { */
580 /* its = 0; */
581 /* useDumbNewton=true; */
582 /* resetSolution = true; */
583 /* std::cout<<"dumb Newton failed"<<std::endl; */
584 /* } */
585 /* if (its == maxIts && useSemiImplicit==true && useDifferenceJacobian==true && useFullPivoting==false) */
586 /* { */
587 /* its = 0; */
588 /* useFullPivoting=true; */
589 /* resetSolution = true; */
590 /* std::cout<<"full pivoting failed"<<std::endl; */
591 /* } */
592 /* if (its == maxIts && useSemiImplicit==true && useDifferenceJacobian==false) */
593 /* { */
594 /* its = 0; */
595 /* useDifferenceJacobian=true; */
596 /* resetSolution = true; */
597 /* std::cout<<"semi implicit BE failed"<<std::endl; */
598 /* } */
599 /* if (its == maxIts && useSemiImplicit==false) */
600 /* { */
601 /* std::cout<<"+"<<" f a "<<std::endl; */
602 /* for(int i=0;i<fhist.size();i++) */
603 /* std::cout<<std::setprecision(16)<<fhist[i]<<" "<<ahist[i]<<std::endl; */
604 /* its = 0; */
605 /* useSemiImplicit = true; */
606 /* //useContinuumTangentModulus=false; */
607 /* resetSolution = true; */
608 /* //std::cout<<"fully implicit BE failed, trying semi-implicit"<<std::endl; */
609 /* } */
610 /* if (resetSolution) */
611 /* { */
612 /* //reset to elastic predictor */
613 /* for(int i=0;i<nSymTen;i++) */
614 /* { */
615 /* r[i]=r0[i]; */
616 /* dr[i]=0.0; */
617 /* Delta_lambda = 0.0; */
618 /* Delta_stress[i] = 0.0; */
619 /* minusDelta_plasticStrain[i] = 0.0; */
620 /* plasticStrain[i] = plasticStrain_last[i]; */
621 /* effectiveStrain[i] = strain_last[i] + Delta_strain[i] - plasticStrain_last[i]; */
622 /* a[i] = 0.0; */
623 /* } */
624 /* elasticStress(C,effectiveStrain,stress); */
625 /* evaluateConstitutiveModel(materialProperties,stress,f,df,r,dr,stress_3); */
626 /* if (useDifferenceJacobian) */
627 /* { */
628 /* differenceJacobian(materialProperties,stress,f,df,r,dr); */
629 /* } */
630 /* if (useSemiImplicit) */
631 /* { */
632 /* for(int i=0;i<nSymTen;i++) */
633 /* { */
634 /* r[i]=r0[i]; */
635 /* for(int j=0;j<nSymTen;j++) */
636 /* dr[i + j*nSymTen]=0.0; */
637 /* } */
638 /* } */
639 /* } */
640 }
641 if (its > 0)//non-zero plastic strain increment
642 {
643 //std::cout<<"plastic it "<<its<<" aNorm "<<aNorm<<" f "<<f<<std::endl;
644 //
645 //stress and Delta_lambda have been set, need to set plasticStrain
646 //
647 for(int i=0;i<nSymTen;i++)
648 {
649 plasticStrain[i] = plasticStrain_last[i] - minusDelta_plasticStrain[i];
650 }
651 //now calculate the algorithm tangent modulus (dstress/dstrain)
652 //
653 //if true or the Newton iteration didn't converge try using the continuum modulus
654 if (its == maxIts)
655 {
656 //failed, reset stress to elastic predictor
657 /* std::cout<<"constitutive relation interation failed"<<std::endl; */
658 /* std::cout<<" f a "<<std::endl; */
659 /* for(int i=0;i<fhist.size();i++) */
660 /* std::cout<<std::setprecision(16)<<fhist[i]<<" "<<ahist[i]<<std::endl; */
661 /* std::cout<<"-------------------------------------"<<std::endl; */
662 std::cout<<"Constiutive relation iteration failed "<<f<<'\t'<<stressNorm<<std::endl;
663 //useOuterPicard=true;
664 useContinuumTangentModulus=true;
665 }
666
667 if(useOuterPicard)
668 {
669 for (int i=0; i<nSymTen; i++)
670 for (int j=0; j<nSymTen; j++)
671 dstress[i*nSymTen+j] = C[i+j*nSymTen];
672 }
673 else if(useContinuumTangentModulus)// || its == maxIts)
674 {
675 double Cr[nSymTen],dfC[nSymTen],dfCr=0.0;
676 for (int i=0; i<nSymTen; i++)
677 {
678 Cr[i] = 0.0;
679 dfC[i] = 0.0;
680 }
681 for (int i=0; i<nSymTen; i++)
682 for (int j=0; j<nSymTen; j++)
683 {
684 Cr[i] += C[i+j*nSymTen]*r[j];
685 dfC[j] += df[i]*C[i+j*nSymTen];
686 dfCr += df[i]*C[i+j*nSymTen]*r[j];
687 }
688 if (fabs(dfCr) < 1.0e-8)
689 dfCr = copysign(1.0e-8,dfCr);
690 for (int i=0; i<nSymTen; i++)
691 for (int j=0; j<nSymTen; j++)
692 {
693 dstress[i*nSymTen+j] = C[i+j*nSymTen] - (Cr[i]*dfC[j])/dfCr;//cek hack, storing dstress in C ordering
694 }
695 }
696 else //should be correct for both Backward Euler and Semi-implicit scheme
697 {
698 double Ainv[nSymTen*nSymTen];
699 for(int i=0;i<nSymTen;i++)
700 {
701 for(int j=0;j<nSymTen;j++)
702 {
703 A[i+j*nSymTen] = Cinv[i+j*nSymTen] + Delta_lambda*dr[i+j*nSymTen];//set to A^{-1}
704 Ainv[i+j*nSymTen] = A[i+j*nSymTen];
705 }
706 Aa[i] = a[i];
707 Ar[i] = r[i];
708 dfA[i] = df[i];
709 }
710 dgetrf_(&N,&N,A,&N,IPIV,&INFO);//factor
711 TRANS='T';
712 dgetrs_(&TRANS,&N,&NRHS,A,&N,IPIV,dfA,&N,&INFO);
713 TRANS='N';
714 dgetrs_(&TRANS,&N,&NRHS,A,&N,IPIV,Aa,&N,&INFO);//back substitute to get products with A
715 dgetrs_(&TRANS,&N,&NRHS,A,&N,IPIV,Ar,&N,&INFO);
716 dfAr = 0.0;
717 for(int i=0;i<nSymTen;i++)
718 {
719 dfAr += df[i]*Ar[i];
720 }
721 dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);
722 if (fabs(dfAr) < 1.0e-8)
723 dfAr = copysign(1.0e-8,dfAr);
724 for (int i=0; i<nSymTen; i++)
725 for (int j=0; j<nSymTen; j++)
726 {
727 dstress[i*nSymTen+j] = A[i+j*nSymTen] - (Ar[i]*dfA[j])/dfAr;//cek hack, storing dstress in C ordering
728 }
729 }
730 }
731 else//0 plastic strain i.e. elastic
732 {
733 for (int i=0; i<nSymTen; i++)
734 for (int j=0; j<nSymTen; j++)
735 dstress[i*nSymTen+j] = C[i+j*nSymTen];
736 }
737 //apply pore pressure
738 for (int i=0; i<nSpace; i++)
739 stress[i] -= pore_pressure;
740 }
741
742 inline void exteriorNumericalStressFlux(const double& pore_pressure_ext,
743 const int& isDOFBoundary_u,
744 const int& isDOFBoundary_v,
745 const int& isDOFBoundary_w,
746 const int& isStressFluxBoundary_u,
747 const int& isStressFluxBoundary_v,
748 const int& isStressFluxBoundary_w,
749 const double& penalty,
750 const double& u,
751 const double& v,
752 const double& w,
753 const double& bc_u,
754 const double& bc_v,
755 const double& bc_w,
756 const double& bc_stressFlux_u,
757 const double& bc_stressFlux_v,
758 const double& bc_stressFlux_w,
759 const double* stress,
760 const double* normal,
761 double& stressFlux_u,
762 double& stressFlux_v,
763 double& stressFlux_w)
764 {
765 //note: pore pressure is already in stress but assumed not in stress flux BC
766 if (isDOFBoundary_u == 1)
767 {
768 stressFlux_u = -(stress[sXX]*normal[X] + stress[sXY]*normal[Y] + stress[sXZ]*normal[Z] - penalty*(u - bc_u));
769 }
770 else if(isStressFluxBoundary_u == 1)
771 {
772 stressFlux_u = -bc_stressFlux_u + normal[X]*pore_pressure_ext;
773 }
774 else
775 {
776 stressFlux_u = 0.0;
777 }
778
779 if (isDOFBoundary_v == 1)
780 {
781 stressFlux_v = -(stress[sYX]*normal[X] + stress[sYY]*normal[Y] + stress[sYZ]*normal[Z] - penalty*(v - bc_v));
782 }
783 else if(isStressFluxBoundary_v == 1)
784 {
785 stressFlux_v = -bc_stressFlux_v + normal[Y]*pore_pressure_ext;
786 }
787 else
788 {
789 stressFlux_v = 0.0;
790 }
791
792 if (isDOFBoundary_w == 1)
793 {
794 stressFlux_w = -(stress[sZX]*normal[X] + stress[sZY]*normal[Y] + stress[sZZ]*normal[Z] - penalty*(w - bc_w));
795 }
796 else if(isStressFluxBoundary_w == 1)
797 {
798 stressFlux_w = -bc_stressFlux_w + normal[Z]*pore_pressure_ext;
799 }
800 else
801 {
802 stressFlux_w = 0.0;
803 }
804 }
805
806 inline void exteriorNumericalStressFluxJacobian(const int& isDOFBoundary_u,
807 const int& isDOFBoundary_v,
808 const int& isDOFBoundary_w,
809 const double* normal,
810 const double* dstress,
811 const double& penalty,
812 const double& disp_trial,
813 const double* disp_grad_trial,
814 double& dstressFlux_u_u,
815 double& dstressFlux_u_v,
816 double& dstressFlux_u_w,
817 double& dstressFlux_v_u,
818 double& dstressFlux_v_v,
819 double& dstressFlux_v_w,
820 double& dstressFlux_w_u,
821 double& dstressFlux_w_v,
822 double& dstressFlux_w_w)
823 {
824 //here we use both the symmetry of the stress tensor and the fact that dstress is w.r.t. the strain in Voigt notation to go directly to derivatives w.r.t. displacement DOF
825 if (isDOFBoundary_u == 1)
826 {
827 //stressFlux_u = -(stress[sXX]*normal[X] + stress[sXY]*normal[Y] + stress[sXZ]*normal[Z] - h_penalty*(u - bc_u));
828 dstressFlux_u_u = -(
829 (dstress[sXX*nSymTen+sXX]*disp_grad_trial[X] + dstress[sXX*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sXX*nSymTen+sXZ]*disp_grad_trial[Z])*normal[X]+
830 (dstress[sXY*nSymTen+sXX]*disp_grad_trial[X] + dstress[sXY*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sXY*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Y]+
831 (dstress[sXZ*nSymTen+sXX]*disp_grad_trial[X] + dstress[sXZ*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sXZ*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Z]
832 -
833 penalty*disp_trial);
834 dstressFlux_u_v = -(
835 (dstress[sXX*nSymTen+sYX]*disp_grad_trial[X] + dstress[sXX*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sXX*nSymTen+sYZ]*disp_grad_trial[Z])*normal[X]+
836 (dstress[sXY*nSymTen+sYX]*disp_grad_trial[X] + dstress[sXY*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sXY*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Y]+
837 (dstress[sXZ*nSymTen+sYX]*disp_grad_trial[X] + dstress[sXZ*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sXZ*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Z]);
838 dstressFlux_u_w = -(
839 (dstress[sXX*nSymTen+sZX]*disp_grad_trial[X] + dstress[sXX*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sXX*nSymTen+sZZ]*disp_grad_trial[Z])*normal[X]+
840 (dstress[sXY*nSymTen+sZX]*disp_grad_trial[X] + dstress[sXY*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sXY*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Y]+
841 (dstress[sXZ*nSymTen+sZX]*disp_grad_trial[X] + dstress[sXZ*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sXZ*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Z]);
842 }
843 else
844 {
845 dstressFlux_u_u = 0.0;
846 dstressFlux_u_v = 0.0;
847 dstressFlux_u_w = 0.0;
848 }
849
850 if (isDOFBoundary_v == 1)
851 {
852 //stressFlux_v = -(stress[sYX]*normal[X] + stress[sYY]*normal[Y] + stress[sYZ]*normal[Z] - h_penalty*(v - bc_v));
853 dstressFlux_v_u = -(
854 (dstress[sYX*nSymTen+sXX]*disp_grad_trial[X] + dstress[sYX*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sYX*nSymTen+sXZ]*disp_grad_trial[Z])*normal[X]+
855 (dstress[sYY*nSymTen+sXX]*disp_grad_trial[X] + dstress[sYY*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sYY*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Y]+
856 (dstress[sYZ*nSymTen+sXX]*disp_grad_trial[X] + dstress[sYZ*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sYZ*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Z]);
857 dstressFlux_v_v = -(
858 (dstress[sYX*nSymTen+sYX]*disp_grad_trial[X] + dstress[sYX*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sYX*nSymTen+sYZ]*disp_grad_trial[Z])*normal[X]+
859 (dstress[sYY*nSymTen+sYX]*disp_grad_trial[X] + dstress[sYY*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sYY*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Y]+
860 (dstress[sYZ*nSymTen+sYX]*disp_grad_trial[X] + dstress[sYZ*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sYZ*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Z]
861 -
862 penalty*disp_trial);
863 dstressFlux_v_w = -(
864 (dstress[sYX*nSymTen+sZX]*disp_grad_trial[X] + dstress[sYX*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sYX*nSymTen+sZZ]*disp_grad_trial[Z])*normal[X]+
865 (dstress[sYY*nSymTen+sZX]*disp_grad_trial[X] + dstress[sYY*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sYY*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Y]+
866 (dstress[sYZ*nSymTen+sZX]*disp_grad_trial[X] + dstress[sYZ*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sYZ*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Z]);
867 }
868 else
869 {
870 dstressFlux_v_u = 0.0;
871 dstressFlux_v_v = 0.0;
872 dstressFlux_v_w = 0.0;
873 }
874
875 if (isDOFBoundary_w == 1)
876 {
877 //stressFlux_w = -(stress[sZX]*normal[X] + stress[sZY]*normal[Y] + stress[sZZ]*normal[Z] - h_penalty*(w - bc_w));
878 dstressFlux_w_u = -(
879 (dstress[sZX*nSymTen+sXX]*disp_grad_trial[X] + dstress[sZX*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sZX*nSymTen+sXZ]*disp_grad_trial[Z])*normal[X]+
880 (dstress[sZY*nSymTen+sXX]*disp_grad_trial[X] + dstress[sZY*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sZY*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Y]+
881 (dstress[sZZ*nSymTen+sXX]*disp_grad_trial[X] + dstress[sZZ*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sZZ*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Z]);
882 dstressFlux_w_v = -(
883 (dstress[sZX*nSymTen+sYX]*disp_grad_trial[X] + dstress[sZX*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sZX*nSymTen+sYZ]*disp_grad_trial[Z])*normal[X]+
884 (dstress[sZY*nSymTen+sYX]*disp_grad_trial[X] + dstress[sZY*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sZY*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Y]+
885 (dstress[sZZ*nSymTen+sYX]*disp_grad_trial[X] + dstress[sZZ*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sZZ*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Z]);
886 dstressFlux_w_w = -(
887 (dstress[sZX*nSymTen+sZX]*disp_grad_trial[X] + dstress[sZX*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sZX*nSymTen+sZZ]*disp_grad_trial[Z])*normal[X]+
888 (dstress[sZY*nSymTen+sZX]*disp_grad_trial[X] + dstress[sZY*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sZY*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Y]+
889 (dstress[sZZ*nSymTen+sZX]*disp_grad_trial[X] + dstress[sZZ*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sZZ*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Z]
890 -
891 penalty*disp_trial);
892 }
893 else
894 {
895 dstressFlux_w_u = 0.0;
896 dstressFlux_w_v = 0.0;
897 dstressFlux_w_w = 0.0;
898 }
899 }
900
901 virtual void calculateResidual(arguments_dict& args)
902 {
903 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
904 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
905 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
906 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
907 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
908 xt::pyarray<double>& disp_trial_ref = args.array<double>("disp_trial_ref");
909 xt::pyarray<double>& disp_grad_trial_ref = args.array<double>("disp_grad_trial_ref");
910 xt::pyarray<double>& disp_test_ref = args.array<double>("disp_test_ref");
911 xt::pyarray<double>& disp_grad_test_ref = args.array<double>("disp_grad_test_ref");
912 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
913 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
914 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
915 xt::pyarray<double>& disp_trial_trace_ref = args.array<double>("disp_trial_trace_ref");
916 xt::pyarray<double>& disp_grad_trial_trace_ref = args.array<double>("disp_grad_trial_trace_ref");
917 xt::pyarray<double>& disp_test_trace_ref = args.array<double>("disp_test_trace_ref");
918 xt::pyarray<double>& disp_grad_test_trace_ref = args.array<double>("disp_grad_test_trace_ref");
919 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
920 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
921 xt::pyarray<double>& ebqe_penalty = args.array<double>("ebqe_penalty");
922 int gravityStep = args.scalar<int>("gravityStep");
923 int nElements_global = args.scalar<int>("nElements_global");
924 xt::pyarray<int>& materialTypes = args.array<int>("materialTypes");
925 int nMaterialProperties = args.scalar<int>("nMaterialProperties");
926 xt::pyarray<double>& materialProperties = args.array<double>("materialProperties");
927 double pore_fluid_unit_weight = args.scalar<double>("pore_fluid_unit_weight");
928 xt::pyarray<double>& pore_pressure_head_dof = args.array<double>("pore_pressure_head_dof");
929 xt::pyarray<double>& q_strain = args.array<double>("q_strain");
930 xt::pyarray<double>& q_strain0 = args.array<double>("q_strain0");
931 xt::pyarray<double>& q_strain_last = args.array<double>("q_strain_last");
932 xt::pyarray<double>& q_plasticStrain = args.array<double>("q_plasticStrain");
933 xt::pyarray<double>& q_plasticStrain_last = args.array<double>("q_plasticStrain_last");
934 xt::pyarray<double>& ebqe_strain = args.array<double>("ebqe_strain");
935 xt::pyarray<double>& ebqe_strain0 = args.array<double>("ebqe_strain0");
936 xt::pyarray<double>& ebqe_strain_last = args.array<double>("ebqe_strain_last");
937 xt::pyarray<double>& ebqe_plasticStrain = args.array<double>("ebqe_plasticStrain");
938 xt::pyarray<double>& ebqe_plasticStrain_last = args.array<double>("ebqe_plasticStrain_last");
939 xt::pyarray<int>& disp_l2g = args.array<int>("disp_l2g");
940 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
941 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
942 xt::pyarray<double>& w_dof = args.array<double>("w_dof");
943 xt::pyarray<double>& bodyForce = args.array<double>("bodyForce");
944 int offset_u = args.scalar<int>("offset_u");
945 int offset_v = args.scalar<int>("offset_v");
946 int offset_w = args.scalar<int>("offset_w");
947 int stride_u = args.scalar<int>("stride_u");
948 int stride_v = args.scalar<int>("stride_v");
949 int stride_w = args.scalar<int>("stride_w");
950 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
951 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
952 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
953 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
954 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
955 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
956 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
957 xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
958 xt::pyarray<int>& isStressFluxBoundary_u = args.array<int>("isStressFluxBoundary_u");
959 xt::pyarray<int>& isStressFluxBoundary_v = args.array<int>("isStressFluxBoundary_v");
960 xt::pyarray<int>& isStressFluxBoundary_w = args.array<int>("isStressFluxBoundary_w");
961 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
962 xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
963 xt::pyarray<double>& ebqe_bc_w_ext = args.array<double>("ebqe_bc_w_ext");
964 xt::pyarray<double>& ebqe_bc_stressFlux_u_ext = args.array<double>("ebqe_bc_stressFlux_u_ext");
965 xt::pyarray<double>& ebqe_bc_stressFlux_v_ext = args.array<double>("ebqe_bc_stressFlux_v_ext");
966 xt::pyarray<double>& ebqe_bc_stressFlux_w_ext = args.array<double>("ebqe_bc_stressFlux_w_ext");
967 //
968 //loop over elements to compute volume integrals and load them into element and global residual
969 //
970 //std::cout<<"nElements_global"<<nElements_global<<std::endl;
971 //std::cout<<"nQuadraturePoints_element"<<nQuadraturePoints_element<<std::endl;
972 const int usePicard = 0;
973 for(int eN=0;eN<nElements_global;eN++)
974 {
975 //declare local storage for element residual and initialize
976 double
977 elementResidual_u[nDOF_test_element],
978 elementResidual_v[nDOF_test_element],
979 elementResidual_w[nDOF_test_element];
980 for (int i=0;i<nDOF_test_element;i++)
981 {
982 elementResidual_u[i]=0.0;
983 elementResidual_v[i]=0.0;
984 elementResidual_w[i]=0.0;
985 }//i
986 //
987 //loop over quadrature points and compute integrands
988 //
989 for(int k=0;k<nQuadraturePoints_element;k++)
990 {
991 //compute indices and declare local storage
992 int eN_k = eN*nQuadraturePoints_element+k,
993 eN_k_nSpace=eN_k*nSpace,
994 eN_nDOF_trial_element = eN*nDOF_trial_element,
995 eN_nDOF_mesh_trial_element = eN*nDOF_mesh_trial_element; //index to a vector at a quadrature point
996 double u=0.0,v=0.0,w=0.0,
997 D[nSpace*nSpace],
998 *grad_u(&D[0]),
999 *grad_v(&D[nSpace]),
1000 *grad_w(&D[2*nSpace]),
1001 jac[nSpace*nSpace],
1002 jacDet,
1003 jacInv[nSpace*nSpace],
1004 disp_grad_trial[nDOF_trial_element*nSpace],
1005 disp_test_dV[nDOF_trial_element],
1006 disp_grad_test_dV[nDOF_test_element*nSpace],
1007 dV,x,y,z,
1008 G[nSpace*nSpace],G_dd_G,tr_G,
1009 stress[ck.nSymTen],strain[ck.nSymTen],
1010 Delta_strain[ck.nSymTen],dstress[ck.nSymTen*ck.nSymTen],pore_pressure=0.0;
1011 //get jacobian, etc for mapping reference element
1012 ck.calculateMapping_element(eN,
1013 k,
1014 mesh_dof.data(),
1015 mesh_l2g.data(),
1016 mesh_trial_ref.data(),
1017 mesh_grad_trial_ref.data(),
1018 jac,
1019 jacDet,
1020 jacInv,
1021 x,y,z);
1022 //get the pressure head
1023 ck.mapping.valFromDOF(pore_pressure_head_dof.data(),&mesh_l2g.data()[eN_nDOF_mesh_trial_element],&mesh_trial_ref.data()[k*nDOF_mesh_trial_element],pore_pressure);
1024 //std::cout<<"pore_pressure \t"<<(5.0-pore_pressure)<<"\t z \t"<<z<<std::endl;
1025 pore_pressure *= pore_fluid_unit_weight;
1026 pore_pressure = fmax(pore_pressure,0.0);//cek hack
1027 //std::cout<<"pore_pressure "<<pore_pressure<<std::endl;
1028 //get the physical integration weight
1029 dV = fabs(jacDet)*dV_ref.data()[k];
1030 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1031 //get the trial function gradients
1032 ck.gradTrialFromRef(&disp_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,disp_grad_trial);
1033 //get the solution
1034 ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],u);
1035 ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],v);
1036 ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],w);
1037 //get the solution gradients
1038 ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_u);
1039 ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_v);
1040 ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_w);
1041 //precalculate test function products with integration weights
1042 for (int j=0;j<nDOF_trial_element;j++)
1043 {
1044 disp_test_dV[j] = disp_test_ref.data()[k*nDOF_trial_element+j]*dV;
1045 for (int I=0;I<nSpace;I++)
1046 {
1047 disp_grad_test_dV[j*nSpace+I] = disp_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
1048 }
1049 }
1050 //save displacement at quadrature points for other models to use
1051 //q_displacement[eN_k_nSpace+0]=u;
1052 //q_displacement[eN_k_nSpace+1]=v;
1053 //q_displacement[eN_k_nSpace+2]=w;
1054
1055 calculateStrain(D,Delta_strain);
1056 //std::cout<<"here "<<materialProperties[0]<<'\t'<<materialProperties[1]<<std::endl<<std::flush;
1057
1058 double C[nSymTen*nSymTen],Cinv[nSymTen*nSymTen];
1059 if (gravityStep)
1060 {
1061 elasticModuli(materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],//E,
1062 materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1],//nu,
1063 C,Cinv);
1064 if (gravityStep == 2)
1065 {
1066 double f,df[nSymTen],r[nSymTen],dr[nSymTen*nSymTen],stress_3_init,E=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],nu=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1];
1067 elasticStress(C,&q_strain0.data()[eN_k*nSymTen],stress);
1068 evaluateConstitutiveModel(&materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1069 stress,f,df,r,dr,stress_3_init);//cek hack, this is just to get stress_3_init
1070 const double n_e = 0.6,
1071 P_a = materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+12],
1072 K_i = 500.0;
1073 if (P_a > 0.0)//cek hack, if not a soil then set P_a=0.0
1074 {
1075 E =K_i*P_a*pow(fmax(stress_3_init,P_a)/P_a,n_e);
1076 elasticModuli(E,nu,C,Cinv);
1077 }
1078 }
1079 for (int i=0; i<nSymTen; i++)
1080 {
1081 strain[i] = Delta_strain[i] + q_strain_last.data()[eN_k*nSymTen+i];
1082 for (int j=0; j<nSymTen; j++)
1083 dstress[i*nSymTen+j] = C[i+j*nSymTen];
1084 }
1085 elasticStress(C,strain,stress);
1086 for (int i=0; i<nSpace; i++)
1087 stress[i] -= pore_pressure;
1088 }
1089 else
1090 evaluateCoefficients(usePicard,
1091 pore_pressure,
1092 &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1093 &q_strain0.data()[eN_k*nSymTen],
1094 &q_strain_last.data()[eN_k*nSymTen],
1095 Delta_strain,
1096 &q_plasticStrain_last.data()[eN_k*nSymTen],
1097 &q_plasticStrain.data()[eN_k*nSymTen],
1098 stress,
1099 dstress);
1100 for (int i=0;i<nSymTen;i++)
1101 {
1102 q_strain.data()[eN_k*nSymTen+i] = q_strain_last.data()[eN_k*nSymTen + i] + Delta_strain[i];
1103 }
1104 //
1105 //moving mesh
1106 //
1107 //omit for now
1108 //
1109 //update element residual
1110 //
1111 for(int i=0;i<nDOF_test_element;i++)
1112 {
1113 int i_nSpace=i*nSpace;
1114 elementResidual_u[i] += ck.Stress_u_weak(stress,&disp_grad_test_dV[i_nSpace]) +
1115 ck.Reaction_weak(-bodyForce.data()[eN_k_nSpace+0],disp_test_dV[i]);
1116 elementResidual_v[i] += ck.Stress_v_weak(stress,&disp_grad_test_dV[i_nSpace]) +
1117 ck.Reaction_weak(-bodyForce.data()[eN_k_nSpace+1],disp_test_dV[i]);
1118 elementResidual_w[i] += ck.Stress_w_weak(stress,&disp_grad_test_dV[i_nSpace]) +
1119 ck.Reaction_weak(-bodyForce.data()[eN_k_nSpace+2],disp_test_dV[i]);
1120 // if (k == nQuadraturePoints_element-1)
1121 // {
1122 // std::cout<<"element residual "<<eN<<'\t'<<i<<std::endl;
1123 // std::cout<<elementResidual_u[i]<<std::endl
1124 // <<elementResidual_v[i]<<std::endl
1125 // <<elementResidual_w[i]<<std::endl;
1126 // }
1127 }//i
1128 }
1129 //
1130 //load element into global residual and save element residual
1131 //
1132 for(int i=0;i<nDOF_test_element;i++)
1133 {
1134 int eN_i=eN*nDOF_test_element+i;
1135
1136 globalResidual.data()[offset_u+stride_u*disp_l2g.data()[eN_i]] += elementResidual_u[i];
1137 globalResidual.data()[offset_v+stride_v*disp_l2g.data()[eN_i]] += elementResidual_v[i];
1138 globalResidual.data()[offset_w+stride_w*disp_l2g.data()[eN_i]] += elementResidual_w[i];
1139 }//i
1140 }//elements
1141 //
1142 //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
1143 //
1144 //ebNE is the Exterior element boundary INdex
1145 //ebN is the element boundary INdex
1146 //eN is the element index
1147 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1148 {
1149 int ebN = exteriorElementBoundariesArray.data()[ebNE],
1150 eN = elementBoundaryElementsArray.data()[ebN*2+0],
1151 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
1152 eN_nDOF_trial_element = eN*nDOF_trial_element,
1153 eN_nDOF_mesh_trial_element = eN*nDOF_mesh_trial_element;
1154 double elementResidual_u[nDOF_test_element],
1155 elementResidual_v[nDOF_test_element],
1156 elementResidual_w[nDOF_test_element];
1157 for (int i=0;i<nDOF_test_element;i++)
1158 {
1159 elementResidual_u[i]=0.0;
1160 elementResidual_v[i]=0.0;
1161 elementResidual_w[i]=0.0;
1162 }
1163 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1164 {
1165 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1166 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1167 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1168 double u_ext=0.0,
1169 v_ext=0.0,
1170 w_ext=0.0,
1171 D[nSpace*nSpace],
1172 *grad_u_ext(&D[0]),
1173 *grad_v_ext(&D[nSpace]),
1174 *grad_w_ext(&D[2*nSpace]),
1175 bc_u_ext=0.0,
1176 bc_v_ext=0.0,
1177 bc_w_ext=0.0,
1178 jac_ext[nSpace*nSpace],
1179 jacDet_ext,
1180 jacInv_ext[nSpace*nSpace],
1181 boundaryJac[nSpace*(nSpace-1)],
1182 metricTensor[(nSpace-1)*(nSpace-1)],
1183 metricTensorDetSqrt,
1184 dS,disp_test_dS[nDOF_test_element],
1185 disp_grad_trial_trace[nDOF_trial_element*nSpace],
1186 normal[3],x_ext,y_ext,z_ext,
1187 G[nSpace*nSpace],G_dd_G,tr_G,h_penalty,
1188 Delta_strain[ck.nSymTen],
1189 strain[ck.nSymTen],stress[ck.nSymTen],dstress[ck.nSymTen*ck.nSymTen],
1190 stressFlux_u,stressFlux_v,stressFlux_w,pore_pressure_ext;
1191 //compute information about mapping from reference element to physical element
1192 ck.calculateMapping_elementBoundary(eN,
1193 ebN_local,
1194 kb,
1195 ebN_local_kb,
1196 mesh_dof.data(),
1197 mesh_l2g.data(),
1198 mesh_trial_trace_ref.data(),
1199 mesh_grad_trial_trace_ref.data(),
1200 boundaryJac_ref.data(),
1201 jac_ext,
1202 jacDet_ext,
1203 jacInv_ext,
1204 boundaryJac,
1205 metricTensor,
1206 metricTensorDetSqrt,
1207 normal_ref.data(),
1208 normal,
1209 x_ext,y_ext,z_ext);
1210 //get the pressure head
1211 ck.mapping.valFromDOF(pore_pressure_head_dof.data(),&mesh_l2g.data()[eN_nDOF_mesh_trial_element],&mesh_trial_trace_ref.data()[ebN_local_kb*nDOF_mesh_trial_element],pore_pressure_ext);
1212 pore_pressure_ext *= pore_fluid_unit_weight;
1213 pore_pressure_ext = fmax(pore_pressure_ext,0.0);
1214 dS = metricTensorDetSqrt*dS_ref.data()[kb];
1215 //get the metric tensor
1216 //cek todo use symmetry
1217 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1218 //compute shape and solution information
1219 //shape
1220 ck.gradTrialFromRef(&disp_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,disp_grad_trial_trace);
1221 //solution and gradients
1222 ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
1223 ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
1224 ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext);
1225 ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_u_ext);
1226 ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_v_ext);
1227 ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_w_ext);
1228 //precalculate test function products with integration weights
1229 for (int j=0;j<nDOF_trial_element;j++)
1230 {
1231 disp_test_dS[j] = disp_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1232 }
1233 //
1234 //load the boundary values
1235 //
1236 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
1237 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
1238 bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext;
1239 //
1240 //calculate the pde coefficients using the solution and the boundary values for the solution
1241 //
1242 calculateStrain(D,Delta_strain);
1243 double C[nSymTen*nSymTen],Cinv[nSymTen*nSymTen];
1244 if (gravityStep)
1245 {
1246 elasticModuli(materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],//E,
1247 materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1],//nu,
1248 C,Cinv);
1249 if (gravityStep == 2)
1250 {
1251 double f,df[nSymTen],r[nSymTen],dr[nSymTen*nSymTen],stress_3_init,E=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],nu=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1];
1252 elasticStress(C,&ebqe_strain0.data()[ebNE_kb*nSymTen],stress);
1253 evaluateConstitutiveModel(&materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1254 stress,f,df,r,dr,stress_3_init);//cek hack, this is just to get stress_3_init
1255 const double n_e = 0.6,
1256 P_a = materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+12],
1257 K_i = 500.0;
1258 if (P_a > 0.0)//cek hack, if not a soil then set P_a=0.0
1259 {
1260 E =K_i*P_a*pow(fmax(stress_3_init,P_a)/P_a,n_e);
1261 elasticModuli(E,nu,C,Cinv);
1262 }
1263 }
1264 for (int i=0; i<nSymTen; i++)
1265 {
1266 strain[i] = Delta_strain[i] + ebqe_strain_last.data()[ebNE_kb*nSymTen+i];
1267 for (int j=0; j<nSymTen; j++)
1268 dstress[i*nSymTen+j] = C[i+j*nSymTen];
1269 }
1270 elasticStress(C,strain,stress);
1271 for (int i=0;i<nSpace;i++)
1272 stress[i] -= pore_pressure_ext;
1273 }
1274 else
1275 evaluateCoefficients(usePicard,
1276 pore_pressure_ext,
1277 &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1278 &ebqe_strain0.data()[ebNE_kb*nSymTen],
1279 &ebqe_strain_last.data()[ebNE_kb*nSymTen],
1280 Delta_strain,
1281 &ebqe_plasticStrain_last.data()[ebNE_kb*nSymTen],
1282 &ebqe_plasticStrain.data()[ebNE_kb*nSymTen],
1283 stress,
1284 dstress);
1285 for (int i=0;i<nSymTen;i++)
1286 {
1287 ebqe_strain.data()[ebNE_kb*nSymTen+i] = ebqe_strain_last.data()[ebNE_kb*nSymTen+i] + Delta_strain[i];
1288 }
1289 //
1290 //calculate the numerical fluxes
1291 //
1292 double E=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1293 nu=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1];
1294 h_penalty=(E/(1.0+nu))*(1.0 + (nu/(1.0-2.0*nu)))*ebqe_penalty.data()[ebNE_kb];
1295 exteriorNumericalStressFlux(pore_pressure_ext,
1296 isDOFBoundary_u.data()[ebNE_kb],
1297 isDOFBoundary_v.data()[ebNE_kb],
1298 isDOFBoundary_w.data()[ebNE_kb],
1299 isStressFluxBoundary_u.data()[ebNE_kb],
1300 isStressFluxBoundary_v.data()[ebNE_kb],
1301 isStressFluxBoundary_w.data()[ebNE_kb],
1302 h_penalty,
1303 u_ext,
1304 v_ext,
1305 w_ext,
1306 bc_u_ext,
1307 bc_v_ext,
1308 bc_w_ext,
1309 ebqe_bc_stressFlux_u_ext.data()[ebNE_kb],
1310 ebqe_bc_stressFlux_v_ext.data()[ebNE_kb],
1311 ebqe_bc_stressFlux_w_ext.data()[ebNE_kb],
1312 stress,
1313 normal,
1314 stressFlux_u,
1315 stressFlux_v,
1316 stressFlux_w);
1317 //
1318 //update residuals
1319 //
1320 for (int i=0;i<nDOF_test_element;i++)
1321 {
1322 elementResidual_u[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_u,disp_test_dS[i]);
1323 elementResidual_v[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_v,disp_test_dS[i]);
1324 elementResidual_w[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_w,disp_test_dS[i]);
1325 }//i
1326 }//kb
1327 //
1328 //update the element and global residual storage
1329 //
1330 for (int i=0;i<nDOF_test_element;i++)
1331 {
1332 int eN_i = eN*nDOF_test_element+i;
1333
1334 globalResidual.data()[offset_u+stride_u*disp_l2g.data()[eN_i]]+=elementResidual_u[i];
1335 globalResidual.data()[offset_v+stride_v*disp_l2g.data()[eN_i]]+=elementResidual_v[i];
1336 globalResidual.data()[offset_w+stride_w*disp_l2g.data()[eN_i]]+=elementResidual_w[i];
1337 }//i
1338 }//ebNE
1339 }
1340
1341 virtual void calculateJacobian(arguments_dict& args)
1342 {
1343 int usePicard = args.scalar<int>("usePicard");
1344 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
1345 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
1346 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
1347 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
1348 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
1349 xt::pyarray<double>& disp_trial_ref = args.array<double>("disp_trial_ref");
1350 xt::pyarray<double>& disp_grad_trial_ref = args.array<double>("disp_grad_trial_ref");
1351 xt::pyarray<double>& disp_test_ref = args.array<double>("disp_test_ref");
1352 xt::pyarray<double>& disp_grad_test_ref = args.array<double>("disp_grad_test_ref");
1353 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
1354 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
1355 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
1356 xt::pyarray<double>& disp_trial_trace_ref = args.array<double>("disp_trial_trace_ref");
1357 xt::pyarray<double>& disp_grad_trial_trace_ref = args.array<double>("disp_grad_trial_trace_ref");
1358 xt::pyarray<double>& disp_test_trace_ref = args.array<double>("disp_test_trace_ref");
1359 xt::pyarray<double>& disp_grad_test_trace_ref = args.array<double>("disp_grad_test_trace_ref");
1360 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
1361 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
1362 xt::pyarray<double>& ebqe_penalty = args.array<double>("ebqe_penalty");
1363 int gravityStep = args.scalar<int>("gravityStep");
1364 int nElements_global = args.scalar<int>("nElements_global");
1365 xt::pyarray<int>& materialTypes = args.array<int>("materialTypes");
1366 int nMaterialProperties = args.scalar<int>("nMaterialProperties");
1367 xt::pyarray<double>& materialProperties = args.array<double>("materialProperties");
1368 double pore_fluid_unit_weight = args.scalar<double>("pore_fluid_unit_weight");
1369 xt::pyarray<double>& pore_pressure_head_dof = args.array<double>("pore_pressure_head_dof");
1370 xt::pyarray<double>& q_strain = args.array<double>("q_strain");
1371 xt::pyarray<double>& q_strain0 = args.array<double>("q_strain0");
1372 xt::pyarray<double>& q_strain_last = args.array<double>("q_strain_last");
1373 xt::pyarray<double>& q_plasticStrain = args.array<double>("q_plasticStrain");
1374 xt::pyarray<double>& q_plasticStrain_last = args.array<double>("q_plasticStrain_last");
1375 xt::pyarray<double>& ebqe_strain = args.array<double>("ebqe_strain");
1376 xt::pyarray<double>& ebqe_strain0 = args.array<double>("ebqe_strain0");
1377 xt::pyarray<double>& ebqe_strain_last = args.array<double>("ebqe_strain_last");
1378 xt::pyarray<double>& ebqe_plasticStrain = args.array<double>("ebqe_plasticStrain");
1379 xt::pyarray<double>& ebqe_plasticStrain_last = args.array<double>("ebqe_plasticStrain_last");
1380 xt::pyarray<int>& disp_l2g = args.array<int>("disp_l2g");
1381 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
1382 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
1383 xt::pyarray<double>& w_dof = args.array<double>("w_dof");
1384 xt::pyarray<double>& bodyForce = args.array<double>("bodyForce");
1385 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
1386 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
1387 xt::pyarray<int>& csrRowIndeces_u_v = args.array<int>("csrRowIndeces_u_v");
1388 xt::pyarray<int>& csrColumnOffsets_u_v = args.array<int>("csrColumnOffsets_u_v");
1389 xt::pyarray<int>& csrRowIndeces_u_w = args.array<int>("csrRowIndeces_u_w");
1390 xt::pyarray<int>& csrColumnOffsets_u_w = args.array<int>("csrColumnOffsets_u_w");
1391 xt::pyarray<int>& csrRowIndeces_v_u = args.array<int>("csrRowIndeces_v_u");
1392 xt::pyarray<int>& csrColumnOffsets_v_u = args.array<int>("csrColumnOffsets_v_u");
1393 xt::pyarray<int>& csrRowIndeces_v_v = args.array<int>("csrRowIndeces_v_v");
1394 xt::pyarray<int>& csrColumnOffsets_v_v = args.array<int>("csrColumnOffsets_v_v");
1395 xt::pyarray<int>& csrRowIndeces_v_w = args.array<int>("csrRowIndeces_v_w");
1396 xt::pyarray<int>& csrColumnOffsets_v_w = args.array<int>("csrColumnOffsets_v_w");
1397 xt::pyarray<int>& csrRowIndeces_w_u = args.array<int>("csrRowIndeces_w_u");
1398 xt::pyarray<int>& csrColumnOffsets_w_u = args.array<int>("csrColumnOffsets_w_u");
1399 xt::pyarray<int>& csrRowIndeces_w_v = args.array<int>("csrRowIndeces_w_v");
1400 xt::pyarray<int>& csrColumnOffsets_w_v = args.array<int>("csrColumnOffsets_w_v");
1401 xt::pyarray<int>& csrRowIndeces_w_w = args.array<int>("csrRowIndeces_w_w");
1402 xt::pyarray<int>& csrColumnOffsets_w_w = args.array<int>("csrColumnOffsets_w_w");
1403 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
1404 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
1405 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
1406 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
1407 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
1408 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
1409 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
1410 xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
1411 xt::pyarray<int>& isStressFluxBoundary_u = args.array<int>("isStressFluxBoundary_u");
1412 xt::pyarray<int>& isStressFluxBoundary_v = args.array<int>("isStressFluxBoundary_v");
1413 xt::pyarray<int>& isStressFluxBoundary_w = args.array<int>("isStressFluxBoundary_w");
1414 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
1415 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.array<int>("csrColumnOffsets_eb_u_v");
1416 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.array<int>("csrColumnOffsets_eb_u_w");
1417 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.array<int>("csrColumnOffsets_eb_v_u");
1418 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.array<int>("csrColumnOffsets_eb_v_v");
1419 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.array<int>("csrColumnOffsets_eb_v_w");
1420 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.array<int>("csrColumnOffsets_eb_w_u");
1421 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.array<int>("csrColumnOffsets_eb_w_v");
1422 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.array<int>("csrColumnOffsets_eb_w_w");
1423 CompKernel<nSpace,nDOF_mesh_trial_element,nDOF_trial_element,nDOF_test_element> ck;
1424 const int nSymTen(ck.nSymTen);
1425 //
1426 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
1427 //
1428 for(int eN=0;eN<nElements_global;eN++)
1429 {
1430 double
1431 elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
1432 elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
1433 elementJacobian_u_w[nDOF_test_element][nDOF_trial_element],
1434 elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
1435 elementJacobian_v_v[nDOF_test_element][nDOF_trial_element],
1436 elementJacobian_v_w[nDOF_test_element][nDOF_trial_element],
1437 elementJacobian_w_u[nDOF_test_element][nDOF_trial_element],
1438 elementJacobian_w_v[nDOF_test_element][nDOF_trial_element],
1439 elementJacobian_w_w[nDOF_test_element][nDOF_trial_element];
1440 for (int i=0;i<nDOF_test_element;i++)
1441 for (int j=0;j<nDOF_trial_element;j++)
1442 {
1443 elementJacobian_u_u[i][j]=0.0;
1444 elementJacobian_u_v[i][j]=0.0;
1445 elementJacobian_u_w[i][j]=0.0;
1446 elementJacobian_v_u[i][j]=0.0;
1447 elementJacobian_v_v[i][j]=0.0;
1448 elementJacobian_v_w[i][j]=0.0;
1449 elementJacobian_w_u[i][j]=0.0;
1450 elementJacobian_w_v[i][j]=0.0;
1451 elementJacobian_w_w[i][j]=0.0;
1452 }
1453 for (int k=0;k<nQuadraturePoints_element;k++)
1454 {
1455 const int eN_k=eN*nQuadraturePoints_element+k,
1456 eN_nDOF_trial_element = eN*nDOF_trial_element, //index to a vector at a quadrature point
1457 eN_nDOF_mesh_trial_element = eN*nDOF_mesh_trial_element; //index to a vector at a quadrature point
1458
1459 //declare local storage
1460 double u=0.0,v=0.0,w=0.0,
1461 D[nSpace*nSpace],
1462 *grad_u(&D[0]),
1463 *grad_v(&D[nSpace]),
1464 *grad_w(&D[2*nSpace]),
1465 jac[nSpace*nSpace],
1466 jacDet,
1467 jacInv[nSpace*nSpace],
1468 disp_grad_trial[nDOF_trial_element*nSpace],
1469 dV,
1470 disp_test_dV[nDOF_test_element],
1471 disp_grad_test_dV[nDOF_test_element*nSpace],
1472 x,y,z,
1473 G[nSpace*nSpace],G_dd_G,tr_G,
1474 Delta_strain[nSymTen],strain[nSymTen],
1475 stress[nSymTen],dstress[nSpace*nSpace*nSpace*nSpace],pore_pressure=0.0;
1476 //get jacobian, etc for mapping reference element
1477 ck.calculateMapping_element(eN,
1478 k,
1479 mesh_dof.data(),
1480 mesh_l2g.data(),
1481 mesh_trial_ref.data(),
1482 mesh_grad_trial_ref.data(),
1483 jac,
1484 jacDet,
1485 jacInv,
1486 x,y,z);
1487 //get the pressure head
1488 ck.mapping.valFromDOF(pore_pressure_head_dof.data(),&mesh_l2g.data()[eN_nDOF_mesh_trial_element],&mesh_trial_ref.data()[k*nDOF_mesh_trial_element],pore_pressure);
1489 pore_pressure *= pore_fluid_unit_weight;
1490 pore_pressure = fmax(pore_pressure,0.0);//cek hack
1491 //get the physical integration weight
1492 dV = fabs(jacDet)*dV_ref.data()[k];
1493 ck.calculateG(jacInv,G,G_dd_G,tr_G);
1494 //get the trial function gradients
1495 ck.gradTrialFromRef(&disp_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,disp_grad_trial);
1496 //get the solution
1497 ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],u);
1498 ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],v);
1499 ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],w);
1500 //get the solution gradients
1501 ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_u);
1502 ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_v);
1503 ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_w);
1504 //precalculate test function products with integration weights
1505 for (int j=0;j<nDOF_trial_element;j++)
1506 {
1507 disp_test_dV[j] = disp_test_ref.data()[k*nDOF_trial_element+j]*dV;
1508 for (int I=0;I<nSpace;I++)
1509 {
1510 disp_grad_test_dV[j*nSpace+I] = disp_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin}
1511 }
1512 }
1513 calculateStrain(D,Delta_strain);
1514 if (gravityStep)
1515 {
1516 double C[nSymTen*nSymTen],Cinv[nSymTen*nSymTen];
1517 elasticModuli(materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],//E,
1518 materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1],//nu,
1519 C,Cinv);
1520 if (gravityStep == 2)
1521 {
1522 double f,df[nSymTen],r[nSymTen],dr[nSymTen*nSymTen],stress_3_init,E=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],nu=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1];
1523 elasticStress(C,&q_strain0.data()[eN_k*nSymTen],stress);
1524 evaluateConstitutiveModel(&materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1525 stress,f,df,r,dr,stress_3_init);//cek hack, this is just to get stress_3_init
1526 const double n_e = 0.6,
1527 P_a = materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+12],
1528 K_i = 500.0;
1529 if (P_a > 0.0)//cek hack, if not a soil then set P_a=0.0
1530 {
1531 E =K_i*P_a*pow(fmax(stress_3_init,P_a)/P_a,n_e);
1532 elasticModuli(E,nu,C,Cinv);
1533 }
1534 }
1535 for (int i=0; i<nSymTen; i++)
1536 {
1537 strain[i] = Delta_strain[i] + q_strain_last.data()[eN_k*nSymTen+i];
1538 for (int j=0; j<nSymTen; j++)
1539 dstress[i*nSymTen+j] = C[i+j*nSymTen];
1540 }
1541 elasticStress(C,strain,stress);
1542 for (int i=0; i<nSpace; i++)
1543 stress[i] -= pore_pressure;
1544 }
1545 else
1546 evaluateCoefficients(usePicard,
1547 pore_pressure,
1548 &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1549 &q_strain0.data()[eN_k*nSymTen],
1550 &q_strain_last.data()[eN_k*nSymTen],
1551 Delta_strain,
1552 &q_plasticStrain_last.data()[eN_k*nSymTen],
1553 &q_plasticStrain.data()[eN_k*nSymTen],
1554 stress,
1555 dstress);
1556 //
1557 //moving mesh
1558 //
1559 //omit for now
1560 //
1561 for(int i=0;i<nDOF_test_element;i++)
1562 {
1563 int i_nSpace = i*nSpace;
1564 for(int j=0;j<nDOF_trial_element;j++)
1565 {
1566 int j_nSpace = j*nSpace;
1567
1568 elementJacobian_u_u[i][j] += ck.StressJacobian_u_u_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1569 elementJacobian_u_v[i][j] += ck.StressJacobian_u_v_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1570 elementJacobian_u_w[i][j] += ck.StressJacobian_u_w_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1571
1572 elementJacobian_v_u[i][j] += ck.StressJacobian_v_u_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1573 elementJacobian_v_v[i][j] += ck.StressJacobian_v_v_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1574 elementJacobian_v_w[i][j] += ck.StressJacobian_v_w_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1575
1576 elementJacobian_w_u[i][j] += ck.StressJacobian_w_u_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1577 elementJacobian_w_v[i][j] += ck.StressJacobian_w_v_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1578 elementJacobian_w_w[i][j] += ck.StressJacobian_w_w_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
1579 }//j
1580 }//i
1581 }//k
1582 //
1583 //load into element Jacobian into global Jacobian
1584 //
1585 for (int i=0;i<nDOF_test_element;i++)
1586 {
1587 int eN_i = eN*nDOF_test_element+i;
1588 for (int j=0;j<nDOF_trial_element;j++)
1589 {
1590 int eN_i_j = eN_i*nDOF_trial_element+j;
1591 // std::cout<<"i "<<i<<"j "<<j<<std::endl
1592 // <<elementJacobian_u_u[i][j]<<'\t'
1593 // <<elementJacobian_u_v[i][j]<<'\t'
1594 // <<elementJacobian_u_w[i][j]<<std::endl
1595 // <<elementJacobian_v_u[i][j]<<'\t'
1596 // <<elementJacobian_v_v[i][j]<<'\t'
1597 // <<elementJacobian_v_w[i][j]<<std::endl
1598 // <<elementJacobian_w_u[i][j]<<'\t'
1599 // <<elementJacobian_w_v[i][j]<<'\t'
1600 // <<elementJacobian_w_w[i][j]<<std::endl;
1601
1602 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_u_u[eN_i_j]] += elementJacobian_u_u[i][j];
1603 globalJacobian.data()[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_u_v[eN_i_j]] += elementJacobian_u_v[i][j];
1604 globalJacobian.data()[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_u_w[eN_i_j]] += elementJacobian_u_w[i][j];
1605
1606 globalJacobian.data()[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_v_u[eN_i_j]] += elementJacobian_v_u[i][j];
1607 globalJacobian.data()[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_v_v[eN_i_j]] += elementJacobian_v_v[i][j];
1608 globalJacobian.data()[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_v_w[eN_i_j]] += elementJacobian_v_w[i][j];
1609
1610 globalJacobian.data()[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_w_u[eN_i_j]] += elementJacobian_w_u[i][j];
1611 globalJacobian.data()[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_w_v[eN_i_j]] += elementJacobian_w_v[i][j];
1612 globalJacobian.data()[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_w_w[eN_i_j]] += elementJacobian_w_w[i][j];
1613 }//j
1614 }//i
1615 }//elements
1616 //
1617 //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
1618 //
1619 for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
1620 {
1621 int ebN = exteriorElementBoundariesArray.data()[ebNE],
1622 eN = elementBoundaryElementsArray.data()[ebN*2+0],
1623 eN_nDOF_trial_element = eN*nDOF_trial_element,
1624 eN_nDOF_mesh_trial_element = eN*nDOF_mesh_trial_element,
1625 ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
1626 for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
1627 {
1628 int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
1629 ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
1630 ebN_local_kb_nSpace = ebN_local_kb*nSpace;
1631
1632 double
1633 u_ext=0.0,
1634 v_ext=0.0,
1635 w_ext=0.0,
1636 D_ext[nSpace*nSpace],
1637 *grad_u_ext(&D_ext[0]),
1638 *grad_v_ext(&D_ext[nSpace]),
1639 *grad_w_ext(&D_ext[2*nSpace]),
1640 fluxJacobian_u_u[nDOF_trial_element],
1641 fluxJacobian_u_v[nDOF_trial_element],
1642 fluxJacobian_u_w[nDOF_trial_element],
1643 fluxJacobian_v_u[nDOF_trial_element],
1644 fluxJacobian_v_v[nDOF_trial_element],
1645 fluxJacobian_v_w[nDOF_trial_element],
1646 fluxJacobian_w_u[nDOF_trial_element],
1647 fluxJacobian_w_v[nDOF_trial_element],
1648 fluxJacobian_w_w[nDOF_trial_element],
1649 jac_ext[nSpace*nSpace],
1650 jacDet_ext,
1651 jacInv_ext[nSpace*nSpace],
1652 boundaryJac[nSpace*(nSpace-1)],
1653 metricTensor[(nSpace-1)*(nSpace-1)],
1654 metricTensorDetSqrt,
1655 disp_grad_trial_trace[nDOF_trial_element*nSpace],
1656 dS,
1657 disp_test_dS[nDOF_test_element],
1658 normal[3],
1659 x_ext,y_ext,z_ext,
1660 G[nSpace*nSpace],G_dd_G,tr_G,h_penalty,
1661 Delta_strain[nSymTen],
1662 strain[nSymTen],
1663 stress[nSymTen],
1664 dstress[nSymTen*nSymTen],pore_pressure_ext;
1665 ck.calculateMapping_elementBoundary(eN,
1666 ebN_local,
1667 kb,
1668 ebN_local_kb,
1669 mesh_dof.data(),
1670 mesh_l2g.data(),
1671 mesh_trial_trace_ref.data(),
1672 mesh_grad_trial_trace_ref.data(),
1673 boundaryJac_ref.data(),
1674 jac_ext,
1675 jacDet_ext,
1676 jacInv_ext,
1677 boundaryJac,
1678 metricTensor,
1679 metricTensorDetSqrt,
1680 normal_ref.data(),
1681 normal,
1682 x_ext,y_ext,z_ext);
1683 //get the pressure head
1684 ck.mapping.valFromDOF(pore_pressure_head_dof.data(),&mesh_l2g.data()[eN_nDOF_mesh_trial_element],&mesh_trial_trace_ref.data()[ebN_local_kb*nDOF_mesh_trial_element],pore_pressure_ext);
1685 pore_pressure_ext *= pore_fluid_unit_weight;
1686 pore_pressure_ext = fmax(pore_pressure_ext,0.0);
1687 dS = metricTensorDetSqrt*dS_ref.data()[kb];
1688 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
1689 //compute shape and solution information
1690 //shape
1691 ck.gradTrialFromRef(&disp_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,disp_grad_trial_trace);
1692 //solution and gradients
1693 ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
1694 ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
1695 ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext);
1696 ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_u_ext);
1697 ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_v_ext);
1698 ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_w_ext);
1699 //precalculate test function products with integration weights
1700 for (int j=0;j<nDOF_trial_element;j++)
1701 {
1702 disp_test_dS[j] = disp_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
1703 }
1704 //
1705 //calculate the internal and external trace of the pde coefficients
1706 //
1707 calculateStrain(D_ext,Delta_strain);
1708 if (gravityStep)
1709 {
1710 double C[nSymTen*nSymTen],Cinv[nSymTen*nSymTen];
1711 elasticModuli(materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],//E,
1712 materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1],//nu,
1713 C,Cinv);
1714 if (gravityStep == 2)
1715 {
1716 double f,df[nSymTen],r[nSymTen],dr[nSymTen*nSymTen],stress_3_init,E=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],nu=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1];
1717 elasticStress(C,&ebqe_strain0.data()[ebNE_kb*nSymTen],stress);
1718 evaluateConstitutiveModel(&materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1719 stress,f,df,r,dr,stress_3_init);//cek hack, this is just to get stress_3_init
1720 const double n_e = 0.6,
1721 P_a = materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+12],
1722 K_i = 500.0;
1723 if (P_a > 0.0)//cek hack, if not a soil then set P_a=0.0
1724 {
1725 E =K_i*P_a*pow(fmax(stress_3_init,P_a)/P_a,n_e);
1726 elasticModuli(E,nu,C,Cinv);
1727 }
1728 }
1729 for (int i=0; i<nSymTen; i++)
1730 {
1731 strain[i] = Delta_strain[i] + ebqe_strain_last.data()[ebNE_kb*nSymTen+i];
1732 for (int j=0; j<nSymTen; j++)
1733 dstress[i*nSymTen+j] = C[i+j*nSymTen];
1734 }
1735 elasticStress(C,strain,stress);
1736 for (int i=0;i<nSpace;i++)
1737 stress[i] -= pore_pressure_ext;
1738 }
1739 else
1740 evaluateCoefficients(usePicard,
1741 pore_pressure_ext,
1742 &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1743 &ebqe_strain0.data()[ebNE_kb*nSymTen],
1744 &ebqe_strain_last.data()[ebNE_kb*nSymTen],
1745 Delta_strain,
1746 &ebqe_plasticStrain_last.data()[ebNE_kb*nSymTen],
1747 &ebqe_plasticStrain.data()[ebNE_kb*nSymTen],
1748 stress,
1749 dstress);
1750 //
1751 //calculate the flux jacobian
1752 //
1753 double E=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
1754 nu=materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties+1];
1755 h_penalty=(E/(1.0+nu))*(1.0 + (nu/(1.0-2.0*nu)))*ebqe_penalty.data()[ebNE_kb];
1756 for (int j=0;j<nDOF_trial_element;j++)
1757 {
1758 int j_nSpace = j*nSpace;
1759
1760 exteriorNumericalStressFluxJacobian(isDOFBoundary_u.data()[ebNE_kb],
1761 isDOFBoundary_v.data()[ebNE_kb],
1762 isDOFBoundary_w.data()[ebNE_kb],
1763 normal,
1764 dstress,
1765 h_penalty,
1766 disp_trial_trace_ref.data()[ebN_local_kb*nDOF_trial_element+j],
1767 &disp_grad_trial_trace[j_nSpace],
1768 fluxJacobian_u_u[j],
1769 fluxJacobian_u_v[j],
1770 fluxJacobian_u_w[j],
1771 fluxJacobian_v_u[j],
1772 fluxJacobian_v_v[j],
1773 fluxJacobian_v_w[j],
1774 fluxJacobian_w_u[j],
1775 fluxJacobian_w_v[j],
1776 fluxJacobian_w_w[j]);
1777 }//j
1778 //
1779 //update the global Jacobian from the flux Jacobian
1780 //
1781 for (int i=0;i<nDOF_test_element;i++)
1782 {
1783 int eN_i = eN*nDOF_test_element+i;
1784 for (int j=0;j<nDOF_trial_element;j++)
1785 {
1786 int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
1787
1788 globalJacobian.data()[csrRowIndeces_u_u[eN_i] + csrColumnOffsets_eb_u_u[ebN_i_j]] += fluxJacobian_u_u[j]*disp_test_dS[i];
1789 globalJacobian.data()[csrRowIndeces_u_v[eN_i] + csrColumnOffsets_eb_u_v[ebN_i_j]] += fluxJacobian_u_v[j]*disp_test_dS[i];
1790 globalJacobian.data()[csrRowIndeces_u_w[eN_i] + csrColumnOffsets_eb_u_w[ebN_i_j]] += fluxJacobian_u_w[j]*disp_test_dS[i];
1791
1792 globalJacobian.data()[csrRowIndeces_v_u[eN_i] + csrColumnOffsets_eb_v_u[ebN_i_j]] += fluxJacobian_v_u[j]*disp_test_dS[i];
1793 globalJacobian.data()[csrRowIndeces_v_v[eN_i] + csrColumnOffsets_eb_v_v[ebN_i_j]] += fluxJacobian_v_v[j]*disp_test_dS[i];
1794 globalJacobian.data()[csrRowIndeces_v_w[eN_i] + csrColumnOffsets_eb_v_w[ebN_i_j]] += fluxJacobian_v_w[j]*disp_test_dS[i];
1795
1796 globalJacobian.data()[csrRowIndeces_w_u[eN_i] + csrColumnOffsets_eb_w_u[ebN_i_j]] += fluxJacobian_w_u[j]*disp_test_dS[i];
1797 globalJacobian.data()[csrRowIndeces_w_v[eN_i] + csrColumnOffsets_eb_w_v[ebN_i_j]] += fluxJacobian_w_v[j]*disp_test_dS[i];
1798 globalJacobian.data()[csrRowIndeces_w_w[eN_i] + csrColumnOffsets_eb_w_w[ebN_i_j]] += fluxJacobian_w_w[j]*disp_test_dS[i];
1799 }//j
1800 }//i
1801 }//kb
1802 }//ebNE
1803 }//computeJacobian
1804 };//ElastoPlastic
1805
1806 inline ElastoPlastic_base* newElastoPlastic(int nSpaceIn,
1807 int nQuadraturePoints_elementIn,
1808 int nDOF_mesh_trial_elementIn,
1809 int nDOF_trial_elementIn,
1810 int nDOF_test_elementIn,
1811 int nQuadraturePoints_elementBoundaryIn,
1812 int CompKernelFlag)
1813 {
1814 return proteus::chooseAndAllocateDiscretization<ElastoPlastic_base,ElastoPlastic,CompKernel>(nSpaceIn,
1815 nQuadraturePoints_elementIn,
1816 nDOF_mesh_trial_elementIn,
1817 nDOF_trial_elementIn,
1818 nDOF_test_elementIn,
1819 nQuadraturePoints_elementBoundaryIn,
1820 CompKernelFlag);
1821 }
1822}//proteus
1823
1824#endif
r
Double r
Definition Headers.h:83
B
Double * B
Definition Headers.h:41
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
R
Double R
Definition Headers.h:82
proteus::ElastoPlastic_base
Definition ElastoPlastic.h:30
proteus::ElastoPlastic_base::~ElastoPlastic_base
virtual ~ElastoPlastic_base()
Definition ElastoPlastic.h:32
proteus::ElastoPlastic_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::ElastoPlastic_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::ElastoPlastic::calculateResidual
virtual void calculateResidual(arguments_dict &args)
Definition ElastoPlastic.h:901
proteus::ElastoPlastic::evaluateCoefficients
void evaluateCoefficients(int usePicard, const double pore_pressure, const double *materialProperties, const double *strain0, const double *strain_last, const double *Delta_strain, const double *plasticStrain_last, double *plasticStrain, double *stress, double *dstress)
Definition ElastoPlastic.h:267
proteus::ElastoPlastic::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition ElastoPlastic.h:49
proteus::ElastoPlastic::ei
EIndex< nSpace > ei
Definition ElastoPlastic.h:51
proteus::ElastoPlastic::elasticStrain
void elasticStrain(const double *Cinv, const double *stress, double *strain)
Definition ElastoPlastic.h:108
proteus::ElastoPlastic::exteriorNumericalStressFlux
void exteriorNumericalStressFlux(const double &pore_pressure_ext, const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isDOFBoundary_w, const int &isStressFluxBoundary_u, const int &isStressFluxBoundary_v, const int &isStressFluxBoundary_w, const double &penalty, const double &u, const double &v, const double &w, const double &bc_u, const double &bc_v, const double &bc_w, const double &bc_stressFlux_u, const double &bc_stressFlux_v, const double &bc_stressFlux_w, const double *stress, const double *normal, double &stressFlux_u, double &stressFlux_v, double &stressFlux_w)
Definition ElastoPlastic.h:742
proteus::ElastoPlastic::elasticModuli
void elasticModuli(const double &E, const double &nu, double *C, double *Cinv)
Definition ElastoPlastic.h:118
proteus::ElastoPlastic::differenceJacobian
void differenceJacobian(const double *materialProperties, const double *stress, const double &f, double *df, const double *r, double *dr)
Definition ElastoPlastic.h:243
proteus::ElastoPlastic::evaluateConstitutiveModel
void evaluateConstitutiveModel(const double *materialProperties, const double *stress, double &f, double *df, double *r, double *dr, double &stress_3)
Definition ElastoPlastic.h:170
proteus::ElastoPlastic::elasticStress
void elasticStress(const double *C, const double *strain, double *stress)
Definition ElastoPlastic.h:98
proteus::ElastoPlastic::exteriorNumericalStressFluxJacobian
void exteriorNumericalStressFluxJacobian(const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isDOFBoundary_w, const double *normal, const double *dstress, const double &penalty, const double &disp_trial, const double *disp_grad_trial, double &dstressFlux_u_u, double &dstressFlux_u_v, double &dstressFlux_u_w, double &dstressFlux_v_u, double &dstressFlux_v_v, double &dstressFlux_v_w, double &dstressFlux_w_u, double &dstressFlux_w_v, double &dstressFlux_w_w)
Definition ElastoPlastic.h:806
proteus::ElastoPlastic::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)
Definition ElastoPlastic.h:1341
proteus::ElastoPlastic::calculateStrain
void calculateStrain(double *D, double *strain)
Definition ElastoPlastic.h:87
df
double df(double C, double b, double a, int q, int r)
Definition analyticalSolutions.c:2209
D
#define D(x, y)
Definition jf.h:9
w
#define w(x)
Definition jf.h:22
mcdp
void mcdp(const double &phi, const double &psi, const double &c, const double *stress, double &f, double *df, double &g, double *r, double *dr, double &stress_3)
Definition mcdp.h:2
proteus
Definition ADR.h:19
proteus::phi
double phi(const double &g, const double &h, const double &hL, const double &hR, const double &uL, const double &uR)
Definition SW2DCV.h:62
proteus::newElastoPlastic
ElastoPlastic_base * newElastoPlastic(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition ElastoPlastic.h:1806
proteus::chooseAndAllocateDiscretization
Model_Base * chooseAndAllocateDiscretization(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:48
proteus::f
double f(const double &g, const double &h, const double &hZ)
Definition SW2DCV.h:58
proteus_blas.h
dgetrf_
int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info)
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
dgetri_
int dgetri_(int *N, double *A, int *LDA, int *IPIV, double *WORK, int *LWORK, int *INFO)
dgetrs_
int dgetrs_(char *trans, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info)
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
proteus::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0