proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
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
10namespace py = pybind11;
11
12namespace 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 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 int //eN_k = eN*nQuadraturePoints_element+k,
355 //eN_k_nSpace=eN_k*nSpace,
356 eN_nDOF_trial_element = eN*nDOF_trial_element;
357 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 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 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 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 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 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 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 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 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 int i_nSpace = i*nSpace;
757 for(int j=0;j<nDOF_trial_element;j++)
758 {
759 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 int eN_i = eN*nDOF_test_element+i;
781 for (int j=0;j<nDOF_trial_element;j++)
782 {
783 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 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 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 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 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 int eN_i = eN*nDOF_test_element+i;
927 for (int j=0;j<nDOF_trial_element;j++)
928 {
929 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
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
proteus::MoveMesh_base
Definition MoveMesh.h:15
proteus::MoveMesh_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::MoveMesh_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::MoveMesh_base::~MoveMesh_base
virtual ~MoveMesh_base()
Definition MoveMesh.h:17
proteus::MoveMesh::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition MoveMesh.h:34
proteus::MoveMesh::HXHX
const int HXHX
Definition MoveMesh.h:49
proteus::MoveMesh::sXY
const int sXY
Definition MoveMesh.h:42
proteus::MoveMesh::sZY
const int sZY
Definition MoveMesh.h:44
proteus::MoveMesh::ZHX
const int ZHX
Definition MoveMesh.h:48
proteus::MoveMesh::Invert3by3
void Invert3by3(double *Amat, double *Ainv)
Definition MoveMesh.h:949
proteus::MoveMesh::XHY
const int XHY
Definition MoveMesh.h:46
proteus::MoveMesh::YHY
const int YHY
Definition MoveMesh.h:47
proteus::MoveMesh::nSymTen
const int nSymTen
Definition MoveMesh.h:45
proteus::MoveMesh::YZ
const int YZ
Definition MoveMesh.h:40
proteus::MoveMesh::sXX
const int sXX
Definition MoveMesh.h:42
proteus::MoveMesh::ZY
const int ZY
Definition MoveMesh.h:41
proteus::MoveMesh::Z
const int Z
Definition MoveMesh.h:38
proteus::MoveMesh::ZX
const int ZX
Definition MoveMesh.h:41
proteus::MoveMesh::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)
Definition MoveMesh.h:594
proteus::MoveMesh::HXHY
const int HXHY
Definition MoveMesh.h:49
proteus::MoveMesh::XX
const int XX
Definition MoveMesh.h:39
proteus::MoveMesh::ck
CompKernelType ck
Definition MoveMesh.h:32
proteus::MoveMesh::XZ
const int XZ
Definition MoveMesh.h:39
proteus::MoveMesh::sYY
const int sYY
Definition MoveMesh.h:43
proteus::MoveMesh::calculateStrain
void calculateStrain(double *D, double *strain)
Definition MoveMesh.h:72
proteus::MoveMesh::YY
const int YY
Definition MoveMesh.h:40
proteus::MoveMesh::XY
const int XY
Definition MoveMesh.h:39
proteus::MoveMesh::ZZ
const int ZZ
Definition MoveMesh.h:41
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::YX
const int YX
Definition MoveMesh.h:40
proteus::MoveMesh::HYHX
const int HYHX
Definition MoveMesh.h:50
proteus::MoveMesh::sXZ
const int sXZ
Definition MoveMesh.h:42
proteus::MoveMesh::evaluateCoefficients
void evaluateCoefficients(double det_J, const double *materialProperties, double *strain, double *stress, double *dstress)
Definition MoveMesh.h:83
proteus::MoveMesh::sYX
const int sYX
Definition MoveMesh.h:43
proteus::MoveMesh::HYHY
const int HYHY
Definition MoveMesh.h:50
proteus::MoveMesh::calculateResidual
virtual void calculateResidual(arguments_dict &args)
Definition MoveMesh.h:280
proteus::MoveMesh::ZHY
const int ZHY
Definition MoveMesh.h:48
proteus::MoveMesh::XHX
const int XHX
Definition MoveMesh.h:46
proteus::MoveMesh::X
const int X
Definition MoveMesh.h:38
proteus::MoveMesh::sZX
const int sZX
Definition MoveMesh.h:44
proteus::MoveMesh::Y
const int Y
Definition MoveMesh.h:38
proteus::MoveMesh::ei
EIndex< nSpace > ei
Definition MoveMesh.h:36
proteus::MoveMesh::YHX
const int YHX
Definition MoveMesh.h:47
proteus::MoveMesh::sZZ
const int sZZ
Definition MoveMesh.h:44
proteus::MoveMesh::sYZ
const int sYZ
Definition MoveMesh.h:43
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
D
#define D(x, y)
Definition jf.h:9
w
#define w(x)
Definition jf.h:22
proteus
Definition ADR.h:19
proteus::chooseAndAllocateDiscretization
Model_Base * chooseAndAllocateDiscretization(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nDOF_v_trial_elementIn, int nDOF_v_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MixedModelFactory.h:48
proteus::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::arguments_dict
Definition ArgumentsDict.h:70
proteus::arguments_dict::scalar
T & scalar(const std::string &key)
proteus::arguments_dict::array
xt::pyarray< T > & array(const std::string &key)
Generated on for proteus by doxygen 1.18.0