proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
MoveMesh2D.h
Go to the documentation of this file.
1#ifndef MOVEMESH2D_H
2#define MOVEMESH2D_H
3#include <cmath>
4#include <iostream>
5#include "CompKernel.h"
6#include "ModelFactory.h"
7#include "ArgumentsDict.h"
8#include "xtensor-python/pyarray.hpp"
9
10namespace py = pybind11;
11
12namespace proteus
13{
14 class MoveMesh2D_base
15 {
16 public:
17 virtual ~MoveMesh2D_base(){}
18 virtual void calculateResidual(arguments_dict& args)=0;
19 virtual void calculateJacobian(arguments_dict& args)=0;
20 };
21
22 template<class CompKernelType,
23 int nSpace,
24 int nQuadraturePoints_element,
25 int nDOF_mesh_trial_element,
26 int nDOF_trial_element,
27 int nDOF_test_element,
28 int nQuadraturePoints_elementBoundary>
29 class MoveMesh2D : public MoveMesh2D_base
30 {
31 public:
32 CompKernelType ck;
33
34 const int nDOF_test_X_trial_element;
35
36 EIndex<nSpace> ei;
37
38 const int X,Y,/* Z, */
39 XX,XY,/* XZ, */
40 YX,YY,/* YZ, */
41 /* ZX,ZY,ZZ, */
42 sXX,sXY,/* sXZ, */
43 sYX,sYY,/* sYZ, */
44 /* sZX,sZY,sZZ, */
45 nSymTen,
46 XHX,/* XHY, */
47 YHX,/* YHY, */
48 /* ZHX,ZHY, */
49 HXHX/*, HXHY, */
50 /* HYHX,HYHY */;
51
52 MoveMesh2D():
53 ck(),
54 nDOF_test_X_trial_element(nDOF_test_element*nDOF_trial_element),
55 X(ei.X),
56 Y(ei.Y),
57 /* Z(ei.Z), */
58 XX(ei.XX),XY(ei.XY),/* XZ(ei.XZ), */
59 YX(ei.YX),YY(ei.YY),/* YZ(ei.YZ), */
60 /* ZX(ei.ZX),ZY(ei.ZY),ZZ(ei.ZZ), */
61 sXX(ei.sXX),sXY(ei.sXY),/* sXZ(ei.sXZ), */
62 sYX(ei.sYX),sYY(ei.sYY),/* sYZ(ei.sYZ), */
63 /* sZX(ei.sZX),sZY(ei.sZY),sZZ(ei.sZZ), */
64 nSymTen(ei.nSymTen),
65 XHX(ei.XHX),/* XHY(ei.XHY), */
66 YHX(ei.YHX),/* YHY(ei.YHY), */
67 /* ZHX(ei.ZHX),ZHY(ei.ZHY), */
68 HXHX(ei.HXHX)/* ,HXHY(ei.HXHY), */
69 /* HYHX(ei.HYHX),HYHY(ei.HYHY) */
70 {}
71
72 inline void calculateStrain(double* D, double* strain)
73 {
74 //Voigt notation from Belytschko, Liu, Moran
75 strain[sXX] = D[XX];//du/dx
76 strain[sYY] = D[YY];//dv/dy
77 /* strain[sZZ] = D[ZZ];//dw/dz */
78 /* strain[sYZ] = D[YZ]+D[ZY];//(dv/dz + dz/dy) */
79 /* strain[sXZ] = D[XZ]+D[ZX];//(du/dz + dw/dx) */
80 strain[sXY] = D[XY]+D[YX];//(du/dy + dv/dx)
81 }
82
83 inline void evaluateCoefficients(double det_J,
84 const double* materialProperties,
85 double* strain,
86 double* stress,
87 double* dstress)
88 {
89 //cek hack/todo need to set E based on reference configuration
90 const double strainTrace=(strain[sXX]+strain[sYY]/* +strain[sZZ] */),
91 E=materialProperties[0]/det_J,//for mesh motion penalize small elements
92 nu=materialProperties[1];
93
94 const double shear = E/(1.0+nu);
95 const double bulk = shear*(nu/(1.0-2.0*nu));
96
97 for (int i=0;i<nSymTen;i++)
98 for (int j=0;j<nSymTen;j++)
99 dstress[i*nSymTen+j] = 0.0;
100 stress[sXX] = shear*strain[sXX] + bulk*strainTrace;
101 dstress[sXX*nSymTen+sXX] = bulk + shear;
102 dstress[sXX*nSymTen+sYY] = bulk;
103 /* dstress[sXX*nSymTen+sZZ] = bulk; */
104
105 stress[sYY] = shear*strain[sYY] + bulk*strainTrace;
106 dstress[sYY*nSymTen+sXX] = bulk;
107 dstress[sYY*nSymTen+sYY] = bulk + shear;
108 /* dstress[sYY*nSymTen+sZZ] = bulk; */
109
110 /* stress[sZZ] = shear*strain[sZZ] + bulk*strainTrace; */
111 /* dstress[sZZ*nSymTen+sXX] = bulk; */
112 /* dstress[sZZ*nSymTen+sYY] = bulk; */
113 /* dstress[sZZ*nSymTen+sZZ] = bulk + shear; */
114
115 /* stress[sYZ] = shear*0.5*strain[sYZ];//the 1/2 comes from the Voigt notation */
116 /* dstress[sYZ*nSymTen+sYZ] = shear*0.5; */
117
118 /* stress[sXZ] = shear*0.5*strain[sXZ];//the 1/2 comes from the Voigt notation */
119 /* dstress[sXZ*nSymTen+sXZ] = shear*0.5; */
120
121 stress[sXY] = shear*0.5*strain[sXY];//the 1/2 comes from the Voigt notation
122 dstress[sXY*nSymTen+sXY] = shear*0.5;
123 }
124
125 inline void exteriorNumericalStressFlux(const int& isDOFBoundary_u,
126 const int& isDOFBoundary_v,
127 /* const int& isDOFBoundary_w, */
128 const int& isStressFluxBoundary_u,
129 const int& isStressFluxBoundary_v,
130 /* const int& isStressFluxBoundary_w, */
131 const double& penalty,
132 const double& u,
133 const double& v,
134 /* const double& w, */
135 const double& bc_u,
136 const double& bc_v,
137 /* const double& bc_w, */
138 const double& bc_stressFlux_u,
139 const double& bc_stressFlux_v,
140 /* const double& bc_stressFlux_w, */
141 const double* stress,
142 const double* normal,
143 double& stressFlux_u,
144 double& stressFlux_v/* , */
145 /* double& stressFlux_w */)
146 {
147 if (isDOFBoundary_u == 1)
148 {
149 stressFlux_u = -(stress[sXX]*normal[X] + stress[sXY]*normal[Y] /* + stress[sXZ]*normal[Z] */ - penalty*(u - bc_u));
150 }
151 else if(isStressFluxBoundary_u == 1)
152 {
153 stressFlux_u = bc_stressFlux_u;
154 }
155 else
156 {
157 stressFlux_u = 0.0;
158 }
159
160 if (isDOFBoundary_v == 1)
161 {
162 stressFlux_v = -(stress[sYX]*normal[X] + stress[sYY]*normal[Y] /* + stress[sYZ]*normal[Z] */ - penalty*(v - bc_v));
163 }
164 else if(isStressFluxBoundary_v == 1)
165 {
166 stressFlux_v = bc_stressFlux_v;
167 }
168 else
169 {
170 stressFlux_v = 0.0;
171 }
172
173 /* if (isDOFBoundary_w == 1) */
174 /* { */
175 /* stressFlux_w = -(stress[sZX]*normal[X] + stress[sZY]*normal[Y] + stress[sZZ]*normal[Z] - penalty*(w - bc_w)); */
176 /* } */
177 /* else if(isStressFluxBoundary_w == 1) */
178 /* { */
179 /* stressFlux_w = bc_stressFlux_w; */
180 /* } */
181 /* else */
182 /* { */
183 /* stressFlux_w = 0.0; */
184 /* } */
185 }
186
187 inline void exteriorNumericalStressFluxJacobian(const int& isDOFBoundary_u,
188 const int& isDOFBoundary_v,
189 /* const int& isDOFBoundary_w, */
190 const double* normal,
191 const double* dstress,
192 const double& penalty,
193 const double& disp_trial,
194 const double* disp_grad_trial,
195 double& dstressFlux_u_u,
196 double& dstressFlux_u_v,
197 /* double& dstressFlux_u_w, */
198 double& dstressFlux_v_u,
199 double& dstressFlux_v_v/* , */
200 /* double& dstressFlux_v_w, */
201 /* double& dstressFlux_w_u, */
202 /* double& dstressFlux_w_v, */
203 /* double& dstressFlux_w_w */)
204 {
205 //here we use both the symmetry of the stress tensor and the fact that dstress is w.r.t. the strain in Voigt notation to go directly to derivatives w.r.t. displacement DOF
206 if (isDOFBoundary_u == 1)
207 {
208 dstressFlux_u_u = -(
209 (dstress[sXX*nSymTen+sXX]*disp_grad_trial[X] + dstress[sXX*nSymTen+sXY]*disp_grad_trial[Y] /* + dstress[sXX*nSymTen+sXZ]*disp_grad_trial[Z] */)*normal[X]+
210 (dstress[sXY*nSymTen+sXX]*disp_grad_trial[X] + dstress[sXY*nSymTen+sXY]*disp_grad_trial[Y] /* + dstress[sXY*nSymTen+sXZ]*disp_grad_trial[Z] */)*normal[Y]+
211 /* (dstress[sXZ*nSymTen+sXX]*disp_grad_trial[X] + dstress[sXZ*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sXZ*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Z] */
212 -
213 penalty*disp_trial);
214 dstressFlux_u_v = -(
215 (dstress[sXX*nSymTen+sYX]*disp_grad_trial[X] + dstress[sXX*nSymTen+sYY]*disp_grad_trial[Y] /* + dstress[sXX*nSymTen+sYZ]*disp_grad_trial[Z] */)*normal[X]+
216 (dstress[sXY*nSymTen+sYX]*disp_grad_trial[X] + dstress[sXY*nSymTen+sYY]*disp_grad_trial[Y] /* + dstress[sXY*nSymTen+sYZ]*disp_grad_trial[Z] */)*normal[Y]/* + */
217 /* (dstress[sXZ*nSymTen+sYX]*disp_grad_trial[X] + dstress[sXZ*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sXZ*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Z] */);
218 /* dstressFlux_u_w = -( */
219 /* (dstress[sXX*nSymTen+sZX]*disp_grad_trial[X] + dstress[sXX*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sXX*nSymTen+sZZ]*disp_grad_trial[Z])*normal[X]+ */
220 /* (dstress[sXY*nSymTen+sZX]*disp_grad_trial[X] + dstress[sXY*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sXY*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Y]+ */
221 /* (dstress[sXZ*nSymTen+sZX]*disp_grad_trial[X] + dstress[sXZ*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sXZ*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Z]); */
222 }
223 else
224 {
225 dstressFlux_u_u = 0.0;
226 dstressFlux_u_v = 0.0;
227 /* dstressFlux_u_w = 0.0; */
228 }
229
230 if (isDOFBoundary_v == 1)
231 {
232 dstressFlux_v_u = -(
233 (dstress[sYX*nSymTen+sXX]*disp_grad_trial[X] + dstress[sYX*nSymTen+sXY]*disp_grad_trial[Y] /* + dstress[sYX*nSymTen+sXZ]*disp_grad_trial[Z] */)*normal[X]+
234 (dstress[sYY*nSymTen+sXX]*disp_grad_trial[X] + dstress[sYY*nSymTen+sXY]*disp_grad_trial[Y] /* + dstress[sYY*nSymTen+sXZ]*disp_grad_trial[Z] */)*normal[Y]/* + */
235 /* (dstress[sYZ*nSymTen+sXX]*disp_grad_trial[X] + dstress[sYZ*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sYZ*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Z] */);
236 dstressFlux_v_v = -(
237 (dstress[sYX*nSymTen+sYX]*disp_grad_trial[X] + dstress[sYX*nSymTen+sYY]*disp_grad_trial[Y] /* + dstress[sYX*nSymTen+sYZ]*disp_grad_trial[Z] */)*normal[X]+
238 (dstress[sYY*nSymTen+sYX]*disp_grad_trial[X] + dstress[sYY*nSymTen+sYY]*disp_grad_trial[Y] /* + dstress[sYY*nSymTen+sYZ]*disp_grad_trial[Z] */)*normal[Y]/* + */
239 /* (dstress[sYZ*nSymTen+sYX]*disp_grad_trial[X] + dstress[sYZ*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sYZ*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Z] */
240 -
241 penalty*disp_trial);
242 /* dstressFlux_v_w = -( */
243 /* (dstress[sYX*nSymTen+sZX]*disp_grad_trial[X] + dstress[sYX*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sYX*nSymTen+sZZ]*disp_grad_trial[Z])*normal[X]+ */
244 /* (dstress[sYY*nSymTen+sZX]*disp_grad_trial[X] + dstress[sYY*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sYY*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Y]+ */
245 /* (dstress[sYZ*nSymTen+sZX]*disp_grad_trial[X] + dstress[sYZ*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sYZ*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Z]); */
246 }
247 else
248 {
249 dstressFlux_v_u = 0.0;
250 dstressFlux_v_v = 0.0;
251 /* dstressFlux_v_w = 0.0; */
252 }
253
254 /* if (isDOFBoundary_w == 1) */
255 /* { */
256 /* dstressFlux_w_u = -( */
257 /* (dstress[sZX*nSymTen+sXX]*disp_grad_trial[X] + dstress[sZX*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sZX*nSymTen+sXZ]*disp_grad_trial[Z])*normal[X]+ */
258 /* (dstress[sZY*nSymTen+sXX]*disp_grad_trial[X] + dstress[sZY*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sZY*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Y]+ */
259 /* (dstress[sZZ*nSymTen+sXX]*disp_grad_trial[X] + dstress[sZZ*nSymTen+sXY]*disp_grad_trial[Y] + dstress[sZZ*nSymTen+sXZ]*disp_grad_trial[Z])*normal[Z]); */
260 /* dstressFlux_w_v = -( */
261 /* (dstress[sZX*nSymTen+sYX]*disp_grad_trial[X] + dstress[sZX*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sZX*nSymTen+sYZ]*disp_grad_trial[Z])*normal[X]+ */
262 /* (dstress[sZY*nSymTen+sYX]*disp_grad_trial[X] + dstress[sZY*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sZY*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Y]+ */
263 /* (dstress[sZZ*nSymTen+sYX]*disp_grad_trial[X] + dstress[sZZ*nSymTen+sYY]*disp_grad_trial[Y] + dstress[sZZ*nSymTen+sYZ]*disp_grad_trial[Z])*normal[Z]); */
264 /* dstressFlux_w_w = -( */
265 /* (dstress[sZX*nSymTen+sZX]*disp_grad_trial[X] + dstress[sZX*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sZX*nSymTen+sZZ]*disp_grad_trial[Z])*normal[X]+ */
266 /* (dstress[sZY*nSymTen+sZX]*disp_grad_trial[X] + dstress[sZY*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sZY*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Y]+ */
267 /* (dstress[sZZ*nSymTen+sZX]*disp_grad_trial[X] + dstress[sZZ*nSymTen+sZY]*disp_grad_trial[Y] + dstress[sZZ*nSymTen+sZZ]*disp_grad_trial[Z])*normal[Z] */
268 /* - */
269 /* penalty*disp_trial); */
270 /* } */
271 /* else */
272 /* { */
273 /* dstressFlux_w_u = 0.0; */
274 /* dstressFlux_w_v = 0.0; */
275 /* dstressFlux_w_w = 0.0; */
276 /* } */
277 }
278
279
280 virtual void calculateResidual(arguments_dict& args)
281 {
282 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
283 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
284 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
285 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
286 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
287 xt::pyarray<double>& disp_trial_ref = args.array<double>("disp_trial_ref");
288 xt::pyarray<double>& disp_grad_trial_ref = args.array<double>("disp_grad_trial_ref");
289 xt::pyarray<double>& disp_test_ref = args.array<double>("disp_test_ref");
290 xt::pyarray<double>& disp_grad_test_ref = args.array<double>("disp_grad_test_ref");
291 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
292 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
293 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
294 xt::pyarray<double>& disp_trial_trace_ref = args.array<double>("disp_trial_trace_ref");
295 xt::pyarray<double>& disp_grad_trial_trace_ref = args.array<double>("disp_grad_trial_trace_ref");
296 xt::pyarray<double>& disp_test_trace_ref = args.array<double>("disp_test_trace_ref");
297 xt::pyarray<double>& disp_grad_test_trace_ref = args.array<double>("disp_grad_test_trace_ref");
298 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
299 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
300 int nElements_global = args.scalar<int>("nElements_global");
301 xt::pyarray<int>& materialTypes = args.array<int>("materialTypes");
302 int nMaterialProperties = args.scalar<int>("nMaterialProperties");
303 xt::pyarray<double>& materialProperties = args.array<double>("materialProperties");
304 xt::pyarray<int>& disp_l2g = args.array<int>("disp_l2g");
305 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
306 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
307 xt::pyarray<double>& w_dof = args.array<double>("w_dof");
308 xt::pyarray<double>& bodyForce = args.array<double>("bodyForce");
309 int offset_u = args.scalar<int>("offset_u");
310 int offset_v = args.scalar<int>("offset_v");
311 int offset_w = args.scalar<int>("offset_w");
312 int stride_u = args.scalar<int>("stride_u");
313 int stride_v = args.scalar<int>("stride_v");
314 int stride_w = args.scalar<int>("stride_w");
315 xt::pyarray<double>& globalResidual = args.array<double>("globalResidual");
316 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
317 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
318 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
319 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
320 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
321 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
322 xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
323 xt::pyarray<int>& isStressFluxBoundary_u = args.array<int>("isStressFluxBoundary_u");
324 xt::pyarray<int>& isStressFluxBoundary_v = args.array<int>("isStressFluxBoundary_v");
325 xt::pyarray<int>& isStressFluxBoundary_w = args.array<int>("isStressFluxBoundary_w");
326 xt::pyarray<double>& ebqe_bc_u_ext = args.array<double>("ebqe_bc_u_ext");
327 xt::pyarray<double>& ebqe_bc_v_ext = args.array<double>("ebqe_bc_v_ext");
328 xt::pyarray<double>& ebqe_bc_w_ext = args.array<double>("ebqe_bc_w_ext");
329 xt::pyarray<double>& ebqe_bc_stressFlux_u_ext = args.array<double>("ebqe_bc_stressFlux_u_ext");
330 xt::pyarray<double>& ebqe_bc_stressFlux_v_ext = args.array<double>("ebqe_bc_stressFlux_v_ext");
331 xt::pyarray<double>& ebqe_bc_stressFlux_w_ext = args.array<double>("ebqe_bc_stressFlux_w_ext");
332 //
333 //loop over elements to compute volume integrals and load them into element and global residual
334 //
335 for(int eN=0;eN<nElements_global;eN++)
336 {
337 //declare local storage for element residual and initialize
338 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[2],x_ext,y_ext,z_ext,
484 G[nSpace*nSpace],G_dd_G,tr_G,h_penalty,
485 strain[ck.nSymTen],stress[ck.nSymTen],dstress[ck.nSymTen*ck.nSymTen],
486 stressFlux_u,stressFlux_v,stressFlux_w;
487 //compute information about mapping from reference element to physical element
488 ck.calculateMapping_elementBoundary(eN,
489 ebN_local,
490 kb,
491 ebN_local_kb,
492 mesh_dof.data(),
493 mesh_l2g.data(),
494 mesh_trial_trace_ref.data(),
495 mesh_grad_trial_trace_ref.data(),
496 boundaryJac_ref.data(),
497 jac_ext,
498 jacDet_ext,
499 jacInv_ext,
500 boundaryJac,
501 metricTensor,
502 metricTensorDetSqrt,
503 normal_ref.data(),
504 normal,
505 x_ext,y_ext,z_ext);
506 dS = metricTensorDetSqrt*dS_ref.data()[kb];
507 //get the metric tensor
508 //cek todo use symmetry
509 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
510 //compute shape and solution information
511 //shape
512 ck.gradTrialFromRef(&disp_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,disp_grad_trial_trace);
513 //solution and gradients
514 ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
515 ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
516 /* ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext); */
517 ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_u_ext);
518 ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_v_ext);
519 /* ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_w_ext); */
520 //precalculate test function products with integration weights
521 for (int j=0;j<nDOF_trial_element;j++)
522 {
523 disp_test_dS[j] = disp_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
524 }
525 //
526 //load the boundary values
527 //
528 bc_u_ext = isDOFBoundary_u.data()[ebNE_kb]*ebqe_bc_u_ext.data()[ebNE_kb]+(1-isDOFBoundary_u.data()[ebNE_kb])*u_ext;
529 bc_v_ext = isDOFBoundary_v.data()[ebNE_kb]*ebqe_bc_v_ext.data()[ebNE_kb]+(1-isDOFBoundary_v.data()[ebNE_kb])*v_ext;
530 /* bc_w_ext = isDOFBoundary_w.data()[ebNE_kb]*ebqe_bc_w_ext.data()[ebNE_kb]+(1-isDOFBoundary_w.data()[ebNE_kb])*w_ext; */
531 //
532 //calculate the pde coefficients using the solution and the boundary values for the solution
533 //
534 calculateStrain(D,strain);
535 evaluateCoefficients(fabs(jacDet_ext),
536 &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
537 strain,
538 stress,
539 dstress);
540 //
541 //calculate the numerical fluxes
542 //
543 //cek debug
544 //ebqe_penalty_ext[ebNE_kb] = 10.0;
545 //
546 ck.calculateGScale(G,normal,h_penalty);
547 //cek hack
548 h_penalty = 100.0/h_penalty;
549 exteriorNumericalStressFlux(isDOFBoundary_u.data()[ebNE_kb],
550 isDOFBoundary_v.data()[ebNE_kb],
551 /* isDOFBoundary_w.data()[ebNE_kb], */
552 isStressFluxBoundary_u.data()[ebNE_kb],
553 isStressFluxBoundary_v.data()[ebNE_kb],
554 /* isStressFluxBoundary_w.data()[ebNE_kb], */
555 h_penalty,
556 u_ext,
557 v_ext,
558 /* w_ext, */
559 bc_u_ext,
560 bc_v_ext,
561 /* bc_w_ext, */
562 ebqe_bc_stressFlux_u_ext.data()[ebNE_kb],
563 ebqe_bc_stressFlux_v_ext.data()[ebNE_kb],
564 /* ebqe_bc_stressFlux_w_ext.data()[ebNE_kb], */
565 stress,
566 normal,
567 stressFlux_u,
568 stressFlux_v/* , */
569 /* stressFlux_w */);
570 //
571 //update residuals
572 //
573 for (int i=0;i<nDOF_test_element;i++)
574 {
575 elementResidual_u[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_u,disp_test_dS[i]);
576 elementResidual_v[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_v,disp_test_dS[i]);
577 /* elementResidual_w[i] += ck.ExteriorElementBoundaryStressFlux(stressFlux_w,disp_test_dS[i]); */
578 }//i
579 }//kb
580 //
581 //update the element and global residual storage
582 //
583 for (int i=0;i<nDOF_test_element;i++)
584 {
585 int eN_i = eN*nDOF_test_element+i;
586
587 globalResidual.data()[offset_u+stride_u*disp_l2g.data()[eN_i]]+=elementResidual_u[i];
588 globalResidual.data()[offset_v+stride_v*disp_l2g.data()[eN_i]]+=elementResidual_v[i];
589 /* globalResidual.data()[offset_w+stride_w*disp_l2g.data()[eN_i]]+=elementResidual_w[i]; */
590 }//i
591 }//ebNE
592 }
593
594 virtual void calculateJacobian(arguments_dict& args)
595 {
596 xt::pyarray<double>& mesh_trial_ref = args.array<double>("mesh_trial_ref");
597 xt::pyarray<double>& mesh_grad_trial_ref = args.array<double>("mesh_grad_trial_ref");
598 xt::pyarray<double>& mesh_dof = args.array<double>("mesh_dof");
599 xt::pyarray<int>& mesh_l2g = args.array<int>("mesh_l2g");
600 xt::pyarray<double>& dV_ref = args.array<double>("dV_ref");
601 xt::pyarray<double>& disp_trial_ref = args.array<double>("disp_trial_ref");
602 xt::pyarray<double>& disp_grad_trial_ref = args.array<double>("disp_grad_trial_ref");
603 xt::pyarray<double>& disp_test_ref = args.array<double>("disp_test_ref");
604 xt::pyarray<double>& disp_grad_test_ref = args.array<double>("disp_grad_test_ref");
605 xt::pyarray<double>& mesh_trial_trace_ref = args.array<double>("mesh_trial_trace_ref");
606 xt::pyarray<double>& mesh_grad_trial_trace_ref = args.array<double>("mesh_grad_trial_trace_ref");
607 xt::pyarray<double>& dS_ref = args.array<double>("dS_ref");
608 xt::pyarray<double>& disp_trial_trace_ref = args.array<double>("disp_trial_trace_ref");
609 xt::pyarray<double>& disp_grad_trial_trace_ref = args.array<double>("disp_grad_trial_trace_ref");
610 xt::pyarray<double>& disp_test_trace_ref = args.array<double>("disp_test_trace_ref");
611 xt::pyarray<double>& disp_grad_test_trace_ref = args.array<double>("disp_grad_test_trace_ref");
612 xt::pyarray<double>& normal_ref = args.array<double>("normal_ref");
613 xt::pyarray<double>& boundaryJac_ref = args.array<double>("boundaryJac_ref");
614 int nElements_global = args.scalar<int>("nElements_global");
615 xt::pyarray<int>& materialTypes = args.array<int>("materialTypes");
616 int nMaterialProperties = args.scalar<int>("nMaterialProperties");
617 xt::pyarray<double>& materialProperties = args.array<double>("materialProperties");
618 xt::pyarray<int>& disp_l2g = args.array<int>("disp_l2g");
619 xt::pyarray<double>& u_dof = args.array<double>("u_dof");
620 xt::pyarray<double>& v_dof = args.array<double>("v_dof");
621 xt::pyarray<double>& w_dof = args.array<double>("w_dof");
622 xt::pyarray<double>& bodyForce = args.array<double>("bodyForce");
623 xt::pyarray<int>& csrRowIndeces_u_u = args.array<int>("csrRowIndeces_u_u");
624 xt::pyarray<int>& csrColumnOffsets_u_u = args.array<int>("csrColumnOffsets_u_u");
625 xt::pyarray<int>& csrRowIndeces_u_v = args.array<int>("csrRowIndeces_u_v");
626 xt::pyarray<int>& csrColumnOffsets_u_v = args.array<int>("csrColumnOffsets_u_v");
627 xt::pyarray<int>& csrRowIndeces_u_w = args.array<int>("csrRowIndeces_u_w");
628 xt::pyarray<int>& csrColumnOffsets_u_w = args.array<int>("csrColumnOffsets_u_w");
629 xt::pyarray<int>& csrRowIndeces_v_u = args.array<int>("csrRowIndeces_v_u");
630 xt::pyarray<int>& csrColumnOffsets_v_u = args.array<int>("csrColumnOffsets_v_u");
631 xt::pyarray<int>& csrRowIndeces_v_v = args.array<int>("csrRowIndeces_v_v");
632 xt::pyarray<int>& csrColumnOffsets_v_v = args.array<int>("csrColumnOffsets_v_v");
633 xt::pyarray<int>& csrRowIndeces_v_w = args.array<int>("csrRowIndeces_v_w");
634 xt::pyarray<int>& csrColumnOffsets_v_w = args.array<int>("csrColumnOffsets_v_w");
635 xt::pyarray<int>& csrRowIndeces_w_u = args.array<int>("csrRowIndeces_w_u");
636 xt::pyarray<int>& csrColumnOffsets_w_u = args.array<int>("csrColumnOffsets_w_u");
637 xt::pyarray<int>& csrRowIndeces_w_v = args.array<int>("csrRowIndeces_w_v");
638 xt::pyarray<int>& csrColumnOffsets_w_v = args.array<int>("csrColumnOffsets_w_v");
639 xt::pyarray<int>& csrRowIndeces_w_w = args.array<int>("csrRowIndeces_w_w");
640 xt::pyarray<int>& csrColumnOffsets_w_w = args.array<int>("csrColumnOffsets_w_w");
641 xt::pyarray<double>& globalJacobian = args.array<double>("globalJacobian");
642 int nExteriorElementBoundaries_global = args.scalar<int>("nExteriorElementBoundaries_global");
643 xt::pyarray<int>& exteriorElementBoundariesArray = args.array<int>("exteriorElementBoundariesArray");
644 xt::pyarray<int>& elementBoundaryElementsArray = args.array<int>("elementBoundaryElementsArray");
645 xt::pyarray<int>& elementBoundaryLocalElementBoundariesArray = args.array<int>("elementBoundaryLocalElementBoundariesArray");
646 xt::pyarray<int>& isDOFBoundary_u = args.array<int>("isDOFBoundary_u");
647 xt::pyarray<int>& isDOFBoundary_v = args.array<int>("isDOFBoundary_v");
648 xt::pyarray<int>& isDOFBoundary_w = args.array<int>("isDOFBoundary_w");
649 xt::pyarray<int>& isStressFluxBoundary_u = args.array<int>("isStressFluxBoundary_u");
650 xt::pyarray<int>& isStressFluxBoundary_v = args.array<int>("isStressFluxBoundary_v");
651 xt::pyarray<int>& isStressFluxBoundary_w = args.array<int>("isStressFluxBoundary_w");
652 xt::pyarray<int>& csrColumnOffsets_eb_u_u = args.array<int>("csrColumnOffsets_eb_u_u");
653 xt::pyarray<int>& csrColumnOffsets_eb_u_v = args.array<int>("csrColumnOffsets_eb_u_v");
654 xt::pyarray<int>& csrColumnOffsets_eb_u_w = args.array<int>("csrColumnOffsets_eb_u_w");
655 xt::pyarray<int>& csrColumnOffsets_eb_v_u = args.array<int>("csrColumnOffsets_eb_v_u");
656 xt::pyarray<int>& csrColumnOffsets_eb_v_v = args.array<int>("csrColumnOffsets_eb_v_v");
657 xt::pyarray<int>& csrColumnOffsets_eb_v_w = args.array<int>("csrColumnOffsets_eb_v_w");
658 xt::pyarray<int>& csrColumnOffsets_eb_w_u = args.array<int>("csrColumnOffsets_eb_w_u");
659 xt::pyarray<int>& csrColumnOffsets_eb_w_v = args.array<int>("csrColumnOffsets_eb_w_v");
660 xt::pyarray<int>& csrColumnOffsets_eb_w_w = args.array<int>("csrColumnOffsets_eb_w_w");
661 CompKernel<nSpace,nDOF_mesh_trial_element,nDOF_trial_element,nDOF_test_element> ck;
662 const int nSymTen(ck.nSymTen);
663 //
664 //loop over elements to compute volume integrals and load them into the element Jacobians and global Jacobian
665 //
666 for(int eN=0;eN<nElements_global;eN++)
667 {
668 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[2],
841 x_ext,y_ext,z_ext,
842 G[nSpace*nSpace],G_dd_G,tr_G,h_penalty,
843 strain[nSymTen],
844 stress[nSymTen],
845 dstress[nSymTen*nSymTen];
846 ck.calculateMapping_elementBoundary(eN,
847 ebN_local,
848 kb,
849 ebN_local_kb,
850 mesh_dof.data(),
851 mesh_l2g.data(),
852 mesh_trial_trace_ref.data(),
853 mesh_grad_trial_trace_ref.data(),
854 boundaryJac_ref.data(),
855 jac_ext,
856 jacDet_ext,
857 jacInv_ext,
858 boundaryJac,
859 metricTensor,
860 metricTensorDetSqrt,
861 normal_ref.data(),
862 normal,
863 x_ext,y_ext,z_ext);
864 dS = metricTensorDetSqrt*dS_ref.data()[kb];
865 ck.calculateG(jacInv_ext,G,G_dd_G,tr_G);
866 //compute shape and solution information
867 //shape
868 ck.gradTrialFromRef(&disp_grad_trial_trace_ref.data()[ebN_local_kb_nSpace*nDOF_trial_element],jacInv_ext,disp_grad_trial_trace);
869 //solution and gradients
870 ck.valFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],u_ext);
871 ck.valFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],v_ext);
872 /* ck.valFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],&disp_trial_trace_ref.data()[ebN_local_kb*nDOF_test_element],w_ext); */
873 ck.gradFromDOF(u_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_u_ext);
874 ck.gradFromDOF(v_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_v_ext);
875 /* ck.gradFromDOF(w_dof.data(),&disp_l2g.data()[eN_nDOF_trial_element],disp_grad_trial_trace,grad_w_ext); */
876 //precalculate test function products with integration weights
877 for (int j=0;j<nDOF_trial_element;j++)
878 {
879 disp_test_dS[j] = disp_test_trace_ref.data()[ebN_local_kb*nDOF_test_element+j]*dS;
880 }
881 //
882 //calculate the internal and external trace of the pde coefficients
883 //
884 calculateStrain(D_ext,strain);
885 evaluateCoefficients(fabs(jacDet_ext),
886 &materialProperties.data()[materialTypes.data()[eN]*nMaterialProperties],
887 strain,
888 stress,
889 dstress);
890 //
891 //calculate the flux jacobian
892 //
893 ck.calculateGScale(G,normal,h_penalty);
894 //cek hack
895 h_penalty = 100.0/h_penalty;
896 // for (int II=0;II<nSymTen;II++)
897 // for (int JJ=0;<nSymTen;JJ++)
898 // dstress[II*nSymTen+JJ] = 0.0;
899 for (int j=0;j<nDOF_trial_element;j++)
900 {
901 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 };//MoveMesh2D
948
949 inline MoveMesh2D_base* newMoveMesh2D(int nSpaceIn,
950 int nQuadraturePoints_elementIn,
951 int nDOF_mesh_trial_elementIn,
952 int nDOF_trial_elementIn,
953 int nDOF_test_elementIn,
954 int nQuadraturePoints_elementBoundaryIn,
955 int CompKernelFlag)
956 {
957 return proteus::chooseAndAllocateDiscretization2D<MoveMesh2D_base,MoveMesh2D,CompKernel>(nSpaceIn,
958 nQuadraturePoints_elementIn,
959 nDOF_mesh_trial_elementIn,
960 nDOF_trial_elementIn,
961 nDOF_test_elementIn,
962 nQuadraturePoints_elementBoundaryIn,
963 CompKernelFlag);
964 }
965}//proteus
966
967#endif
u
Double u
Definition Headers.h:89
z
Double * z
Definition Headers.h:49
v
Double v
Definition Headers.h:95
proteus::MoveMesh2D_base
Definition MoveMesh2D.h:15
proteus::MoveMesh2D_base::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)=0
proteus::MoveMesh2D_base::calculateResidual
virtual void calculateResidual(arguments_dict &args)=0
proteus::MoveMesh2D_base::~MoveMesh2D_base
virtual ~MoveMesh2D_base()
Definition MoveMesh2D.h:17
proteus::MoveMesh2D::ei
EIndex< nSpace > ei
Definition MoveMesh2D.h:36
proteus::MoveMesh2D::exteriorNumericalStressFluxJacobian
void exteriorNumericalStressFluxJacobian(const int &isDOFBoundary_u, const int &isDOFBoundary_v, const double *normal, const double *dstress, const double &penalty, const double &disp_trial, const double *disp_grad_trial, double &dstressFlux_u_u, double &dstressFlux_u_v, double &dstressFlux_v_u, double &dstressFlux_v_v)
Definition MoveMesh2D.h:187
proteus::MoveMesh2D::calculateStrain
void calculateStrain(double *D, double *strain)
Definition MoveMesh2D.h:72
proteus::MoveMesh2D::calculateJacobian
virtual void calculateJacobian(arguments_dict &args)
Definition MoveMesh2D.h:594
proteus::MoveMesh2D::calculateResidual
virtual void calculateResidual(arguments_dict &args)
Definition MoveMesh2D.h:280
proteus::MoveMesh2D::exteriorNumericalStressFlux
void exteriorNumericalStressFlux(const int &isDOFBoundary_u, const int &isDOFBoundary_v, const int &isStressFluxBoundary_u, const int &isStressFluxBoundary_v, const double &penalty, const double &u, const double &v, const double &bc_u, const double &bc_v, const double &bc_stressFlux_u, const double &bc_stressFlux_v, const double *stress, const double *normal, double &stressFlux_u, double &stressFlux_v)
Definition MoveMesh2D.h:125
proteus::MoveMesh2D::ck
CompKernelType ck
Definition MoveMesh2D.h:32
proteus::MoveMesh2D::evaluateCoefficients
void evaluateCoefficients(double det_J, const double *materialProperties, double *strain, double *stress, double *dstress)
Definition MoveMesh2D.h:83
proteus::MoveMesh2D::nDOF_test_X_trial_element
const int nDOF_test_X_trial_element
Definition MoveMesh2D.h:34
D
#define D(x, y)
Definition jf.h:9
proteus
Definition ADR.h:19
proteus::newMoveMesh2D
MoveMesh2D_base * newMoveMesh2D(int nSpaceIn, int nQuadraturePoints_elementIn, int nDOF_mesh_trial_elementIn, int nDOF_trial_elementIn, int nDOF_test_elementIn, int nQuadraturePoints_elementBoundaryIn, int CompKernelFlag)
Definition MoveMesh2D.h:949
proteus::chooseAndAllocateDiscretization2D
Model_Base * chooseAndAllocateDiscretization2D(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:309
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