proteus  1.8.1
C/C++/Fortran libraries
MoveMesh.h
Go to the documentation of this file.
1 #ifndef MOVEMESH_H
2 #define MOVEMESH_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 MoveMesh_base
15  {
16  public:
17  virtual ~MoveMesh_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 MoveMesh : public MoveMesh_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  MoveMesh():
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[3],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[3],
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 
948 
949  inline void Invert3by3(double* Amat,double* Ainv)
950  {
951 
952  Ainv[0*3+0] = Amat[1*3+1] * Amat[2*3+2] - Amat[2*3+1] * Amat[1*3+2];
953  Ainv[0*3+1] = Amat[2*3+1] * Amat[0*3+2] - Amat[0*3+1] * Amat[2*3+2];
954  Ainv[0*3+2] = Amat[0*3+1] * Amat[1*3+2] - Amat[0*3+2] * Amat[1*3+1];
955 
956  double tmp = 1.0 / ( Ainv[0*3+0] * Amat[0*3+0]
957  + Ainv[0*3+1] * Amat[1*3+0]
958  + Ainv[0*3+2] * Amat[2*3+0] );
959  Ainv[0*3+0] = Ainv[0*3+0] * tmp;
960  Ainv[0*3+1] = Ainv[0*3+1] * tmp;
961  Ainv[0*3+2] = Ainv[0*3+2] * tmp;
962 
963  Ainv[1*3+0] = (Amat[1*3+2] * Amat[2*3+0] - Amat[1*3+0] * Amat[2*3+2]) * tmp;
964  Ainv[1*3+1] = (Amat[0*3+0] * Amat[2*3+2] - Amat[2*3+0] * Amat[0*3+2]) * tmp;
965  Ainv[1*3+2] = (Amat[1*3+0] * Amat[0*3+2] - Amat[0*3+0] * Amat[1*3+2]) * tmp;
966 
967  Ainv[2*3+0] = (Amat[1*3+0] * Amat[2*3+1] - Amat[1*3+1] * Amat[2*3+0]) * tmp;
968  Ainv[2*3+1] = (Amat[2*3+0] * Amat[0*3+1] - Amat[0*3+0] * Amat[2*3+1]) * tmp;
969  Ainv[2*3+2] = (Amat[0*3+0] * Amat[1*3+1] - Amat[0*3+1] * Amat[1*3+0]) * tmp;
970 
971  }
972 
973  /* //---------------------------------------------------------------------------- */
974  /* virtual void moveRigidBody(double mass, double* inertiaRef, */
975  /* double* force, double* moment, */
976  /* double* disp0, double* disp1, */
977  /* double* vel0, double* vel1, */
978  /* double* rot0, double* rot1, */
979  /* double* angVel0, double* angVel1, */
980  /* double deltaT, */
981  /* int* linConstraints, int* angConstraints, */
982  /* double linRelaxFac, double angRelaxFac, */
983  /* double linNorm, double angNorm, */
984  /* int iterMax) */
985  /* { */
986  /* int i,j,k,l, iter; */
987  /* double inertia0[3*3], inertia1[3*3], res[3], res2[3*3], Ainv[3*3], angVelTen[3*3],RotLHS[3*3]; */
988 
989  /* //std::cout<<" lin const = "<<linConstraints[0]<<" "<<linConstraints[1]<<" "<<linConstraints[2]<<std::endl; */
990  /* //std::cout<<" ang const = "<<angConstraints[0]<<" "<<angConstraints[1]<<" "<<angConstraints[2]<<std::endl; */
991 
992  /* //----------------------------------------------- */
993  /* // Linear balance */
994  /* // Integrated using Crank-Nicolson (midpoint) */
995  /* //----------------------------------------------- */
996  /* for (i = 0; i < 3; i++) */
997  /* { */
998  /* if(linConstraints[i] >= 1) */
999  /* { */
1000  /* vel1[i] = 0.0; */
1001  /* disp1[i] = disp0[i]; */
1002  /* res[i] = 0.0; */
1003  /* } */
1004  /* else */
1005  /* { */
1006  /* res[i] = force[i] - (mass/deltaT)*(vel1[i]-vel0[i]); */
1007  /* vel1[i] += linRelaxFac*(deltaT/mass)*res[i]; */
1008 
1009  /* disp1[i] = disp0[i] + 0.5*deltaT*(vel1[i] + vel0[i]); */
1010  /* } */
1011  /* } */
1012 
1013  /* linNorm = 0.0; */
1014  /* for (i = 0; i < 3; i++) */
1015  /* { */
1016  /* linNorm += res[i]*res[i]; */
1017  /* } */
1018  /* linNorm = sqrt(linNorm); */
1019 
1020  /* //--------------------------------------- */
1021  /* // Get inertia tensor - at time level 0 */
1022  /* //--------------------------------------- */
1023  /* for (i = 0; i < 3; i++) */
1024  /* { */
1025  /* for (j = 0; j < 3; j++) */
1026  /* { */
1027  /* inertia0[i*3+j] = 0.0; */
1028  /* for (k = 0; k < 3; k++) */
1029  /* { */
1030  /* for (l = 0; l < 3; l++) */
1031  /* { */
1032  /* inertia0[i*3+j] += rot0[i*3+k]*inertiaRef[k*3+l]*rot0[j*3+l]; */
1033  /* } */
1034  /* } */
1035  /* } */
1036  /* } */
1037 
1038  /* //--------------------------------------- */
1039  /* for (iter = 0; iter<iterMax; iter++) */
1040  /* { */
1041  /* //--------------------------------------- */
1042  /* // Get inertia tensor - at time level 1 */
1043  /* //--------------------------------------- */
1044  /* for (i = 0; i < 3; i++) */
1045  /* { */
1046  /* for (j = 0; j < 3; j++) */
1047  /* { */
1048  /* inertia1[i*3+j] = 0.0; */
1049  /* for (k = 0; k < 3; k++) */
1050  /* { */
1051  /* for (l = 0; l < 3; l++) */
1052  /* { */
1053  /* inertia1[i*3+j] += rot1[i*3+k]*inertiaRef[k*3+l]*rot1[j*3+l]; */
1054  /* } */
1055  /* } */
1056  /* } */
1057  /* } */
1058  /* //--------------------------------------- */
1059  /* // Get angular residuals */
1060  /* //--------------------------------------- */
1061  /* for (i = 0; i < 3; i++) */
1062  /* { */
1063  /* if(angConstraints[i] >= 1) */
1064  /* { */
1065  /* res[i] =0.0; */
1066 
1067  /* for (j = 0; j < 3; j++) */
1068  /* { */
1069  /* inertia1[i*3+j] = 0.0; */
1070  /* inertia1[j*3+i] = 0.0; */
1071  /* } */
1072  /* inertia1[i*3+i] = 1.0; */
1073  /* } */
1074  /* else */
1075  /* { */
1076  /* res[i] = moment[i]; */
1077  /* for (j = 0; j < 3; j++) */
1078  /* { */
1079  /* res[i] -= (inertia1[i*3+j]*angVel1[j] - inertia0[i*3+j]*angVel0[j])/deltaT; */
1080  /* } */
1081  /* } */
1082  /* } */
1083 
1084  /* angNorm = 0.0; */
1085  /* for (i = 0; i < 3; i++) */
1086  /* { */
1087  /* angNorm += res[i]*res[i]; */
1088  /* } */
1089  /* angNorm = sqrt(angNorm); */
1090 
1091  /* //--------------------------------------- */
1092  /* // Compute angular velocity */
1093  /* //--------------------------------------- */
1094  /* Invert3by3(inertia1,Ainv); */
1095 
1096  /* for (i = 0; i < 3; i++) */
1097  /* { */
1098  /* for (j = 0; j < 3; j++) */
1099  /* { */
1100  /* angVel1[i] += (angRelaxFac*deltaT)*Ainv[i*3+j]*res[j]; */
1101  /* } */
1102  /* } */
1103 
1104  /* //------------------------------------------------- */
1105  /* // Compute rotation matrix */
1106  /* // */
1107  /* // Integrating using Crank-Nicolson (midpoint) */
1108  /* // According to Hughes-Winget the rotation matrix */
1109  /* // remains orthonormal (a proper rotation matrix) */
1110  /* // */
1111  /* //------------------------------------------------- */
1112  /* angVelTen[0*3+0] = 0.0; */
1113  /* angVelTen[0*3+1] = -0.5*(angVel0[2]+angVel1[2]); */
1114  /* angVelTen[0*3+2] = 0.5*(angVel0[1]+angVel1[1]); */
1115 
1116  /* angVelTen[1*3+0] = 0.5*(angVel0[2]+angVel1[2]); */
1117  /* angVelTen[1*3+1] = 0.0; */
1118  /* angVelTen[1*3+2] = -0.5*(angVel0[0]+angVel1[0]); */
1119 
1120  /* angVelTen[2*3+0] = -0.5*(angVel0[1]+angVel1[1]); */
1121  /* angVelTen[2*3+1] = 0.5*(angVel0[0]+angVel1[0]); */
1122  /* angVelTen[2*3+2] = 0.0; */
1123 
1124  /* for (i = 0; i < 3; i++) */
1125  /* { */
1126  /* for (j = 0; j < 3; j++) */
1127  /* { */
1128  /* res2[i*3+j] = - (rot1[i*3+j]-rot0[i*3+j])/deltaT; */
1129 
1130  /* for (k = 0; k < 3; k++) */
1131  /* { */
1132  /* res2[i*3+j] += angVelTen[i*3+k]*0.5*(rot1[k*3+j]+rot0[k*3+j]); */
1133  /* } */
1134  /* } */
1135  /* } */
1136 
1137  /* for (i = 0; i < 3; i++) */
1138  /* { */
1139  /* for (j = 0; j < 3; j++) */
1140  /* { */
1141  /* RotLHS[i*3+j] = -0.5*angVelTen[i*3+j]; */
1142  /* } */
1143  /* RotLHS[i*3+i] += 1.0/deltaT; */
1144  /* } */
1145 
1146  /* Invert3by3(RotLHS,Ainv); */
1147 
1148  /* for (i = 0; i < 3; i++) */
1149  /* { */
1150  /* for (j = 0; j < 3; j++) */
1151  /* { */
1152  /* for (k = 0; k < 3; k++) */
1153  /* { */
1154  /* rot1[i*3+j] += Ainv[i*3+k]*res2[k*3+j]; */
1155  /* } */
1156  /* } */
1157  /* } */
1158 
1159  /* } // End of newton loop */
1160 
1161  /* //std::cout<<" lin const = "<<linConstraints[0]<<" "<<linConstraints[1]<<" "<<linConstraints[2]<<std::endl; */
1162  /* //std::cout<<" ang const = "<<angConstraints[0]<<" "<<angConstraints[1]<<" "<<angConstraints[2]<<std::endl; */
1163  /* } */
1164 
1165  }; //MoveMesh
1166 
1167  inline MoveMesh_base* newMoveMesh(int nSpaceIn,
1168  int nQuadraturePoints_elementIn,
1169  int nDOF_mesh_trial_elementIn,
1170  int nDOF_trial_elementIn,
1171  int nDOF_test_elementIn,
1172  int nQuadraturePoints_elementBoundaryIn,
1173  int CompKernelFlag)
1174  {
1175  return proteus::chooseAndAllocateDiscretization<MoveMesh_base,MoveMesh,CompKernel>(nSpaceIn,
1176  nQuadraturePoints_elementIn,
1177  nDOF_mesh_trial_elementIn,
1178  nDOF_trial_elementIn,
1179  nDOF_test_elementIn,
1180  nQuadraturePoints_elementBoundaryIn,
1181  CompKernelFlag);
1182  }
1183 }//proteus
1184 
1185 #endif
proteus::MoveMesh::Invert3by3
void Invert3by3(double *Amat, double *Ainv)
Definition: MoveMesh.h:949
proteus::MoveMesh::ZHY
const int ZHY
Definition: MoveMesh.h:48
proteus::MoveMesh::ZY
const int ZY
Definition: MoveMesh.h:41
proteus::MoveMesh::calculateStrain
void calculateStrain(double *D, double *strain)
Definition: MoveMesh.h:72
proteus::MoveMesh::YZ
const int YZ
Definition: MoveMesh.h:40
w
#define w(x)
Definition: jf.h:22
proteus::MoveMesh::YX
const int YX
Definition: MoveMesh.h:40
proteus::MoveMesh::evaluateCoefficients
void evaluateCoefficients(double det_J, const double *materialProperties, double *strain, double *stress, double *dstress)
Definition: MoveMesh.h:83
proteus::MoveMesh::HYHY
const int HYHY
Definition: MoveMesh.h:50
D
#define D(x, y)
Definition: jf.h:9
proteus::MoveMesh::XHX
const int XHX
Definition: MoveMesh.h:46
proteus::MoveMesh::sYX
const int sYX
Definition: MoveMesh.h:43
proteus::MoveMesh
Definition: MoveMesh.h:30
proteus::MoveMesh::HXHY
const int HXHY
Definition: MoveMesh.h:49
proteus::MoveMesh::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)
Definition: MoveMesh.h:594
proteus::MoveMesh_base::~MoveMesh_base
virtual ~MoveMesh_base()
Definition: MoveMesh.h:17
proteus::MoveMesh::X
const int X
Definition: MoveMesh.h:38
proteus::MoveMesh::HXHX
const int HXHX
Definition: MoveMesh.h:49
ModelFactory.h
CompKernel.h
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::MoveMesh::XX
const int XX
Definition: MoveMesh.h:39
proteus::MoveMesh::XZ
const int XZ
Definition: MoveMesh.h:39
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
proteus::MoveMesh::Z
const int Z
Definition: MoveMesh.h:38
proteus::MoveMesh_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::MoveMesh::ZHX
const int ZHX
Definition: MoveMesh.h:48
v
Double v
Definition: Headers.h:95
proteus::MoveMesh::XHY
const int XHY
Definition: MoveMesh.h:46
proteus::MoveMesh::nSymTen
const int nSymTen
Definition: MoveMesh.h:45
proteus::MoveMesh::sYY
const int sYY
Definition: MoveMesh.h:43
proteus::MoveMesh::YHY
const int YHY
Definition: MoveMesh.h:47
z
Double * z
Definition: Headers.h:49
proteus::MoveMesh::HYHX
const int HYHX
Definition: MoveMesh.h:50
u
Double u
Definition: Headers.h:89
proteus::MoveMesh::ei
EIndex< nSpace > ei
Definition: MoveMesh.h:36
proteus::MoveMesh::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition: MoveMesh.h:34
proteus::MoveMesh::ZZ
const int ZZ
Definition: MoveMesh.h:41
proteus::MoveMesh::YY
const int YY
Definition: MoveMesh.h:40
CompKernel
Definition: CompKernel.h:1018
proteus::MoveMesh::sZX
const int sZX
Definition: MoveMesh.h:44
proteus::MoveMesh::XY
const int XY
Definition: MoveMesh.h:39
proteus::MoveMesh_base
Definition: MoveMesh.h:15
proteus::MoveMesh::YHX
const int YHX
Definition: MoveMesh.h:47
proteus
Definition: ADR.h:17
proteus::MoveMesh::sXY
const int sXY
Definition: MoveMesh.h:42
proteus::newMoveMesh
MoveMesh_base * newMoveMesh(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition: MoveMesh.h:1167
proteus::MoveMesh::sXX
const int sXX
Definition: MoveMesh.h:42
proteus::MoveMesh::sXZ
const int sXZ
Definition: MoveMesh.h:42
proteus::arguments_dict
Definition: ArgumentsDict.h:70
proteus::MoveMesh::sYZ
const int sYZ
Definition: MoveMesh.h:43
proteus::MoveMesh_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::MoveMesh::calculateResidual
virtual void calculateResidual(arguments_dict &args)
Definition: MoveMesh.h:280
proteus::MoveMesh::sZY
const int sZY
Definition: MoveMesh.h:44
proteus::MoveMesh::exteriorNumericalStressFlux
void exteriorNumericalStressFlux(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: MoveMesh.h:125
proteus::MoveMesh::ZX
const int ZX
Definition: MoveMesh.h:41
proteus::MoveMesh::sZZ
const int sZZ
Definition: MoveMesh.h:44
proteus::MoveMesh::MoveMesh
MoveMesh()
Definition: MoveMesh.h:52
ArgumentsDict.h
proteus::MoveMesh::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: MoveMesh.h:187
EIndex< nSpace >
proteus::MoveMesh::Y
const int Y
Definition: MoveMesh.h:38
proteus::MoveMesh::ck
CompKernelType ck
Definition: MoveMesh.h:32
Generated on Fri Jul 1 2022 11:17:39 for proteus by doxygen 1.8.20