proteus  1.8.1
C/C++/Fortran libraries
MoveMesh2D.h
Go to the documentation of this file.
1 #ifndef MOVEMESH2D_H
2 #define MOVEMESH2D_H
3 #include <cmath>
4 #include <iostream>
5 #include "CompKernel.h"
6 #include "ModelFactory.h"
7 #include "ArgumentsDict.h"
8 #include "xtensor-python/pyarray.hpp"
9 
10 namespace py = pybind11;
11 
12 namespace proteus
13 {
14  class MoveMesh2D_base
15  {
16  public:
17  virtual ~MoveMesh2D_base(){}
18  virtual void calculateResidual(arguments_dict& args)=0;
19  virtual void calculateJacobian(arguments_dict& args)=0;
20  };
21 
22  template<class CompKernelType,
23  int nSpace,
24  int nQuadraturePoints_element,
25  int nDOF_mesh_trial_element,
26  int nDOF_trial_element,
27  int nDOF_test_element,
28  int nQuadraturePoints_elementBoundary>
29  class MoveMesh2D : public MoveMesh2D_base
30  {
31  public:
32  CompKernelType ck;
33 
34  const int nDOF_test_X_trial_element;
35 
36  EIndex<nSpace> ei;
37 
38  const int X,Y,/* Z, */
39  XX,XY,/* XZ, */
40  YX,YY,/* YZ, */
41  /* ZX,ZY,ZZ, */
42  sXX,sXY,/* sXZ, */
43  sYX,sYY,/* sYZ, */
44  /* sZX,sZY,sZZ, */
45  nSymTen,
46  XHX,/* XHY, */
47  YHX,/* YHY, */
48  /* ZHX,ZHY, */
49  HXHX/*, HXHY, */
50  /* HYHX,HYHY */;
51 
52  MoveMesh2D():
53  ck(),
54  nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
55  X(ei.X),
56  Y(ei.Y),
57  /* Z(ei.Z), */
58  XX(ei.XX),XY(ei.XY),/* XZ(ei.XZ), */
59  YX(ei.YX),YY(ei.YY),/* YZ(ei.YZ), */
60  /* ZX(ei.ZX),ZY(ei.ZY),ZZ(ei.ZZ), */
61  sXX(ei.sXX),sXY(ei.sXY),/* sXZ(ei.sXZ), */
62  sYX(ei.sYX),sYY(ei.sYY),/* sYZ(ei.sYZ), */
63  /* sZX(ei.sZX),sZY(ei.sZY),sZZ(ei.sZZ), */
64  nSymTen(ei.nSymTen),
65  XHX(ei.XHX),/* XHY(ei.XHY), */
66  YHX(ei.YHX),/* YHY(ei.YHY), */
67  /* ZHX(ei.ZHX),ZHY(ei.ZHY), */
68  HXHX(ei.HXHX)/* ,HXHY(ei.HXHY), */
69  /* HYHX(ei.HYHX),HYHY(ei.HYHY) */
70  {}
71 
72  inline void calculateStrain(double* D, double* strain)
73  {
74  //Voigt notation from Belytschko, Liu, Moran
75  strain[sXX] = D[XX];//du/dx
76  strain[sYY] = D[YY];//dv/dy
77  /* strain[sZZ] = D[ZZ];//dw/dz */
78  /* strain[sYZ] = D[YZ]+D[ZY];//(dv/dz + dz/dy) */
79  /* strain[sXZ] = D[XZ]+D[ZX];//(du/dz + dw/dx) */
80  strain[sXY] = D[XY]+D[YX];//(du/dy + dv/dx)
81  }
82 
83  inline void evaluateCoefficients(double det_J,
84  const double* materialProperties,
85  double* strain,
86  double* stress,
87  double* dstress)
88  {
89  //cek hack/todo need to set E based on reference configuration
90  const double strainTrace=(strain[sXX]+strain[sYY]/* +strain[sZZ] */),
91  E=materialProperties[0]/det_J,//for mesh motion penalize small elements
92  nu=materialProperties[1];
93 
94  const double shear = E/(1.0+nu);
95  const double bulk = shear*(nu/(1.0-2.0*nu));
96 
97  for (int i=0;i<nSymTen;i++)
98  for (int j=0;j<nSymTen;j++)
99  dstress[i*nSymTen+j] = 0.0;
100  stress[sXX] = shear*strain[sXX] + bulk*strainTrace;
101  dstress[sXX*nSymTen+sXX] = bulk + shear;
102  dstress[sXX*nSymTen+sYY] = bulk;
103  /* dstress[sXX*nSymTen+sZZ] = bulk; */
104 
105  stress[sYY] = shear*strain[sYY] + bulk*strainTrace;
106  dstress[sYY*nSymTen+sXX] = bulk;
107  dstress[sYY*nSymTen+sYY] = bulk + shear;
108  /* dstress[sYY*nSymTen+sZZ] = bulk; */
109 
110  /* stress[sZZ] = shear*strain[sZZ] + bulk*strainTrace; */
111  /* dstress[sZZ*nSymTen+sXX] = bulk; */
112  /* dstress[sZZ*nSymTen+sYY] = bulk; */
113  /* dstress[sZZ*nSymTen+sZZ] = bulk + shear; */
114 
115  /* stress[sYZ] = shear*0.5*strain[sYZ];//the 1/2 comes from the Voigt notation */
116  /* dstress[sYZ*nSymTen+sYZ] = shear*0.5; */
117 
118  /* stress[sXZ] = shear*0.5*strain[sXZ];//the 1/2 comes from the Voigt notation */
119  /* dstress[sXZ*nSymTen+sXZ] = shear*0.5; */
120 
121  stress[sXY] = shear*0.5*strain[sXY];//the 1/2 comes from the Voigt notation
122  dstress[sXY*nSymTen+sXY] = shear*0.5;
123  }
124 
125  inline void exteriorNumericalStressFlux(const int& isDOFBoundary_u,
126  const int& isDOFBoundary_v,
127  /* const int& isDOFBoundary_w, */
128  const int& isStressFluxBoundary_u,
129  const int& isStressFluxBoundary_v,
130  /* const int& isStressFluxBoundary_w, */
131  const double& penalty,
132  const double& u,
133  const double& v,
134  /* const double& w, */
135  const double& bc_u,
136  const double& bc_v,
137  /* const double& bc_w, */
138  const double& bc_stressFlux_u,
139  const double& bc_stressFlux_v,
140  /* const double& bc_stressFlux_w, */
141  const double* stress,
142  const double* normal,
143  double& stressFlux_u,
144  double& stressFlux_v/* , */
145  /* double& stressFlux_w */)
146  {
147  if (isDOFBoundary_u == 1)
148  {
149  stressFlux_u = -(stress[sXX]*normal[X] + stress[sXY]*normal[Y] /* + stress[sXZ]*normal[Z] */ - penalty*(u - bc_u));
150  }
151  else if(isStressFluxBoundary_u == 1)
152  {
153  stressFlux_u = bc_stressFlux_u;
154  }
155  else
156  {
157  stressFlux_u = 0.0;
158  }
159 
160  if (isDOFBoundary_v == 1)
161  {
162  stressFlux_v = -(stress[sYX]*normal[X] + stress[sYY]*normal[Y] /* + stress[sYZ]*normal[Z] */ - penalty*(v - bc_v));
163  }
164  else if(isStressFluxBoundary_v == 1)
165  {
166  stressFlux_v = bc_stressFlux_v;
167  }
168  else
169  {
170  stressFlux_v = 0.0;
171  }
172 
173  /* if (isDOFBoundary_w == 1) */
174  /* { */
175  /* stressFlux_w = -(stress[sZX]*normal[X] + stress[sZY]*normal[Y] + stress[sZZ]*normal[Z] - penalty*(w - bc_w)); */
176  /* } */
177  /* else if(isStressFluxBoundary_w == 1) */
178  /* { */
179  /* stressFlux_w = bc_stressFlux_w; */
180  /* } */
181  /* else */
182  /* { */
183  /* stressFlux_w = 0.0; */
184  /* } */
185  }
186 
187  inline void exteriorNumericalStressFluxJacobian(const int& isDOFBoundary_u,
188  const int& isDOFBoundary_v,
189  /* const int& isDOFBoundary_w, */
190  const double* normal,
191  const double* dstress,
192  const double& penalty,
193  const double& disp_trial,
194  const double* disp_grad_trial,
195  double& dstressFlux_u_u,
196  double& dstressFlux_u_v,
197  /* double& dstressFlux_u_w, */
198  double& dstressFlux_v_u,
199  double& dstressFlux_v_v/* , */
200  /* double& dstressFlux_v_w, */
201  /* double& dstressFlux_w_u, */
202  /* double& dstressFlux_w_v, */
203  /* double& dstressFlux_w_w */)
204  {
205  //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
206  if (isDOFBoundary_u == 1)
207  {
208  dstressFlux_u_u = -(
209  (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]+
210  (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]+
211  /* (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] */
212  -
213  penalty*disp_trial);
214  dstressFlux_u_v = -(
215  (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]+
216  (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]/* + */
217  /* (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] */);
218  /* dstressFlux_u_w = -( */
219  /* (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]+ */
220  /* (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]+ */
221  /* (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]); */
222  }
223  else
224  {
225  dstressFlux_u_u = 0.0;
226  dstressFlux_u_v = 0.0;
227  /* dstressFlux_u_w = 0.0; */
228  }
229 
230  if (isDOFBoundary_v == 1)
231  {
232  dstressFlux_v_u = -(
233  (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]+
234  (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]/* + */
235  /* (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] */);
236  dstressFlux_v_v = -(
237  (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]+
238  (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]/* + */
239  /* (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] */
240  -
241  penalty*disp_trial);
242  /* dstressFlux_v_w = -( */
243  /* (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]+ */
244  /* (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]+ */
245  /* (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]); */
246  }
247  else
248  {
249  dstressFlux_v_u = 0.0;
250  dstressFlux_v_v = 0.0;
251  /* dstressFlux_v_w = 0.0; */
252  }
253 
254  /* if (isDOFBoundary_w == 1) */
255  /* { */
256  /* dstressFlux_w_u = -( */
257  /* (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]+ */
258  /* (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]+ */
259  /* (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]); */
260  /* dstressFlux_w_v = -( */
261  /* (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]+ */
262  /* (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]+ */
263  /* (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]); */
264  /* dstressFlux_w_w = -( */
265  /* (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]+ */
266  /* (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]+ */
267  /* (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] */
268  /* - */
269  /* penalty*disp_trial); */
270  /* } */
271  /* else */
272  /* { */
273  /* dstressFlux_w_u = 0.0; */
274  /* dstressFlux_w_v = 0.0; */
275  /* dstressFlux_w_w = 0.0; */
276  /* } */
277  }
278 
279 
280  virtual void calculateResidual(arguments_dict& args)
281  {
282  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
283  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
284  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
285  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
286  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
287  xt::pyarray<double>& disp_trial_ref = args.array<double>("disp_trial_ref");
288  xt::pyarray<double>& disp_grad_trial_ref = args.array<double>("disp_grad_trial_ref");
289  xt::pyarray<double>& disp_test_ref = args.array<double>("disp_test_ref");
290  xt::pyarray<double>& disp_grad_test_ref = args.array<double>("disp_grad_test_ref");
291  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
292  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
293  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
294  xt::pyarray<double>& disp_trial_trace_ref = args.array<double>("disp_trial_trace_ref");
295  xt::pyarray<double>& disp_grad_trial_trace_ref = args.array<double>("disp_grad_trial_trace_ref");
296  xt::pyarray<double>& disp_test_trace_ref = args.array<double>("disp_test_trace_ref");
297  xt::pyarray<double>& disp_grad_test_trace_ref = args.array<double>("disp_grad_test_trace_ref");
298  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
299  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
300  int nElements_global = args.scalar<int>("nElements_global");
301  xt::pyarray<int>& materialTypes = args.array<int>("materialTypes");
302  int nMaterialProperties = args.scalar<int>("nMaterialProperties");
303  xt::pyarray<double>& materialProperties = args.array<double>("materialProperties");
304  xt::pyarray<int>& disp_l2g = args.array<int>("disp_l2g");
305  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
306  xt::pyarray<double>& v_dof = args.array<double>("v_dof");
307  xt::pyarray<double>& w_dof = args.array<double>("w_dof");
308  xt::pyarray<double>& bodyForce = args.array<double>("bodyForce");
309  int offset_u = args.scalar<int>("offset_u");
310  int offset_v = args.scalar<int>("offset_v");
311  int offset_w = args.scalar<int>("offset_w");
312  int stride_u = args.scalar<int>("stride_u");
313  int stride_v = args.scalar<int>("stride_v");
314  int stride_w = args.scalar<int>("stride_w");
315  xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
316  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
317  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
318  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
319  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
320  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
321  xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
322  xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
323  xt::pyarray<int>& isStressFluxBoundary_u = args.array<int>("isStressFluxBoundary_u");
324  xt::pyarray<int>& isStressFluxBoundary_v = args.array<int>("isStressFluxBoundary_v");
325  xt::pyarray<int>& isStressFluxBoundary_w = args.array<int>("isStressFluxBoundary_w");
326  xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
327  xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
328  xt::pyarray<double>& ebqe_bc_w_ext = args.array<double>("ebqe_bc_w_ext");
329  xt::pyarray<double>& ebqe_bc_stressFlux_u_ext = args.array<double>("ebqe_bc_stressFlux_u_ext");
330  xt::pyarray<double>& ebqe_bc_stressFlux_v_ext = args.array<double>("ebqe_bc_stressFlux_v_ext");
331  xt::pyarray<double>& ebqe_bc_stressFlux_w_ext = args.array<double>("ebqe_bc_stressFlux_w_ext");
332  //
333  //loop over elements to compute volume integrals and load them into element and global residual
334  //
335  for(int eN=0;eN<nElements_global;eN++)
336  {
337  //declare local storage for element residual and initialize
338  register double
339  elementResidual_u[nDOF_test_element],
340  elementResidual_v[nDOF_test_element]/* , */
341  /* elementResidual_w[nDOF_test_element] */;
342  for (int i=0;i<nDOF_test_element;i++)
343  {
344  elementResidual_u[i]=0.0;
345  elementResidual_v[i]=0.0;
346  /* elementResidual_w[i]=0.0; */
347  }//i
348  //
349  //loop over quadrature points and compute integrands
350  //
351  for(int k=0;k<nQuadraturePoints_element;k++)
352  {
353  //compute indices and declare local storage
354  register int //eN_k = eN*nQuadraturePoints_element+k,
355  //eN_k_nSpace=eN_k*nSpace,
356  eN_nDOF_trial_element = eN*nDOF_trial_element;
357  register double u=0.0,v=0.0,/* w=0.0, */
358  D[nSpace*nSpace],
359  *grad_u(&D[0]),
360  *grad_v(&D[nSpace]),
361  /* *grad_w(&D[2*nSpace]), */
362  jac[nSpace*nSpace],
363  jacDet,
364  jacInv[nSpace*nSpace],
365  disp_grad_trial[nDOF_trial_element*nSpace],
366  disp_test_dV[nDOF_trial_element],
367  disp_grad_test_dV[nDOF_test_element*nSpace],
368  dV,x,y,z,
369  G[nSpace*nSpace],G_dd_G,tr_G,
370  strain[ck.nSymTen],stress[ck.nSymTen],dstress[ck.nSymTen*ck.nSymTen];
371  //get jacobian, etc for mapping reference element
372  ck.calculateMapping_element(eN,
373  k,
374  mesh_dof.data(),
375  mesh_l2g.data(),
376  mesh_trial_ref.data(),
377  mesh_grad_trial_ref.data(),
378  jac,
379  jacDet,
380  jacInv,
381  x,y,z);
382  //get the physical integration weight
383  dV = fabs(jacDet)*dV_ref.data()[k];
384  ck.calculateG(jacInv,G,G_dd_G,tr_G);
385  //get the trial function gradients
386  ck.gradTrialFromRef(&disp_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,disp_grad_trial);
387  //get the solution
388  ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],u);
389  ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],v);
390  /* ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],w); */
391  //get the solution gradients
392  ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_u);
393  ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_v);
394  /* ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_w); */
395  //precalculate test function products with integration weights
396  for (int j=0;j<nDOF_trial_element;j++)
397  {
398  disp_test_dV[j] = disp_test_ref.data()[k*nDOF_trial_element+j]*dV;
399  for (int I=0;I<nSpace;I++)
400  {
401  disp_grad_test_dV[j*nSpace+I] = disp_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin
402  }
403  }
404  //save displacement at quadrature points for other models to use
405  //q_displacement[eN_k_nSpace+0]=u;
406  //q_displacement[eN_k_nSpace+1]=v;
407  //q_displacement[eN_k_nSpace+2]=w;
408 
409  calculateStrain(D,strain);
410  evaluateCoefficients(fabs(jacDet),
411  &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
412  strain,
413  stress,
414  dstress);
415  //
416  //update element residual
417  //
418  for(int i=0;i<nDOF_test_element;i++)
419  {
420  register int i_nSpace=i*nSpace;
421 
422  elementResidual_u[i] += ck.Stress_u_weak(stress,&disp_grad_test_dV[i_nSpace]);
423  elementResidual_v[i] += ck.Stress_v_weak(stress,&disp_grad_test_dV[i_nSpace]);
424  /* elementResidual_w[i] += ck.Stress_w_weak(stress,&disp_grad_test_dV[i_nSpace]); */
425  }//i
426  }
427  //
428  //load element into global residual and save element residual
429  //
430  for(int i=0;i<nDOF_test_element;i++)
431  {
432  register int eN_i=eN*nDOF_test_element+i;
433 
434  globalResidual.data()[offset_u+stride_u*disp_l2g.data()[eN_i]] += elementResidual_u[i];
435  globalResidual.data()[offset_v+stride_v*disp_l2g.data()[eN_i]] += elementResidual_v[i];
436  /* globalResidual.data()[offset_w+stride_w*disp_l2g.data()[eN_i]] += elementResidual_w[i]; */
437  }//i
438  }//elements
439  //
440  //loop over exterior element boundaries to calculate surface integrals and load into element and global residuals
441  //
442  //ebNE is the Exterior element boundary INdex
443  //ebN is the element boundary INdex
444  //eN is the element index
445  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
446  {
447  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
448  eN = elementBoundaryElementsArray.data()[ebN*2+0],
449  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0],
450  eN_nDOF_trial_element = eN*nDOF_trial_element;
451  register double elementResidual_u[nDOF_test_element],
452  elementResidual_v[nDOF_test_element]/* , */
453  /* elementResidual_w[nDOF_test_element] */;
454  for (int i=0;i<nDOF_test_element;i++)
455  {
456  elementResidual_u[i]=0.0;
457  elementResidual_v[i]=0.0;
458  /* elementResidual_w[i]=0.0; */
459  }
460  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
461  {
462  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
463  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
464  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
465  register double u_ext=0.0,
466  v_ext=0.0,
467  /* w_ext=0.0, */
468  D[nSpace*nSpace],
469  *grad_u_ext(&D[0]),
470  *grad_v_ext(&D[nSpace]),
471  /* *grad_w_ext(&D[2*nSpace]), */
472  bc_u_ext=0.0,
473  bc_v_ext=0.0,
474  bc_w_ext=0.0,
475  jac_ext[nSpace*nSpace],
476  jacDet_ext,
477  jacInv_ext[nSpace*nSpace],
478  boundaryJac[nSpace*(nSpace-1)],
479  metricTensor[(nSpace-1)*(nSpace-1)],
480  metricTensorDetSqrt,
481  dS,disp_test_dS[nDOF_test_element],
482  disp_grad_trial_trace[nDOF_trial_element*nSpace],
483  normal[2],x_ext,y_ext,z_ext,
484  G[nSpace*nSpace],G_dd_G,tr_G,h_penalty,
485  strain[ck.nSymTen],stress[ck.nSymTen],dstress[ck.nSymTen*ck.nSymTen],
486  stressFlux_u,stressFlux_v,stressFlux_w;
487  //compute information about mapping from reference element to physical element
488  ck.calculateMapping_elementBoundary(eN,
489  ebN_local,
490  kb,
491  ebN_local_kb,
492  mesh_dof.data(),
493  mesh_l2g.data(),
494  mesh_trial_trace_ref.data(),
495  mesh_grad_trial_trace_ref.data(),
496  boundaryJac_ref.data(),
497  jac_ext,
498  jacDet_ext,
499  jacInv_ext,
500  boundaryJac,
501  metricTensor,
502  metricTensorDetSqrt,
503  normal_ref.data(),
504  normal,
505  x_ext,y_ext,z_ext);
506  dS = metricTensorDetSqrt*dS_ref.data()[kb];
507  //get the metric tensor
508  //cek todo use symmetry
509  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
510  //compute shape and solution information
511  //shape
512  ck.gradTrialFromRef(&disp_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,disp_grad_trial_trace);
513  //solution and gradients
514  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);
515  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);
516  /* 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); */
517  ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_u_ext);
518  ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_v_ext);
519  /* ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_w_ext); */
520  //precalculate test function products with integration weights
521  for (int j=0;j<nDOF_trial_element;j++)
522  {
523  disp_test_dS[j] = disp_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
524  }
525  //
526  //load the boundary values
527  //
528  bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
529  bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
530  /* bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext; */
531  //
532  //calculate the pde coefficients using the solution and the boundary values for the solution
533  //
534  calculateStrain(D,strain);
535  evaluateCoefficients(fabs(jacDet_ext),
536  &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
537  strain,
538  stress,
539  dstress);
540  //
541  //calculate the numerical fluxes
542  //
543  //cek debug
544  //ebqe_penalty_ext[ebNE_kb] = 10.0;
545  //
546  ck.calculateGScale(G,normal,h_penalty);
547  //cek hack
548  h_penalty = 100.0/h_penalty;
549  exteriorNumericalStressFlux(isDOFBoundary_u.data()[ebNE_kb],
550  isDOFBoundary_v.data()[ebNE_kb],
551  /* isDOFBoundary_w.data()[ebNE_kb], */
552  isStressFluxBoundary_u.data()[ebNE_kb],
553  isStressFluxBoundary_v.data()[ebNE_kb],
554  /* isStressFluxBoundary_w.data()[ebNE_kb], */
555  h_penalty,
556  u_ext,
557  v_ext,
558  /* w_ext, */
559  bc_u_ext,
560  bc_v_ext,
561  /* bc_w_ext, */
562  ebqe_bc_stressFlux_u_ext.data()[ebNE_kb],
563  ebqe_bc_stressFlux_v_ext.data()[ebNE_kb],
564  /* ebqe_bc_stressFlux_w_ext.data()[ebNE_kb], */
565  stress,
566  normal,
567  stressFlux_u,
568  stressFlux_v/* , */
569  /* stressFlux_w */);
570  //
571  //update residuals
572  //
573  for (int i=0;i<nDOF_test_element;i++)
574  {
575  elementResidual_u[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_u,disp_test_dS[i]);
576  elementResidual_v[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_v,disp_test_dS[i]);
577  /* elementResidual_w[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_w,disp_test_dS[i]); */
578  }//i
579  }//kb
580  //
581  //update the element and global residual storage
582  //
583  for (int i=0;i<nDOF_test_element;i++)
584  {
585  int eN_i = eN*nDOF_test_element+i;
586 
587  globalResidual.data()[offset_u+stride_u*disp_l2g.data()[eN_i]]+=elementResidual_u[i];
588  globalResidual.data()[offset_v+stride_v*disp_l2g.data()[eN_i]]+=elementResidual_v[i];
589  /* globalResidual.data()[offset_w+stride_w*disp_l2g.data()[eN_i]]+=elementResidual_w[i]; */
590  }//i
591  }//ebNE
592  }
593 
594  virtual void calculateJacobian(arguments_dict& args)
595  {
596  xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
597  xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
598  xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
599  xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
600  xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
601  xt::pyarray<double>& disp_trial_ref = args.array<double>("disp_trial_ref");
602  xt::pyarray<double>& disp_grad_trial_ref = args.array<double>("disp_grad_trial_ref");
603  xt::pyarray<double>& disp_test_ref = args.array<double>("disp_test_ref");
604  xt::pyarray<double>& disp_grad_test_ref = args.array<double>("disp_grad_test_ref");
605  xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
606  xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
607  xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
608  xt::pyarray<double>& disp_trial_trace_ref = args.array<double>("disp_trial_trace_ref");
609  xt::pyarray<double>& disp_grad_trial_trace_ref = args.array<double>("disp_grad_trial_trace_ref");
610  xt::pyarray<double>& disp_test_trace_ref = args.array<double>("disp_test_trace_ref");
611  xt::pyarray<double>& disp_grad_test_trace_ref = args.array<double>("disp_grad_test_trace_ref");
612  xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
613  xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
614  int nElements_global = args.scalar<int>("nElements_global");
615  xt::pyarray<int>& materialTypes = args.array<int>("materialTypes");
616  int nMaterialProperties = args.scalar<int>("nMaterialProperties");
617  xt::pyarray<double>& materialProperties = args.array<double>("materialProperties");
618  xt::pyarray<int>& disp_l2g = args.array<int>("disp_l2g");
619  xt::pyarray<double>& u_dof = args.array<double>("u_dof");
620  xt::pyarray<double>& v_dof = args.array<double>("v_dof");
621  xt::pyarray<double>& w_dof = args.array<double>("w_dof");
622  xt::pyarray<double>& bodyForce = args.array<double>("bodyForce");
623  xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
624  xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
625  xt::pyarray<int>& csrRowIndeces_u_v = args.array<int>("csrRowIndeces_u_v");
626  xt::pyarray<int>& csrColumnOffsets_u_v = args.array<int>("csrColumnOffsets_u_v");
627  xt::pyarray<int>& csrRowIndeces_u_w = args.array<int>("csrRowIndeces_u_w");
628  xt::pyarray<int>& csrColumnOffsets_u_w = args.array<int>("csrColumnOffsets_u_w");
629  xt::pyarray<int>& csrRowIndeces_v_u = args.array<int>("csrRowIndeces_v_u");
630  xt::pyarray<int>& csrColumnOffsets_v_u = args.array<int>("csrColumnOffsets_v_u");
631  xt::pyarray<int>& csrRowIndeces_v_v = args.array<int>("csrRowIndeces_v_v");
632  xt::pyarray<int>& csrColumnOffsets_v_v = args.array<int>("csrColumnOffsets_v_v");
633  xt::pyarray<int>& csrRowIndeces_v_w = args.array<int>("csrRowIndeces_v_w");
634  xt::pyarray<int>& csrColumnOffsets_v_w = args.array<int>("csrColumnOffsets_v_w");
635  xt::pyarray<int>& csrRowIndeces_w_u = args.array<int>("csrRowIndeces_w_u");
636  xt::pyarray<int>& csrColumnOffsets_w_u = args.array<int>("csrColumnOffsets_w_u");
637  xt::pyarray<int>& csrRowIndeces_w_v = args.array<int>("csrRowIndeces_w_v");
638  xt::pyarray<int>& csrColumnOffsets_w_v = args.array<int>("csrColumnOffsets_w_v");
639  xt::pyarray<int>& csrRowIndeces_w_w = args.array<int>("csrRowIndeces_w_w");
640  xt::pyarray<int>& csrColumnOffsets_w_w = args.array<int>("csrColumnOffsets_w_w");
641  xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
642  int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
643  xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
644  xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
645  xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
646  xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
647  xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
648  xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
649  xt::pyarray<int>& isStressFluxBoundary_u = args.array<int>("isStressFluxBoundary_u");
650  xt::pyarray<int>& isStressFluxBoundary_v = args.array<int>("isStressFluxBoundary_v");
651  xt::pyarray<int>& isStressFluxBoundary_w = args.array<int>("isStressFluxBoundary_w");
652  xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
653  xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.array<int>("csrColumnOffsets_eb_u_v");
654  xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.array<int>("csrColumnOffsets_eb_u_w");
655  xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.array<int>("csrColumnOffsets_eb_v_u");
656  xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.array<int>("csrColumnOffsets_eb_v_v");
657  xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.array<int>("csrColumnOffsets_eb_v_w");
658  xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.array<int>("csrColumnOffsets_eb_w_u");
659  xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.array<int>("csrColumnOffsets_eb_w_v");
660  xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.array<int>("csrColumnOffsets_eb_w_w");
661  CompKernel<nSpace,nDOF_mesh_trial_element,nDOF_trial_element,nDOF_test_element> ck;
662  const int nSymTen(ck.nSymTen);
663  //
664  //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
665  //
666  for(int eN=0;eN<nElements_global;eN++)
667  {
668  register double
669  elementJacobian_u_u[nDOF_test_element][nDOF_trial_element],
670  elementJacobian_u_v[nDOF_test_element][nDOF_trial_element],
671  /* elementJacobian_u_w[nDOF_test_element][nDOF_trial_element], */
672  elementJacobian_v_u[nDOF_test_element][nDOF_trial_element],
673  elementJacobian_v_v[nDOF_test_element][nDOF_trial_element]/* , */
674  /* elementJacobian_v_w[nDOF_test_element][nDOF_trial_element], */
675  /* elementJacobian_w_u[nDOF_test_element][nDOF_trial_element], */
676  /* elementJacobian_w_v[nDOF_test_element][nDOF_trial_element], */
677  /* elementJacobian_w_w[nDOF_test_element][nDOF_trial_element] */;
678  for (int i=0;i<nDOF_test_element;i++)
679  for (int j=0;j<nDOF_trial_element;j++)
680  {
681  elementJacobian_u_u[i][j]=0.0;
682  elementJacobian_u_v[i][j]=0.0;
683  /* elementJacobian_u_w[i][j]=0.0; */
684  elementJacobian_v_u[i][j]=0.0;
685  elementJacobian_v_v[i][j]=0.0;
686  /* elementJacobian_v_w[i][j]=0.0; */
687  /* elementJacobian_w_u[i][j]=0.0; */
688  /* elementJacobian_w_v[i][j]=0.0; */
689  /* elementJacobian_w_w[i][j]=0.0; */
690  }
691  for (int k=0;k<nQuadraturePoints_element;k++)
692  {
693  const int
694  eN_nDOF_trial_element = eN*nDOF_trial_element; //index to a vector at a quadrature point
695 
696  //declare local storage
697  register double u=0.0,v=0.0,/* w=0.0, */
698  D[nSpace*nSpace],
699  *grad_u(&D[0]),
700  *grad_v(&D[nSpace]),
701  /* *grad_w(&D[2*nSpace]), */
702  jac[nSpace*nSpace],
703  jacDet,
704  jacInv[nSpace*nSpace],
705  disp_grad_trial[nDOF_trial_element*nSpace],
706  dV,
707  disp_test_dV[nDOF_test_element],
708  disp_grad_test_dV[nDOF_test_element*nSpace],
709  x,y,z,
710  G[nSpace*nSpace],G_dd_G,tr_G,
711  strain[nSymTen],stress[nSymTen],dstress[nSpace*nSpace*nSpace*nSpace];
712  //get jacobian, etc for mapping reference element
713  ck.calculateMapping_element(eN,
714  k,
715  mesh_dof.data(),
716  mesh_l2g.data(),
717  mesh_trial_ref.data(),
718  mesh_grad_trial_ref.data(),
719  jac,
720  jacDet,
721  jacInv,
722  x,y,z);
723  //get the physical integration weight
724  dV = fabs(jacDet)*dV_ref.data()[k];
725  ck.calculateG(jacInv,G,G_dd_G,tr_G);
726  //get the trial function gradients
727  ck.gradTrialFromRef(&disp_grad_trial_ref.data()[k*nDOF_trial_element*nSpace],jacInv,disp_grad_trial);
728  //get the solution
729  ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],u);
730  ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],v);
731  /* ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_ref.data()[k*nDOF_trial_element],w); */
732  //get the solution gradients
733  ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_u);
734  ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_v);
735  /* ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial,grad_w); */
736  //precalculate test function products with integration weights
737  for (int j=0;j<nDOF_trial_element;j++)
738  {
739  disp_test_dV[j] = disp_test_ref.data()[k*nDOF_trial_element+j]*dV;
740  for (int I=0;I<nSpace;I++)
741  {
742  disp_grad_test_dV[j*nSpace+I] = disp_grad_trial[j*nSpace+I]*dV;//cek warning won't work for Petrov-Galerkin}
743  }
744  }
745  calculateStrain(D,strain);
746  evaluateCoefficients(fabs(jacDet),
747  &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
748  strain,
749  stress,
750  dstress);
751  //
752  //omit for now
753  //
754  for(int i=0;i<nDOF_test_element;i++)
755  {
756  register int i_nSpace = i*nSpace;
757  for(int j=0;j<nDOF_trial_element;j++)
758  {
759  register int j_nSpace = j*nSpace;
760 
761  elementJacobian_u_u[i][j] += ck.StressJacobian_u_u_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
762  elementJacobian_u_v[i][j] += ck.StressJacobian_u_v_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
763  /* elementJacobian_u_w[i][j] += ck.StressJacobian_u_w_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]); */
764 
765  elementJacobian_v_u[i][j] += ck.StressJacobian_v_u_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
766  elementJacobian_v_v[i][j] += ck.StressJacobian_v_v_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]);
767  /* elementJacobian_v_w[i][j] += ck.StressJacobian_v_w_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]); */
768 
769  /* elementJacobian_w_u[i][j] += ck.StressJacobian_w_u_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]); */
770  /* elementJacobian_w_v[i][j] += ck.StressJacobian_w_v_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]); */
771  /* elementJacobian_w_w[i][j] += ck.StressJacobian_w_w_weak(dstress,&disp_grad_trial[j_nSpace],&disp_grad_test_dV[i_nSpace]); */
772  }//j
773  }//i
774  }//k
775  //
776  //load into element Jacobian into global Jacobian
777  //
778  for (int i=0;i<nDOF_test_element;i++)
779  {
780  register int eN_i = eN*nDOF_test_element+i;
781  for (int j=0;j<nDOF_trial_element;j++)
782  {
783  register int eN_i_j = eN_i*nDOF_trial_element+j;
784 
785  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_u_u.data()[eN_i_j]] += elementJacobian_u_u[i][j];
786  globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_u_v.data()[eN_i_j]] += elementJacobian_u_v[i][j];
787  /* globalJacobian.data()[csrRowIndeces_u_w.data()[eN_i] + csrColumnOffsets_u_w.data()[eN_i_j]] += elementJacobian_u_w[i][j]; */
788 
789  globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_v_u.data()[eN_i_j]] += elementJacobian_v_u[i][j];
790  globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_v_v.data()[eN_i_j]] += elementJacobian_v_v[i][j];
791  /* globalJacobian.data()[csrRowIndeces_v_w.data()[eN_i] + csrColumnOffsets_v_w.data()[eN_i_j]] += elementJacobian_v_w[i][j]; */
792 
793  /* globalJacobian.data()[csrRowIndeces_w_u.data()[eN_i] + csrColumnOffsets_w_u.data()[eN_i_j]] += elementJacobian_w_u[i][j]; */
794  /* globalJacobian.data()[csrRowIndeces_w_v.data()[eN_i] + csrColumnOffsets_w_v.data()[eN_i_j]] += elementJacobian_w_v[i][j]; */
795  /* globalJacobian.data()[csrRowIndeces_w_w.data()[eN_i] + csrColumnOffsets_w_w.data()[eN_i_j]] += elementJacobian_w_w[i][j]; */
796  }//j
797  }//i
798  }//elements
799  //
800  //loop over exterior element boundaries to compute the surface integrals and load them into the global Jacobian
801  //
802  for (int ebNE = 0; ebNE < nExteriorElementBoundaries_global; ebNE++)
803  {
804  register int ebN = exteriorElementBoundariesArray.data()[ebNE],
805  eN = elementBoundaryElementsArray.data()[ebN*2+0],
806  eN_nDOF_trial_element = eN*nDOF_trial_element,
807  ebN_local = elementBoundaryLocalElementBoundariesArray.data()[ebN*2+0];
808  for (int kb=0;kb<nQuadraturePoints_elementBoundary;kb++)
809  {
810  register int ebNE_kb = ebNE*nQuadraturePoints_elementBoundary+kb,
811  ebN_local_kb = ebN_local*nQuadraturePoints_elementBoundary+kb,
812  ebN_local_kb_nSpace = ebN_local_kb*nSpace;
813 
814  register double
815  u_ext=0.0,
816  v_ext=0.0,
817  /* w_ext=0.0, */
818  D_ext[nSpace*nSpace],
819  *grad_u_ext(&D_ext[0]),
820  *grad_v_ext(&D_ext[nSpace]),
821  /* *grad_w_ext(&D_ext[2*nSpace]), */
822  fluxJacobian_u_u[nDOF_trial_element],
823  fluxJacobian_u_v[nDOF_trial_element],
824  /* fluxJacobian_u_w[nDOF_trial_element], */
825  fluxJacobian_v_u[nDOF_trial_element],
826  fluxJacobian_v_v[nDOF_trial_element],
827  /* fluxJacobian_v_w[nDOF_trial_element], */
828  /* fluxJacobian_w_u[nDOF_trial_element], */
829  /* fluxJacobian_w_v[nDOF_trial_element], */
830  /* fluxJacobian_w_w[nDOF_trial_element], */
831  jac_ext[nSpace*nSpace],
832  jacDet_ext,
833  jacInv_ext[nSpace*nSpace],
834  boundaryJac[nSpace*(nSpace-1)],
835  metricTensor[(nSpace-1)*(nSpace-1)],
836  metricTensorDetSqrt,
837  disp_grad_trial_trace[nDOF_trial_element*nSpace],
838  dS,
839  disp_test_dS[nDOF_test_element],
840  normal[2],
841  x_ext,y_ext,z_ext,
842  G[nSpace*nSpace],G_dd_G,tr_G,h_penalty,
843  strain[nSymTen],
844  stress[nSymTen],
845  dstress[nSymTen*nSymTen];
846  ck.calculateMapping_elementBoundary(eN,
847  ebN_local,
848  kb,
849  ebN_local_kb,
850  mesh_dof.data(),
851  mesh_l2g.data(),
852  mesh_trial_trace_ref.data(),
853  mesh_grad_trial_trace_ref.data(),
854  boundaryJac_ref.data(),
855  jac_ext,
856  jacDet_ext,
857  jacInv_ext,
858  boundaryJac,
859  metricTensor,
860  metricTensorDetSqrt,
861  normal_ref.data(),
862  normal,
863  x_ext,y_ext,z_ext);
864  dS = metricTensorDetSqrt*dS_ref.data()[kb];
865  ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
866  //compute shape and solution information
867  //shape
868  ck.gradTrialFromRef(&disp_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,disp_grad_trial_trace);
869  //solution and gradients
870  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);
871  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);
872  /* 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); */
873  ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_u_ext);
874  ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_v_ext);
875  /* ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_w_ext); */
876  //precalculate test function products with integration weights
877  for (int j=0;j<nDOF_trial_element;j++)
878  {
879  disp_test_dS[j] = disp_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
880  }
881  //
882  //calculate the internal and external trace of the pde coefficients
883  //
884  calculateStrain(D_ext,strain);
885  evaluateCoefficients(fabs(jacDet_ext),
886  &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
887  strain,
888  stress,
889  dstress);
890  //
891  //calculate the flux jacobian
892  //
893  ck.calculateGScale(G,normal,h_penalty);
894  //cek hack
895  h_penalty = 100.0/h_penalty;
896  // for (int II=0;II<nSymTen;II++)
897  // for (int JJ=0;<nSymTen;JJ++)
898  // dstress[II*nSymTen+JJ] = 0.0;
899  for (int j=0;j<nDOF_trial_element;j++)
900  {
901  register int j_nSpace = j*nSpace;
902 
903  exteriorNumericalStressFluxJacobian(isDOFBoundary_u.data()[ebNE_kb],
904  isDOFBoundary_v.data()[ebNE_kb],
905  /* isDOFBoundary_w.data()[ebNE_kb], */
906  normal,
907  dstress,
908  h_penalty,
909  disp_trial_trace_ref.data()[ebN_local_kb*nDOF_trial_element+j],
910  &disp_grad_trial_trace[j_nSpace],
911  fluxJacobian_u_u[j],
912  fluxJacobian_u_v[j],
913  /* fluxJacobian_u_w[j], */
914  fluxJacobian_v_u[j],
915  fluxJacobian_v_v[j]/* , */
916  /* fluxJacobian_v_w[j], */
917  /* fluxJacobian_w_u[j], */
918  /* fluxJacobian_w_v[j], */
919  /* fluxJacobian_w_w[j] */);
920  }//j
921  //
922  //update the global Jacobian from the flux Jacobian
923  //
924  for (int i=0;i<nDOF_test_element;i++)
925  {
926  register int eN_i = eN*nDOF_test_element+i;
927  for (int j=0;j<nDOF_trial_element;j++)
928  {
929  register int ebN_i_j = ebN*4*nDOF_test_X_trial_element + i*nDOF_trial_element + j;
930 
931  globalJacobian.data()[csrRowIndeces_u_u.data()[eN_i] + csrColumnOffsets_eb_u_u.data()[ebN_i_j]] += fluxJacobian_u_u[j]*disp_test_dS[i];
932  globalJacobian.data()[csrRowIndeces_u_v.data()[eN_i] + csrColumnOffsets_eb_u_v.data()[ebN_i_j]] += fluxJacobian_u_v[j]*disp_test_dS[i];
933  /* globalJacobian.data()[csrRowIndeces_u_w.data()[eN_i] + csrColumnOffsets_eb_u_w.data()[ebN_i_j]] += fluxJacobian_u_w[j]*disp_test_dS[i]; */
934 
935  globalJacobian.data()[csrRowIndeces_v_u.data()[eN_i] + csrColumnOffsets_eb_v_u.data()[ebN_i_j]] += fluxJacobian_v_u[j]*disp_test_dS[i];
936  globalJacobian.data()[csrRowIndeces_v_v.data()[eN_i] + csrColumnOffsets_eb_v_v.data()[ebN_i_j]] += fluxJacobian_v_v[j]*disp_test_dS[i];
937  /* globalJacobian.data()[csrRowIndeces_v_w.data()[eN_i] + csrColumnOffsets_eb_v_w.data()[ebN_i_j]] += fluxJacobian_v_w[j]*disp_test_dS[i]; */
938 
939  /* globalJacobian.data()[csrRowIndeces_w_u.data()[eN_i] + csrColumnOffsets_eb_w_u.data()[ebN_i_j]] += fluxJacobian_w_u[j]*disp_test_dS[i]; */
940  /* globalJacobian.data()[csrRowIndeces_w_v.data()[eN_i] + csrColumnOffsets_eb_w_v.data()[ebN_i_j]] += fluxJacobian_w_v[j]*disp_test_dS[i]; */
941  /* globalJacobian.data()[csrRowIndeces_w_w.data()[eN_i] + csrColumnOffsets_eb_w_w.data()[ebN_i_j]] += fluxJacobian_w_w[j]*disp_test_dS[i]; */
942  }//j
943  }//i
944  }//kb
945  }//ebNE
946  }//computeJacobian
947  };//MoveMesh2D
948 
949  inline MoveMesh2D_base* newMoveMesh2D(int nSpaceIn,
950  int nQuadraturePoints_elementIn,
951  int nDOF_mesh_trial_elementIn,
952  int nDOF_trial_elementIn,
953  int nDOF_test_elementIn,
954  int nQuadraturePoints_elementBoundaryIn,
955  int CompKernelFlag)
956  {
957  return proteus::chooseAndAllocateDiscretization2D<MoveMesh2D_base,MoveMesh2D,CompKernel>(nSpaceIn,
958  nQuadraturePoints_elementIn,
959  nDOF_mesh_trial_elementIn,
960  nDOF_trial_elementIn,
961  nDOF_test_elementIn,
962  nQuadraturePoints_elementBoundaryIn,
963  CompKernelFlag);
964  }
965 }//proteus
966 
967 #endif
proteus::MoveMesh2D::evaluateCoefficients
void evaluateCoefficients(double det_J, const double *materialProperties, double *strain, double *stress, double *dstress)
Definition: MoveMesh2D.h:83
proteus::MoveMesh2D::exteriorNumericalStressFlux
void exteriorNumericalStressFlux(const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isStressFluxBoundary_u, const int &isStressFluxBoundary_v, const double &penalty, const double &u, const double &v, const double &bc_u, const double &bc_v, const double &bc_stressFlux_u, const double &bc_stressFlux_v, const double *stress, const double *normal, double &stressFlux_u, double &stressFlux_v)
Definition: MoveMesh2D.h:125
proteus::MoveMesh2D::ck
CompKernelType ck
Definition: MoveMesh2D.h:32
proteus::MoveMesh2D_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::MoveMesh2D::XX
const int XX
Definition: MoveMesh2D.h:39
proteus::MoveMesh2D::YX
const int YX
Definition: MoveMesh2D.h:40
proteus::MoveMesh2D::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: MoveMesh2D.h:34
proteus::MoveMesh2D
Definition: MoveMesh2D.h:30
proteus::MoveMesh2D::sXX
const int sXX
Definition: MoveMesh2D.h:42
proteus::MoveMesh2D::calculateStrain
void calculateStrain(double *D, double *strain)
Definition: MoveMesh2D.h:72
D
#define D(x, y)
Definition: jf.h:9
proteus::MoveMesh2D_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::MoveMesh2D::HXHX
const int HXHX
Definition: MoveMesh2D.h:50
proteus::MoveMesh2D::X
const int X
Definition: MoveMesh2D.h:38
proteus::newMoveMesh2D
MoveMesh2D_base * newMoveMesh2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: MoveMesh2D.h:949
proteus::MoveMesh2D::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)
Definition: MoveMesh2D.h:594
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::MoveMesh2D::XHX
const int XHX
Definition: MoveMesh2D.h:46
proteus::MoveMesh2D::sYY
const int sYY
Definition: MoveMesh2D.h:43
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::MoveMesh2D::Y
const int Y
Definition: MoveMesh2D.h:38
v
Double v
Definition: Headers.h:95
proteus::MoveMesh2D::exteriorNumericalStressFluxJacobian
void exteriorNumericalStressFluxJacobian(const int &isDOFBoundary_u, const int &isDOFBoundary_v, 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_v_u, double &dstressFlux_v_v)
Definition: MoveMesh2D.h:187
proteus::MoveMesh2D::MoveMesh2D
MoveMesh2D()
Definition: MoveMesh2D.h:52
z
Double * z
Definition: Headers.h:49
u
Double u
Definition: Headers.h:89
proteus::MoveMesh2D_base
Definition: MoveMesh2D.h:15
CompKernel
Definition: CompKernel.h:1018
proteus::MoveMesh2D::XY
const int XY
Definition: MoveMesh2D.h:39
proteus
Definition: ADR.h:17
proteus::MoveMesh2D::ei
EIndex< nSpace > ei
Definition: MoveMesh2D.h:36
proteus::MoveMesh2D::YY
const int YY
Definition: MoveMesh2D.h:40
proteus::MoveMesh2D::sXY
const int sXY
Definition: MoveMesh2D.h:42
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::MoveMesh2D_base::~MoveMesh2D_base
virtual ~MoveMesh2D_base()
Definition: MoveMesh2D.h:17
proteus::MoveMesh2D::nSymTen
const int nSymTen
Definition: MoveMesh2D.h:45
proteus::MoveMesh2D::sYX
const int sYX
Definition: MoveMesh2D.h:43
proteus::MoveMesh2D::YHX
const int YHX
Definition: MoveMesh2D.h:47
ArgumentsDict.h
EIndex< nSpace >
proteus::MoveMesh2D::calculateResidual
virtual void calculateResidual(arguments_dict &args)
Definition: MoveMesh2D.h:280
Generated on Fri Jul 1 2022 11:17:39 for proteus by doxygen 1.8.20