proteus  1.8.1
C/C++/Fortran libraries
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 
25 namespace py = pybind11;
26 
27 namespace 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  register 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  register 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  register 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  register 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  register 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  register 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  register 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  register 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  register 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  register 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  register 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  register int i_nSpace = i*nSpace;
1564  for(int j=0;j<nDOF_trial_element;j++)
1565  {
1566  register 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  register int eN_i = eN*nDOF_test_element+i;
1588  for (int j=0;j<nDOF_trial_element;j++)
1589  {
1590  register 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  register 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  register 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  register 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  register 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  register int eN_i = eN*nDOF_test_element+i;
1784  for (int j=0;j<nDOF_trial_element;j++)
1785  {
1786  register 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
proteus::ElastoPlastic::YX
const int YX
Definition: ElastoPlastic.h:55
proteus::ElastoPlastic::sXX
const int sXX
Definition: ElastoPlastic.h:57
proteus::ElastoPlastic::sYY
const int sYY
Definition: ElastoPlastic.h:58
druckerPrager.h
w
#define w(x)
Definition: jf.h:22
B
Double * B
Definition: Headers.h:41
proteus::ElastoPlastic::sXZ
const int sXZ
Definition: ElastoPlastic.h:57
D
#define D(x, y)
Definition: jf.h:9
proteus::ElastoPlastic::YY
const int YY
Definition: ElastoPlastic.h:55
dgesc2_
int dgesc2_(int *n, double *a, int *lda, double *rhs, int *ipiv, int *jpiv, double *scale)
proteus::ElastoPlastic::X
const int X
Definition: ElastoPlastic.h:53
proteus::ElastoPlastic::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)
Definition: ElastoPlastic.h:1341
mohrCoulomb.h
proteus::ElastoPlastic::ElastoPlastic
ElastoPlastic()
Definition: ElastoPlastic.h:67
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::ElastoPlastic::elasticModuli
void elasticModuli(const double &E, const double &nu, double *C, double *Cinv)
Definition: ElastoPlastic.h:118
proteus::ElastoPlastic::sYZ
const int sYZ
Definition: ElastoPlastic.h:58
proteus::ElastoPlastic::nSymTen
const int nSymTen
Definition: ElastoPlastic.h:60
df
double df(double C, double b, double a, int q, int r)
Definition: analyticalSolutions.c:2209
proteus::ElastoPlastic::XX
const int XX
Definition: ElastoPlastic.h:54
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::ElastoPlastic_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
dgetrf_
int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info)
dgetri_
int dgetri_(int *N, double *A, int *LDA, int *IPIV, double *WORK, int *LWORK, int *INFO)
R
Double R
Definition: Headers.h:82
proteus::ElastoPlastic::sZZ
const int sZZ
Definition: ElastoPlastic.h:59
proteus::ElastoPlastic::sXY
const int sXY
Definition: ElastoPlastic.h:57
proteus::ElastoPlastic::ZY
const int ZY
Definition: ElastoPlastic.h:56
proteus::ElastoPlastic::calculateStrain
void calculateStrain(double *D, double *strain)
Definition: ElastoPlastic.h:87
proteus::ElastoPlastic::XHX
const int XHX
Definition: ElastoPlastic.h:61
proteus::ElastoPlastic::YHY
const int YHY
Definition: ElastoPlastic.h:62
proteus::ElastoPlastic::XZ
const int XZ
Definition: ElastoPlastic.h:54
v
Double v
Definition: Headers.h:95
proteus::ElastoPlastic::ZX
const int ZX
Definition: ElastoPlastic.h:56
proteus::ElastoPlastic::sZX
const int sZX
Definition: ElastoPlastic.h:59
proteus::ElastoPlastic
Definition: ElastoPlastic.h:45
proteus::ElastoPlastic::Y
const int Y
Definition: ElastoPlastic.h:53
mohrCoulomb2.h
z
Double * z
Definition: Headers.h:49
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::ElastoPlastic::HXHX
const int HXHX
Definition: ElastoPlastic.h:64
u
Double u
Definition: Headers.h:89
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
dgetc2_
int dgetc2_(int *n, double *a, int *lda, int *ipiv, int *jpiv, int *info)
proteus::ElastoPlastic::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: ElastoPlastic.h:49
tresca.h
proteus::ElastoPlastic_base::~ElastoPlastic_base
virtual ~ElastoPlastic_base()
Definition: ElastoPlastic.h:32
mcdp.h
proteus::ElastoPlastic::HYHY
const int HYHY
Definition: ElastoPlastic.h:65
proteus::ElastoPlastic::sZY
const int sZY
Definition: ElastoPlastic.h:59
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::HYHX
const int HYHX
Definition: ElastoPlastic.h:65
proteus::ElastoPlastic::YHX
const int YHX
Definition: ElastoPlastic.h:62
CompKernel
Definition: CompKernel.h:1018
proteus::ElastoPlastic::ZHX
const int ZHX
Definition: ElastoPlastic.h:63
proteus::ElastoPlastic::HXHY
const int HXHY
Definition: ElastoPlastic.h:64
proteus::ElastoPlastic::Z
const int Z
Definition: ElastoPlastic.h:53
proteus::ElastoPlastic::calculateResidual
virtual void calculateResidual(arguments_dict &args)
Definition: ElastoPlastic.h:901
proteus::ElastoPlastic_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
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
Definition: ADR.h:17
dgetrs_
int dgetrs_(char *trans, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info)
proteus::ElastoPlastic_base
Definition: ElastoPlastic.h:30
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::f
double f(const double &g, const double &h, const double &hZ)
Definition: SW2DCV.h:58
proteus::ElastoPlastic::YZ
const int YZ
Definition: ElastoPlastic.h:55
r
Double r
Definition: Headers.h:83
proteus::ElastoPlastic::ei
EIndex< nSpace > ei
Definition: ElastoPlastic.h:51
proteus::ElastoPlastic::ZHY
const int ZHY
Definition: ElastoPlastic.h:63
proteus::ElastoPlastic::ck
CompKernelType ck
Definition: ElastoPlastic.h:47
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::ElastoPlastic::XHY
const int XHY
Definition: ElastoPlastic.h:61
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::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::elasticStress
void elasticStress(const double *C, const double *strain, double *stress)
Definition: ElastoPlastic.h:98
vonMises.h
proteus::ElastoPlastic::sYX
const int sYX
Definition: ElastoPlastic.h:58
proteus::ElastoPlastic::ZZ
const int ZZ
Definition: ElastoPlastic.h:56
proteus::ElastoPlastic::XY
const int XY
Definition: ElastoPlastic.h:54
EIndex< nSpace >
proteus_blas.h
proteus::ElastoPlastic::elasticStrain
void elasticStrain(const double *Cinv, const double *stress, double *strain)
Definition: ElastoPlastic.h:108
Generated on Fri Jul 1 2022 11:17:29 for proteus by doxygen 1.8.20