proteus 1.9.0
C/C++/Fortran libraries
Loading...
Searching...
No Matches
CompKernel.h
Go to the documentation of this file.
1#ifndef COMPKERNEL_H
2#define COMPKERNEL_H
3#include <cmath>
4//#include "xtensor-python/pyarray.hpp"
5//#include "xtensor-python/pyvectorize.hpp"
9template<int NSPACE>
10class EIndex
11{
12};
13
14template<>
15class EIndex<3>
16{
17public:
18 const int X,Y,Z,
32 X(0),
33 Y(1),
34 Z(2),
35 XX(0),XY(1),XZ(2),
36 YX(3),YY(4),YZ(5),
37 ZX(6),ZY(7),ZZ(8),
38 sXX(0),sXY(5),sXZ(4),
39 sYX(5),sYY(1),sYZ(3),
40 sZX(4),sZY(3),sZZ(2),
41 nSymTen(6),
42 XHX(0),XHY(1),
43 YHX(2),YHY(3),
44 ZHX(4),ZHY(5),
45 HXHX(0),HXHY(1),
46 HYHX(2),HYHY(3)
47 {}
48};
49
50template<>
51class EIndex<2>
52{
53public:
54 const int X,Y,
64 X(0),
65 Y(1),
66 XX(0),XY(1),
67 YX(2),YY(3),
68 sXX(0),sXY(2),
69 sYX(2),sYY(1),
70 nSymTen(3),
71 XHX(0),XHY(1),
72 YHX(2),YHY(3),
73 HXHX(0)
74 {}
75};
76
77template<>
78class EIndex<1>
79{
80public:
81 const int X,
88 X(0),
89 XX(0),
90 sXX(0),
91 nSymTen(1),
92 XHX(0),
93 HXHX(0)
94 {}
95};
96
97//I separated the space mapping part of the kernel so I could partially specialize the template on NSPACE
98template<int NSPACE, int NDOF_MESH_TRIAL_ELEMENT>
100{
101public:
102 inline void calculateMapping_element(const int eN,
103 const int k,
104 double* mesh_dof,
105 int* mesh_l2g,
106 //xt::pyarray<double>& mesh_trial_ref,
107 double* mesh_trial_ref,
108 double* mesh_grad_trial_ref,
109 double* jac,
110 double& jacDet,
111 double* jacInv,
112 double& x,
113 double& y,
114 double& z);
115 inline void calculateH_element(const int eN,
116 const int k,
117 double* h_dof,
118 int* mesh_l2g,
119 //xt::pyarray<double>& mesh_trial_ref,
120 double* mesh_trial_ref,
121 double& h);
122 inline void calculateMappingVelocity_element(const int eN,
123 const int k,
124 double* mesh_velocity_dof,
125 int* mesh_l2g,
126 //xt::pyarray<double>& mesh_trial_ref,
127 double* mesh_trial_ref,
128 double& xt,
129 double& yt,
130 double& zt);
131 inline void calculateMapping_elementBoundary(const int eN,
132 const int ebN_local,
133 const int kb,
134 const int ebN_local_kb,
135 double* mesh_dof,
136 int* mesh_l2g,
137 double* mesh_trial_trace_ref,
138 double* mesh_grad_trial_trace_ref,
139 double* boundaryJac_ref,
140 double* jac,
141 double& jacDet,
142 double* jacInv,
143 double* boundaryJac,
144 double* metricTensor,
145 double& metricTensorDetSqrt,
146 double* normal_ref,
147 double* normal,
148 double& x,
149 double& y,
150 double& z);
152 const int ebN_local,
153 const int kb,
154 const int ebN_local_kb,
155 double* mesh_velocity_dof,
156 int* mesh_l2g,
157 double* mesh_trial_trace_ref,
158 double& xt,
159 double& yt,
160 double& zt,
161 double* normal,
162 double* boundaryJac,
163 double* metricTensor,
164 double& metricTensorDetSqrt);
165 inline void valFromDOF(const double* dof,const int* l2g_element,const double* trial_ref,double& val)
166 {
167 val=0.0;
168 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
169 val+=dof[l2g_element[j]]*trial_ref[j];
170 }
171
172 inline void gradFromDOF(const double* dof,const int* l2g_element,const double* grad_trial,double* grad)
173 {
174 for(int I=0;I<NSPACE;I++)
175 grad[I] = 0.0;
176 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
177 for(int I=0;I<NSPACE;I++)
178 grad[I] += dof[l2g_element[j]]*grad_trial[j*NSPACE+I];
179 }
180
181 inline void hessFromDOF(const double* dof,const int* l2g_element,const double* hess_trial,double* hess)
182 {
183 const int NSPACE2=NSPACE*NSPACE;
184 for(int I=0;I<NSPACE;I++)
185 for(int J=0;J<NSPACE;J++)
186 hess[I*NSPACE+J] = 0.0;
187 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
188 for(int I=0;I<NSPACE;I++)
189 for(int J=0;J<NSPACE;J++)
190 hess[I*NSPACE+J] += dof[l2g_element[j]]*hess_trial[j*NSPACE2+I*NSPACE+J];
191 }
192};
193
194//specialization for 3D
195template<int NDOF_MESH_TRIAL_ELEMENT>
196class CompKernelSpaceMapping<3,NDOF_MESH_TRIAL_ELEMENT>
197{
198public:
199 const int X,Y,Z,
212 ;
214 X(0),
215 Y(1),
216 Z(2),
217 XX(0),XY(1),XZ(2),
218 YX(3),YY(4),YZ(5),
219 ZX(6),ZY(7),ZZ(8),
220 sXX(0),sXY(5),sXZ(4),
221 sYX(5),sYY(1),sYZ(3),
222 sZX(4),sZY(3),sZZ(2),
223 nSymTen(6),
224 XHX(0),XHY(1),
225 YHX(2),YHY(3),
226 ZHX(4),ZHY(5),
227 HXHX(0),HXHY(1),
228 HYHX(2),HYHY(3)
229 {}
230 inline void calculateMapping_element(const int eN,
231 const int k,
232 double* mesh_dof,
233 int* mesh_l2g,
234 //xt::pyarray<double>& mesh_trial_ref,
235 double* mesh_trial_ref,
236 double* mesh_grad_trial_ref,
237 double* jac,
238 double& jacDet,
239 double* jacInv,
240 double& x,
241 double& y,
242 double& z)
243 {
244 double Grad_x[3],Grad_y[3],Grad_z[3],oneOverJacDet;
245
246 //
247 //mapping of reference element to physical element
248 //
249 x=0.0;y=0.0;z=0.0;
250 for (int I=0;I<3;I++)
251 {
252 Grad_x[I]=0.0;Grad_y[I]=0.0;Grad_z[I]=0.0;
253 }
254 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
255 {
256 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
257 //x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref(k, j);
258 //y += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref(k, j);
259 //z += mesh_dof[mesh_l2g[eN_j]*3+2]*mesh_trial_ref(k, j);
260 x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
261 y += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
262 z += mesh_dof[mesh_l2g[eN_j]*3+2]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
263 for (int I=0;I<3;I++)
264 {
265 Grad_x[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT*3+j*3+I];
266 Grad_y[I] += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT*3+j*3+I];
267 Grad_z[I] += mesh_dof[mesh_l2g[eN_j]*3+2]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT*3+j*3+I];
268 }
269 }
270 jac[XX] = Grad_x[X];
271 jac[XY] = Grad_x[Y];
272 jac[XZ] = Grad_x[Z];
273 jac[YX] = Grad_y[X];
274 jac[YY] = Grad_y[Y];
275 jac[YZ] = Grad_y[Z];
276 jac[ZX] = Grad_z[X];
277 jac[ZY] = Grad_z[Y];
278 jac[ZZ] = Grad_z[Z];
279 jacDet =
280 jac[XX]*(jac[YY]*jac[ZZ] - jac[YZ]*jac[ZY]) -
281 jac[XY]*(jac[YX]*jac[ZZ] - jac[YZ]*jac[ZX]) +
282 jac[XZ]*(jac[YX]*jac[ZY] - jac[YY]*jac[ZX]);
283 oneOverJacDet = 1.0/jacDet;
284 jacInv[XX] = oneOverJacDet*(jac[YY]*jac[ZZ] - jac[YZ]*jac[ZY]);
285 jacInv[YX] = oneOverJacDet*(jac[YZ]*jac[ZX] - jac[YX]*jac[ZZ]);
286 jacInv[ZX] = oneOverJacDet*(jac[YX]*jac[ZY] - jac[YY]*jac[ZX]);
287 jacInv[XY] = oneOverJacDet*(jac[ZY]*jac[XZ] - jac[ZZ]*jac[XY]);
288 jacInv[YY] = oneOverJacDet*(jac[ZZ]*jac[XX] - jac[ZX]*jac[XZ]);
289 jacInv[ZY] = oneOverJacDet*(jac[ZX]*jac[XY] - jac[ZY]*jac[XX]);
290 jacInv[XZ] = oneOverJacDet*(jac[XY]*jac[YZ] - jac[XZ]*jac[YY]);
291 jacInv[YZ] = oneOverJacDet*(jac[XZ]*jac[YX] - jac[XX]*jac[YZ]);
292 jacInv[ZZ] = oneOverJacDet*(jac[XX]*jac[YY] - jac[XY]*jac[YX]);
293 }
294
295 inline void calculateH_element(const int eN,
296 const int k,
297 double* h_dof,
298 int* mesh_l2g,
299 //xt::pyarray<double>& mesh_trial_ref,
300 double* mesh_trial_ref,
301 double& h)
302 {
303 h=0.0;
304 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
305 {
306 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
307 //h += h_dof[mesh_l2g[eN_j]]*mesh_trial_ref(k, j);
308 h += h_dof[mesh_l2g[eN_j]]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
309 }
310 }
311
312 inline void calculateMappingVelocity_element(const int eN,
313 const int k,
314 double* mesh_velocity_dof,
315 int* mesh_l2g,
316 //xt::pyarray<double>& mesh_trial_ref,
317 double* mesh_trial_ref,
318 double& xt,
319 double& yt,
320 double& zt)
321 {
322 //
323 //time derivative of mapping of reference element to physical element
324 //
325 xt=0.0;yt=0.0;zt=0.0;
326 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
327 {
328 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
329 //xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref(k, j);
330 //yt += mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref(k, j);
331 //zt += mesh_velocity_dof[mesh_l2g[eN_j]*3+2]*mesh_trial_ref(k, j);
332 xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
333 yt += mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
334 zt += mesh_velocity_dof[mesh_l2g[eN_j]*3+2]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
335 }
336 }
337
338 inline void calculateMapping_elementBoundary(const int eN,
339 const int ebN_local,
340 const int kb,
341 const int ebN_local_kb,
342 double* mesh_dof,
343 int* mesh_l2g,
344 double* mesh_trial_trace_ref,
345 double* mesh_grad_trial_trace_ref,
346 double* boundaryJac_ref,
347 double* jac,
348 double& jacDet,
349 double* jacInv,
350 double* boundaryJac,
351 double* metricTensor,
352 double& metricTensorDetSqrt,
353 double* normal_ref,
354 double* normal,
355 double& x,
356 double& y,
357 double& z)
358 {
359 const int ebN_local_kb_nSpace = ebN_local_kb*3,
360 ebN_local_kb_nSpace_nSpacem1 = ebN_local_kb*3*2;
361
362 double Grad_x_ext[3],Grad_y_ext[3],Grad_z_ext[3],oneOverJacDet,norm_normal=0.0;
363 //
364 //calculate mapping from the reference element to the physical element
365 //
366 x=0.0;y=0.0;z=0.0;
367 for (int I=0;I<3;I++)
368 {
369 Grad_x_ext[I] = 0.0;
370 Grad_y_ext[I] = 0.0;
371 Grad_z_ext[I] = 0.0;
372 }
373 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
374 {
375 int eN_j = eN*NDOF_MESH_TRIAL_ELEMENT+j;
376 int ebN_local_kb_j = ebN_local_kb*NDOF_MESH_TRIAL_ELEMENT+j;
377 int ebN_local_kb_j_nSpace = ebN_local_kb_j*3;
378 x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_trace_ref[ebN_local_kb_j];
379 y += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_trace_ref[ebN_local_kb_j];
380 z += mesh_dof[mesh_l2g[eN_j]*3+2]*mesh_trial_trace_ref[ebN_local_kb_j];
381 for (int I=0;I<3;I++)
382 {
383 Grad_x_ext[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I];
384 Grad_y_ext[I] += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I];
385 Grad_z_ext[I] += mesh_dof[mesh_l2g[eN_j]*3+2]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I];
386 }
387 }
388 //Space Mapping Jacobian
389 jac[XX] = Grad_x_ext[X];
390 jac[XY] = Grad_x_ext[Y];
391 jac[XZ] = Grad_x_ext[Z];
392 jac[YX] = Grad_y_ext[X];
393 jac[YY] = Grad_y_ext[Y];
394 jac[YZ] = Grad_y_ext[Z];
395 jac[ZX] = Grad_z_ext[X];
396 jac[ZY] = Grad_z_ext[Y];
397 jac[ZZ] = Grad_z_ext[Z];
398 jacDet =
399 jac[XX]*(jac[YY]*jac[ZZ] - jac[YZ]*jac[ZY]) -
400 jac[XY]*(jac[YX]*jac[ZZ] - jac[YZ]*jac[ZX]) +
401 jac[XZ]*(jac[YX]*jac[ZY] - jac[YY]*jac[ZX]);
402 oneOverJacDet = 1.0/jacDet;
403 jacInv[XX] = oneOverJacDet*(jac[YY]*jac[ZZ] - jac[YZ]*jac[ZY]);
404 jacInv[YX] = oneOverJacDet*(jac[YZ]*jac[ZX] - jac[YX]*jac[ZZ]);
405 jacInv[ZX] = oneOverJacDet*(jac[YX]*jac[ZY] - jac[YY]*jac[ZX]);
406 jacInv[XY] = oneOverJacDet*(jac[ZY]*jac[XZ] - jac[ZZ]*jac[XY]);
407 jacInv[YY] = oneOverJacDet*(jac[ZZ]*jac[XX] - jac[ZX]*jac[XZ]);
408 jacInv[ZY] = oneOverJacDet*(jac[ZX]*jac[XY] - jac[ZY]*jac[XX]);
409 jacInv[XZ] = oneOverJacDet*(jac[XY]*jac[YZ] - jac[XZ]*jac[YY]);
410 jacInv[YZ] = oneOverJacDet*(jac[XZ]*jac[YX] - jac[XX]*jac[YZ]);
411 jacInv[ZZ] = oneOverJacDet*(jac[XX]*jac[YY] - jac[XY]*jac[YX]);
412 //normal
413 norm_normal=0.0;
414 for (int I=0;I<3;I++)
415 normal[I] = 0.0;
416 for (int I=0;I<3;I++)
417 {
418 for (int J=0;J<3;J++)
419 {
420 normal[I] += jacInv[J*3+I]*normal_ref[ebN_local_kb_nSpace+J];
421 }
422 norm_normal+=normal[I]*normal[I];
423 }
424 norm_normal = sqrt(norm_normal);
425 for (int I=0;I<3;I++)
426 {
427 normal[I] /= norm_normal;
428 }
429 //metric tensor and determinant
430 boundaryJac[XHX] = jac[XX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHX]+jac[XY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHX]+jac[XZ]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+ZHX];
431 boundaryJac[XHY] = jac[XX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHY]+jac[XY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHY]+jac[XZ]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+ZHY];
432 boundaryJac[YHX] = jac[YX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHX]+jac[YY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHX]+jac[YZ]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+ZHX];
433 boundaryJac[YHY] = jac[YX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHY]+jac[YY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHY]+jac[YZ]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+ZHY];
434 boundaryJac[ZHX] = jac[ZX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHX]+jac[ZY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHX]+jac[ZZ]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+ZHX];
435 boundaryJac[ZHY] = jac[ZX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHY]+jac[ZY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHY]+jac[ZZ]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+ZHY];
436
437 metricTensor[HXHX] = boundaryJac[XHX]*boundaryJac[XHX]+boundaryJac[YHX]*boundaryJac[YHX]+boundaryJac[ZHX]*boundaryJac[ZHX];
438 metricTensor[HXHY] = boundaryJac[XHX]*boundaryJac[XHY]+boundaryJac[YHX]*boundaryJac[YHY]+boundaryJac[ZHX]*boundaryJac[ZHY];
439 metricTensor[HYHX] = boundaryJac[XHY]*boundaryJac[XHX]+boundaryJac[YHY]*boundaryJac[YHX]+boundaryJac[ZHY]*boundaryJac[ZHX];
440 metricTensor[HYHY] = boundaryJac[XHY]*boundaryJac[XHY]+boundaryJac[YHY]*boundaryJac[YHY]+boundaryJac[ZHY]*boundaryJac[ZHY];
441
442 metricTensorDetSqrt=sqrt(metricTensor[HXHX]*metricTensor[HYHY]- metricTensor[HXHY]*metricTensor[HYHX]);
443 }
444
446 const int ebN_local,
447 const int kb,
448 const int ebN_local_kb,
449 double* mesh_velocity_dof,
450 int* mesh_l2g,
451 double* mesh_trial_trace_ref,
452 double& xt,
453 double& yt,
454 double& zt,
455 double* normal,
456 double* boundaryJac,
457 double* metricTensor,
458 double& metricTensorDetSqrt)
459 {
460 //const int ebN_local_kb_nSpace = ebN_local_kb*3,
461 // ebN_local_kb_nSpace_nSpacem1 = ebN_local_kb*3*2;
462 //
463 //calculate velocity of mapping from the reference element to the physical element
464 //
465 xt=0.0;yt=0.0;zt=0.0;
466 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
467 {
468 int eN_j = eN*NDOF_MESH_TRIAL_ELEMENT+j;
469 int ebN_local_kb_j = ebN_local_kb*NDOF_MESH_TRIAL_ELEMENT+j;
470 xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_trace_ref[ebN_local_kb_j];
471 yt += mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_trace_ref[ebN_local_kb_j];
472 zt += mesh_velocity_dof[mesh_l2g[eN_j]*3+2]*mesh_trial_trace_ref[ebN_local_kb_j];
473 }
474 //modify the metricTensorDetSqrt to include the effect of the moving domain
475 //it's not exactly the Sqrt(Det(G_tr_G)) now, see notes
476 //just do it brute force
477 double
478 Gy_tr_Gy_00 = 1.0 + xt*xt + yt*yt + zt*zt,
479 Gy_tr_Gy_01 = boundaryJac[XHX]*xt+boundaryJac[YHX]*yt+boundaryJac[ZHX]*zt,
480 Gy_tr_Gy_02 = boundaryJac[XHY]*xt+boundaryJac[YHY]*yt+boundaryJac[ZHY]*zt,
481 Gy_tr_Gy_10 = Gy_tr_Gy_01,
482 Gy_tr_Gy_20 = Gy_tr_Gy_02,
483 Gy_tr_Gy_11 = metricTensor[HXHX],
484 Gy_tr_Gy_12 = metricTensor[HXHY],
485 Gy_tr_Gy_21 = metricTensor[HYHX],
486 Gy_tr_Gy_22 = metricTensor[HYHY],
487 xt_dot_n = xt*normal[X]+yt*normal[Y]+zt*normal[Z];
488 metricTensorDetSqrt=sqrt((Gy_tr_Gy_00*Gy_tr_Gy_11*Gy_tr_Gy_22 +
489 Gy_tr_Gy_01*Gy_tr_Gy_12*Gy_tr_Gy_20 +
490 Gy_tr_Gy_02*Gy_tr_Gy_10*Gy_tr_Gy_21 -
491 Gy_tr_Gy_20*Gy_tr_Gy_11*Gy_tr_Gy_02 -
492 Gy_tr_Gy_21*Gy_tr_Gy_12*Gy_tr_Gy_00 -
493 Gy_tr_Gy_22*Gy_tr_Gy_10*Gy_tr_Gy_01) / (1.0+xt_dot_n*xt_dot_n));
494 }
495 inline void valFromDOF(const double* dof,const int* l2g_element,const double* trial_ref,double& val)
496 {
497 val=0.0;
498 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
499 val+=dof[l2g_element[j]]*trial_ref[j];
500 }
501
502 inline void gradFromDOF(const double* dof,const int* l2g_element,const double* grad_trial,double* grad)
503 {
504 for(int I=0;I<3;I++)
505 grad[I] = 0.0;
506 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
507 for(int I=0;I<3;I++)
508 grad[I] += dof[l2g_element[j]]*grad_trial[j*3+I];
509 }
510
511 inline void hessFromDOF(const double* dof,const int* l2g_element,const double* hess_trial,double* hess)
512 {
513 for(int I=0;I<3;I++)
514 for(int J=0;J<3;J++)
515 hess[I*3+J] = 0.0;
516 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
517 for(int I=0;I<3;I++)
518 for(int J=0;J<3;J++)
519 hess[I*3+J] += dof[l2g_element[j]]*hess_trial[j*9+I*3+J];
520 }
521
522 inline void calculateMapping_element(const int eN,
523 const int k,
524 double* mesh_dof,
525 int* mesh_l2g,
526 //xt::pyarray<double>& mesh_trial_ref,
527 double* mesh_trial_ref,
528 double* mesh_grad_trial_ref,
529 double* jac,
530 double& jacDet,
531 double* jacInv,
532 double& x,
533 double& y)
534 {
536 k,
537 mesh_dof,
538 mesh_l2g,
539 mesh_trial_ref,
540 mesh_grad_trial_ref,
541 jac,
542 jacDet,
543 jacInv,
544 x,
545 y,
546 0.0);
547 }
548
549 inline void calculateMappingVelocity_element(const int eN,
550 const int k,
551 double* mesh_velocity_dof,
552 int* mesh_l2g,
553 //xt::pyarray<double>& mesh_trial_ref,
554 double* mesh_trial_ref,
555 double& xt,
556 double& yt)
557 {
559 k,
560 mesh_velocity_dof,
561 mesh_l2g,
562 mesh_trial_ref,
563 xt,
564 yt,
565 0.0);
566 }
567 inline void calculateMapping_elementBoundary(const int eN,
568 const int ebN_local,
569 const int kb,
570 const int ebN_local_kb,
571 double* mesh_dof,
572 int* mesh_l2g,
573 double* mesh_trial_trace_ref,
574 double* mesh_grad_trial_trace_ref,
575 double* boundaryJac_ref,
576 double* jac,
577 double& jacDet,
578 double* jacInv,
579 double* boundaryJac,
580 double* metricTensor,
581 double& metricTensorDetSqrt,
582 double* normal_ref,
583 double* normal,
584 double& x,
585 double& y)
586 {
588 ebN_local,
589 kb,
590 ebN_local_kb,
591 mesh_dof,
592 mesh_l2g,
593 mesh_trial_trace_ref,
594 mesh_grad_trial_trace_ref,
595 boundaryJac_ref,
596 jac,
597 jacDet,
598 jacInv,
599 boundaryJac,
600 metricTensor,
601 metricTensorDetSqrt,
602 normal_ref,
603 normal,
604 x,
605 y,
606 0.0);
607 }
609 const int ebN_local,
610 const int kb,
611 const int ebN_local_kb,
612 double* mesh_velocity_dof,
613 int* mesh_l2g,
614 double* mesh_trial_trace_ref,
615 double& xt,
616 double& yt,
617 double* normal,
618 double* boundaryJac,
619 double* metricTensor,
620 double& metricTensorDetSqrt)
621 {
623 ebN_local,
624 kb,
625 ebN_local_kb,
626 mesh_velocity_dof,
627 mesh_l2g,
628 mesh_trial_trace_ref,
629 xt,
630 yt,
631 0.0,
632 normal,
633 boundaryJac,
634 metricTensor,
635 metricTensorDetSqrt);
636 }
637
638};
639
640//specialization for 2D
641template<int NDOF_MESH_TRIAL_ELEMENT>
642class CompKernelSpaceMapping<2,NDOF_MESH_TRIAL_ELEMENT>
643{
644public:
645 const int X,Y,
655 X(0),
656 Y(1),
657 XX(0),XY(1),
658 YX(2),YY(3),
659 sXX(0),sXY(2),
660 sYX(2),sYY(1),
661 nSymTen(3),
662 XHX(0),
663 YHX(1),
664 HXHX(0)
665 {}
666 inline void calculateMapping_element(const int eN,
667 const int k,
668 double* mesh_dof,
669 int* mesh_l2g,
670 //xt::pyarray<double>& mesh_trial_ref,
671 double* mesh_trial_ref,
672 double* mesh_grad_trial_ref,
673 double* jac,
674 double& jacDet,
675 double* jacInv,
676 double& x,
677 double& y)
678 {
679 double Grad_x[2],Grad_y[2],oneOverJacDet;
680
681 //
682 //mapping of reference element to physical element
683 //
684 x=0.0;y=0.0;
685 for (int I=0;I<2;I++)
686 {
687 Grad_x[I]=0.0;Grad_y[I]=0.0;
688 }
689 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
690 {
691 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
692 /* x += mesh_dof[mesh_l2g[eN_j]*2+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j]; */
693 /* y += mesh_dof[mesh_l2g[eN_j]*2+1]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j]; */
694 /* for (int I=0;I<2;I++) */
695 /* { */
696 /* Grad_x[I] += mesh_dof[mesh_l2g[eN_j]*2+0]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT*2+j*2+I]; */
697 /* Grad_y[I] += mesh_dof[mesh_l2g[eN_j]*2+1]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT*2+j*2+I]; */
698 /* } */
699 //x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref(k, j);
700 //y += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref(k, j);
701 x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
702 y += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
703 for (int I=0;I<2;I++)
704 {
705 Grad_x[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT*2+j*2+I];
706 Grad_y[I] += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT*2+j*2+I];
707 }
708 }
709 jac[XX] = Grad_x[X];
710 jac[XY] = Grad_x[Y];
711 jac[YX] = Grad_y[X];
712 jac[YY] = Grad_y[Y];
713 jacDet = jac[XX]*jac[YY] - jac[XY]*jac[YX];
714 oneOverJacDet = 1.0/jacDet;
715 jacInv[XX] = oneOverJacDet*jac[YY];
716 jacInv[XY] = -oneOverJacDet*jac[XY];
717 jacInv[YX] = -oneOverJacDet*jac[YX];
718 jacInv[YY] = oneOverJacDet*jac[XX];
719 }
720
721 inline void calculateH_element(const int eN,
722 const int k,
723 double* h_dof,
724 int* mesh_l2g,
725 //xt::pyarray<double>& mesh_trial_ref,
726 double* mesh_trial_ref,
727 double& h)
728 {
729 h=0.0;
730 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
731 {
732 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
733 //h += h_dof[mesh_l2g[eN_j]]*mesh_trial_ref(k, j);
734 h += h_dof[mesh_l2g[eN_j]]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
735 }
736 }
737
738 inline void calculateMappingVelocity_element(const int eN,
739 const int k,
740 double* mesh_velocity_dof,
741 int* mesh_l2g,
742 //xt::pyarray<double>& mesh_trial_ref,
743 double* mesh_trial_ref,
744 double& xt,
745 double& yt)
746 {
747 //
748 //time derivative of mapping of reference element to physical element
749 //
750 xt=0.0;yt=0.0;
751 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
752 {
753 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
754 /* xt += mesh_velocity_dof[mesh_l2g[eN_j]*2+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j]; */
755 /* yt += mesh_velocity_dof[mesh_l2g[eN_j]*2+1]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j]; */
756 //xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref(k, j);
757 //yt += mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref(k, j);
758 xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
759 yt += mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
760 }
761 }
762
763 inline void calculateMapping_elementBoundary(const int eN,
764 const int ebN_local,
765 const int kb,
766 const int ebN_local_kb,
767 double* mesh_dof,
768 int* mesh_l2g,
769 double* mesh_trial_trace_ref,
770 double* mesh_grad_trial_trace_ref,
771 double* boundaryJac_ref,
772 double* jac,
773 double& jacDet,
774 double* jacInv,
775 double* boundaryJac,
776 double* metricTensor,
777 double& metricTensorDetSqrt,
778 double* normal_ref,
779 double* normal,
780 double& x,
781 double& y)
782 {
783 const int ebN_local_kb_nSpace = ebN_local_kb*2,
784 ebN_local_kb_nSpace_nSpacem1 = ebN_local_kb*2*1;
785
786 double Grad_x_ext[2],Grad_y_ext[2],oneOverJacDet,norm_normal=0.0;
787 //
788 //calculate mapping from the reference element to the physical element
789 //
790 x=0.0;y=0.0;
791 for (int I=0;I<2;I++)
792 {
793 Grad_x_ext[I] = 0.0;
794 Grad_y_ext[I] = 0.0;
795 }
796 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
797 {
798 int eN_j = eN*NDOF_MESH_TRIAL_ELEMENT+j;
799 int ebN_local_kb_j = ebN_local_kb*NDOF_MESH_TRIAL_ELEMENT+j;
800 int ebN_local_kb_j_nSpace = ebN_local_kb_j*2;
801 /* x += mesh_dof[mesh_l2g[eN_j]*2+0]*mesh_trial_trace_ref[ebN_local_kb_j]; */
802 /* y += mesh_dof[mesh_l2g[eN_j]*2+1]*mesh_trial_trace_ref[ebN_local_kb_j]; */
803 /* for (int I=0;I<2;I++) */
804 /* { */
805 /* Grad_x_ext[I] += mesh_dof[mesh_l2g[eN_j]*2+0]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I]; */
806 /* Grad_y_ext[I] += mesh_dof[mesh_l2g[eN_j]*2+1]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I]; */
807 /* } */
808 x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_trace_ref[ebN_local_kb_j];
809 y += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_trace_ref[ebN_local_kb_j];
810 for (int I=0;I<2;I++)
811 {
812 Grad_x_ext[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I];
813 Grad_y_ext[I] += mesh_dof[mesh_l2g[eN_j]*3+1]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I];
814 }
815 }
816 //Space Mapping Jacobian
817 jac[XX] = Grad_x_ext[X];
818 jac[XY] = Grad_x_ext[Y];
819 jac[YX] = Grad_y_ext[X];
820 jac[YY] = Grad_y_ext[Y];
821 jacDet = jac[XX]*jac[YY] - jac[XY]*jac[YX];
822 oneOverJacDet = 1.0/jacDet;
823 jacInv[XX] = oneOverJacDet*jac[YY];
824 jacInv[XY] = -oneOverJacDet*jac[XY];
825 jacInv[YX] = -oneOverJacDet*jac[YX];
826 jacInv[YY] = oneOverJacDet*jac[XX];
827 //normal
828 norm_normal=0.0;
829 for (int I=0;I<2;I++)
830 normal[I] = 0.0;
831 for (int I=0;I<2;I++)
832 {
833 for (int J=0;J<2;J++)
834 {
835 normal[I] += jacInv[J*2+I]*normal_ref[ebN_local_kb_nSpace+J];
836 }
837 norm_normal+=normal[I]*normal[I];
838 }
839 norm_normal = sqrt(norm_normal);
840 for (int I=0;I<2;I++)
841 {
842 normal[I] /= norm_normal;
843 }
844 //metric tensor and determinant
845 boundaryJac[XHX] = jac[XX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHX]+jac[XY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHX];
846 boundaryJac[YHX] = jac[YX]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+XHX]+jac[YY]*boundaryJac_ref[ebN_local_kb_nSpace_nSpacem1+YHX];
847
848 metricTensor[HXHX] = boundaryJac[XHX]*boundaryJac[XHX]+boundaryJac[YHX]*boundaryJac[YHX];
849
850 metricTensorDetSqrt=sqrt(metricTensor[HXHX]);
851 }
852
854 const int ebN_local,
855 const int kb,
856 const int ebN_local_kb,
857 double* mesh_velocity_dof,
858 int* mesh_l2g,
859 double* mesh_trial_trace_ref,
860 double& xt,
861 double& yt,
862 double* normal,
863 double* boundaryJac,
864 double* metricTensor,
865 double& metricTensorDetSqrt)
866 {
867 //const int ebN_local_kb_nSpace = ebN_local_kb*2,
868 // ebN_local_kb_nSpace_nSpacem1 = ebN_local_kb*2*2;
869 //
870 //calculate velocity of mapping from the reference element to the physical element
871 //
872 xt=0.0;yt=0.0;
873 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
874 {
875 int eN_j = eN*NDOF_MESH_TRIAL_ELEMENT+j;
876 int ebN_local_kb_j = ebN_local_kb*NDOF_MESH_TRIAL_ELEMENT+j;
877 /* xt += mesh_velocity_dof[mesh_l2g[eN_j]*2+0]*mesh_trial_trace_ref[ebN_local_kb_j]; */
878 /* yt += mesh_velocity_dof[mesh_l2g[eN_j]*2+1]*mesh_trial_trace_ref[ebN_local_kb_j]; */
879 xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_trace_ref[ebN_local_kb_j];
880 yt += mesh_velocity_dof[mesh_l2g[eN_j]*3+1]*mesh_trial_trace_ref[ebN_local_kb_j];
881 }
882 //modify the metricTensorDetSqrt to include the effect of the moving domain
883 //it's not exactly the Sqrt(Det(G_tr_G)) now, see notes
884 //just do it brute force
885 double
886 Gy_tr_Gy_00 = 1.0 + xt*xt + yt*yt,
887 Gy_tr_Gy_01 = boundaryJac[XHX]*xt+boundaryJac[YHX]*yt,
888 Gy_tr_Gy_10 = Gy_tr_Gy_01,
889 Gy_tr_Gy_11 = metricTensor[HXHX],
890 xt_dot_n = xt*normal[X]+yt*normal[Y];
891 metricTensorDetSqrt=sqrt((Gy_tr_Gy_00*Gy_tr_Gy_11 - Gy_tr_Gy_01*Gy_tr_Gy_10) / (1.0+xt_dot_n*xt_dot_n));
892 }
893 inline void valFromDOF(const double* dof,const int* l2g_element,const double* trial_ref,double& val)
894 {
895 val=0.0;
896 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
897 val+=dof[l2g_element[j]]*trial_ref[j];
898 }
899
900 inline void gradFromDOF(const double* dof,const int* l2g_element,const double* grad_trial,double* grad)
901 {
902 for(int I=0;I<2;I++)
903 grad[I] = 0.0;
904 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
905 for(int I=0;I<2;I++)
906 grad[I] += dof[l2g_element[j]]*grad_trial[j*2+I];
907 }
908
909 inline void hessFromDOF(const double* dof,const int* l2g_element,const double* hess_trial,double* hess)
910 {
911 for(int I=0;I<2;I++)
912 for(int J=0;J<2;J++)
913 hess[I*2+J] = 0.0;
914 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
915 for(int I=0;I<2;I++)
916 for(int J=0;J<2;J++)
917 hess[I*2+J] += dof[l2g_element[j]]*hess_trial[j*4+I*2+J];
918 }
919
920 inline void calculateMapping_element(const int eN,
921 const int k,
922 double* mesh_dof,
923 int* mesh_l2g,
924 //xt::pyarray<double>& mesh_trial_ref,
925 double* mesh_trial_ref,
926 double* mesh_grad_trial_ref,
927 double* jac,
928 double& jacDet,
929 double* jacInv,
930 double& x,
931 double& y,
932 double& z)
933 {
935 k,
936 mesh_dof,
937 mesh_l2g,
938 mesh_trial_ref,
939 mesh_grad_trial_ref,
940 jac,
941 jacDet,
942 jacInv,
943 x,
944 y);
945 }
946 inline void calculateMappingVelocity_element(const int eN,
947 const int k,
948 double* mesh_velocity_dof,
949 int* mesh_l2g,
950 //xt::pyarray<double>& mesh_trial_ref,
951 double* mesh_trial_ref,
952 double& xt,
953 double& yt,
954 double& zt)
955 {
957 k,
958 mesh_velocity_dof,
959 mesh_l2g,
960 mesh_trial_ref,
961 xt,
962 yt);
963 }
964 inline void calculateMapping_elementBoundary(const int eN,
965 const int ebN_local,
966 const int kb,
967 const int ebN_local_kb,
968 double* mesh_dof,
969 int* mesh_l2g,
970 double* mesh_trial_trace_ref,
971 double* mesh_grad_trial_trace_ref,
972 double* boundaryJac_ref,
973 double* jac,
974 double& jacDet,
975 double* jacInv,
976 double* boundaryJac,
977 double* metricTensor,
978 double& metricTensorDetSqrt,
979 double* normal_ref,
980 double* normal,
981 double& x,
982 double& y,
983 double& z)
984 {
986 ebN_local,
987 kb,
988 ebN_local_kb,
989 mesh_dof,
990 mesh_l2g,
991 mesh_trial_trace_ref,
992 mesh_grad_trial_trace_ref,
993 boundaryJac_ref,
994 jac,
995 jacDet,
996 jacInv,
997 boundaryJac,
998 metricTensor,
999 metricTensorDetSqrt,
1000 normal_ref,
1001 normal,
1002 x,
1003 y);
1004 //A 2D mapping has no z component, but z is an OUTPUT parameter the caller will
1005 //read, so it must be assigned: leaving it alone hands back whatever the caller's
1006 //stack slot happened to hold.
1007 z=0.0;
1008 }
1010 const int ebN_local,
1011 const int kb,
1012 const int ebN_local_kb,
1013 double* mesh_velocity_dof,
1014 int* mesh_l2g,
1015 double* mesh_trial_trace_ref,
1016 double& xt,
1017 double& yt,
1018 double& zt,
1019 double* normal,
1020 double* boundaryJac,
1021 double* metricTensor,
1022 double& metricTensorDetSqrt)
1023 {
1025 ebN_local,
1026 kb,
1027 ebN_local_kb,
1028 mesh_velocity_dof,
1029 mesh_l2g,
1030 mesh_trial_trace_ref,
1031 xt,
1032 yt,
1033 normal,
1034 boundaryJac,
1035 metricTensor,
1036 metricTensorDetSqrt);
1037 }
1038};
1039
1040//specialization for 1D
1041template<int NDOF_MESH_TRIAL_ELEMENT>
1042class CompKernelSpaceMapping<1,NDOF_MESH_TRIAL_ELEMENT>
1043{
1044public:
1045 const int X,
1052 X(0),
1053 XX(0),
1054 sXX(0),
1055 nSymTen(1),
1056 XHX(0),
1057 HXHX(0)
1058 {}
1059 inline void calculateMapping_element(const int eN,
1060 const int k,
1061 double* mesh_dof,
1062 int* mesh_l2g,
1063 //xt::pyarray<double>& mesh_trial_ref,
1064 double* mesh_trial_ref,
1065 double* mesh_grad_trial_ref,
1066 double* jac,
1067 double& jacDet,
1068 double* jacInv,
1069 double& x)
1070 {
1071 double Grad_x[1],Grad_y[1],oneOverJacDet;
1072
1073 //
1074 //mapping of reference element to physical element
1075 //
1076 x=0.0;
1077 for (int I=0;I<1;I++)
1078 {
1079 Grad_x[I]=0.0;
1080 }
1081 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1082 {
1083 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
1084 x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
1085 for (int I=0;I<1;I++)
1086 {
1087 Grad_x[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j+I];
1088 }
1089 }
1090 jac[XX] = Grad_x[X];
1091 jacDet = jac[XX];
1092 oneOverJacDet = 1.0/jacDet;
1093 jacInv[XX] = oneOverJacDet;
1094 }
1095
1096 inline void calculateH_element(const int eN,
1097 const int k,
1098 double* h_dof,
1099 int* mesh_l2g,
1100 //xt::pyarray<double>& mesh_trial_ref,
1101 double* mesh_trial_ref,
1102 double& h)
1103 {
1104 h=0.0;
1105 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1106 {
1107 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
1108 h += h_dof[mesh_l2g[eN_j]]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
1109 }
1110 }
1111
1112 inline void calculateMappingVelocity_element(const int eN,
1113 const int k,
1114 double* mesh_velocity_dof,
1115 int* mesh_l2g,
1116 //xt::pyarray<double>& mesh_trial_ref,
1117 double* mesh_trial_ref,
1118 double& xt)
1119 {
1120 //
1121 //time derivative of mapping of reference element to physical element
1122 //
1123 xt=0.0;
1124 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1125 {
1126 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
1127 xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
1128 }
1129 }
1130
1131 inline void calculateMapping_elementBoundary(const int eN,
1132 const int ebN_local,
1133 const int kb,
1134 const int ebN_local_kb,
1135 double* mesh_dof,
1136 int* mesh_l2g,
1137 double* mesh_trial_trace_ref,
1138 double* mesh_grad_trial_trace_ref,
1139 double* boundaryJac_ref,
1140 double* jac,
1141 double& jacDet,
1142 double* jacInv,
1143 double* boundaryJac,
1144 double* metricTensor,
1145 double& metricTensorDetSqrt,
1146 double* normal_ref,
1147 double* normal,
1148 double& x)
1149 {
1150 const int ebN_local_kb_nSpace = ebN_local_kb*1;
1151
1152 double Grad_x_ext[1],oneOverJacDet,norm_normal=0.0;
1153 //
1154 //calculate mapping from the reference element to the physical element
1155 //
1156 x=0.0;
1157 for (int I=0;I<1;I++)
1158 {
1159 Grad_x_ext[I] = 0.0;
1160 }
1161 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1162 {
1163 int eN_j = eN*NDOF_MESH_TRIAL_ELEMENT+j;
1164 int ebN_local_kb_j = ebN_local_kb*NDOF_MESH_TRIAL_ELEMENT+j;
1165 int ebN_local_kb_j_nSpace = ebN_local_kb_j*1;
1166 x += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_trace_ref[ebN_local_kb_j];
1167 for (int I=0;I<1;I++)
1168 {
1169 Grad_x_ext[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I];
1170 }
1171 }
1172 //Space Mapping Jacobian
1173 jac[XX] = Grad_x_ext[X];
1174 jacDet = jac[XX];
1175 oneOverJacDet = 1.0/jacDet;
1176 jacInv[XX] = oneOverJacDet;
1177 //normal
1178 norm_normal=0.0;
1179 for (int I=0;I<1;I++)
1180 normal[I] = 0.0;
1181 for (int I=0;I<1;I++)
1182 {
1183 for (int J=0;J<1;J++)
1184 {
1185 normal[I] += jacInv[J+I]*normal_ref[ebN_local_kb_nSpace+J];
1186 }
1187 norm_normal+=normal[I]*normal[I];
1188 }
1189 norm_normal = sqrt(norm_normal);
1190 for (int I=0;I<1;I++)
1191 {
1192 normal[I] /= norm_normal;
1193 }
1194 //metric tensor and determinant
1195 boundaryJac[XHX] = 1.0;
1196
1197 metricTensor[HXHX] = 1.0;
1198
1199 metricTensorDetSqrt=1.0;
1200 }
1201
1203 const int ebN_local,
1204 const int kb,
1205 const int ebN_local_kb,
1206 double* mesh_velocity_dof,
1207 int* mesh_l2g,
1208 double* mesh_trial_trace_ref,
1209 double& xt,
1210 double* normal,
1211 double* boundaryJac,
1212 double* metricTensor,
1213 double& metricTensorDetSqrt)
1214 {
1215 //
1216 //calculate velocity of mapping from the reference element to the physical element
1217 //
1218 xt=0.0;
1219 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1220 {
1221 int eN_j = eN*NDOF_MESH_TRIAL_ELEMENT+j;
1222 int ebN_local_kb_j = ebN_local_kb*NDOF_MESH_TRIAL_ELEMENT+j;
1223 xt += mesh_velocity_dof[mesh_l2g[eN_j]*3+0]*mesh_trial_trace_ref[ebN_local_kb_j];
1224 }
1225 //modify the metricTensorDetSqrt to include the effect of the moving domain
1226 //it's not exactly the Sqrt(Det(G_tr_G)) now, see notes
1227 //just do it brute force
1228 double
1229 Gy_tr_Gy_00 = 1.0 + xt*xt,
1230 Gy_tr_Gy_01 = boundaryJac[XHX]*xt,
1231 Gy_tr_Gy_10 = Gy_tr_Gy_01,
1232 Gy_tr_Gy_11 = metricTensor[HXHX],
1233 xt_dot_n = xt*normal[X];
1234 metricTensorDetSqrt=sqrt((Gy_tr_Gy_00*Gy_tr_Gy_11 - Gy_tr_Gy_01*Gy_tr_Gy_10) / (1.0+xt_dot_n*xt_dot_n));
1235 }
1236 inline void valFromDOF(const double* dof,const int* l2g_element,const double* trial_ref,double& val)
1237 {
1238 val=0.0;
1239 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1240 val+=dof[l2g_element[j]]*trial_ref[j];
1241 }
1242
1243 inline void gradFromDOF(const double* dof,const int* l2g_element,const double* grad_trial,double* grad)
1244 {
1245 for(int I=0;I<1;I++)
1246 grad[I] = 0.0;
1247 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1248 for(int I=0;I<1;I++)
1249 grad[I] += dof[l2g_element[j]]*grad_trial[j+I];
1250 }
1251
1252 inline void hessFromDOF(const double* dof,const int* l2g_element,const double* hess_trial,double* hess)
1253 {
1254 for(int I=0;I<1;I++)
1255 for(int J=0;J<1;J++)
1256 hess[I+J] = 0.0;
1257 for (int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1258 for(int I=0;I<1;I++)
1259 for(int J=0;J<1;J++)
1260 hess[I+J] += dof[l2g_element[j]]*hess_trial[j+I+J];
1261 }
1262
1263 inline void calculateMapping_element(const int eN,
1264 const int k,
1265 double* mesh_dof,
1266 int* mesh_l2g,
1267 //xt::pyarray<double>& mesh_trial_ref,
1268 double* mesh_trial_ref,
1269 double* mesh_grad_trial_ref,
1270 double* jac,
1271 double& jacDet,
1272 double* jacInv,
1273 double& x,
1274 double& y,
1275 double& z)
1276 {
1278 k,
1279 mesh_dof,
1280 mesh_l2g,
1281 mesh_trial_ref,
1282 mesh_grad_trial_ref,
1283 jac,
1284 jacDet,
1285 jacInv,
1286 x);
1287 }
1288 inline void calculateMappingVelocity_element(const int eN,
1289 const int k,
1290 double* mesh_velocity_dof,
1291 int* mesh_l2g,
1292 //xt::pyarray<double>& mesh_trial_ref,
1293 double* mesh_trial_ref,
1294 double& xt,
1295 double& yt,
1296 double& zt)
1297 {
1299 k,
1300 mesh_velocity_dof,
1301 mesh_l2g,
1302 mesh_trial_ref,
1303 xt);
1304 }
1305 inline void calculateMapping_elementBoundary(const int eN,
1306 const int ebN_local,
1307 const int kb,
1308 const int ebN_local_kb,
1309 double* mesh_dof,
1310 int* mesh_l2g,
1311 double* mesh_trial_trace_ref,
1312 double* mesh_grad_trial_trace_ref,
1313 double* boundaryJac_ref,
1314 double* jac,
1315 double& jacDet,
1316 double* jacInv,
1317 double* boundaryJac,
1318 double* metricTensor,
1319 double& metricTensorDetSqrt,
1320 double* normal_ref,
1321 double* normal,
1322 double& x,
1323 double& y,
1324 double& z)
1325 {
1327 ebN_local,
1328 kb,
1329 ebN_local_kb,
1330 mesh_dof,
1331 mesh_l2g,
1332 mesh_trial_trace_ref,
1333 mesh_grad_trial_trace_ref,
1334 boundaryJac_ref,
1335 jac,
1336 jacDet,
1337 jacInv,
1338 boundaryJac,
1339 metricTensor,
1340 metricTensorDetSqrt,
1341 normal_ref,
1342 normal,
1343 x);
1344 //likewise for a 1D mapping: y and z have no counterpart in the mapping but are
1345 //outputs the caller will read.
1346 y=0.0;
1347 z=0.0;
1348 }
1350 const int ebN_local,
1351 const int kb,
1352 const int ebN_local_kb,
1353 double* mesh_velocity_dof,
1354 int* mesh_l2g,
1355 double* mesh_trial_trace_ref,
1356 double& xt,
1357 double& yt,
1358 double& zt,
1359 double* normal,
1360 double* boundaryJac,
1361 double* metricTensor,
1362 double& metricTensorDetSqrt)
1363 {
1365 ebN_local,
1366 kb,
1367 ebN_local_kb,
1368 mesh_velocity_dof,
1369 mesh_l2g,
1370 mesh_trial_trace_ref,
1371 xt,
1372 normal,
1373 boundaryJac,
1374 metricTensor,
1375 metricTensorDetSqrt);
1376 }
1377};
1378
1379
1380template<int NSPACE, int NDOF_MESH_TRIAL_ELEMENT, int NDOF_TRIAL_ELEMENT, int NDOF_TEST_ELEMENT>
1382{
1383public:
1385 const int X,
1417 inline void calculateG(double* jacInv,double* G,double& G_dd_G, double& tr_G)
1418 {
1419 //get the metric tensor
1420 //cek todo use symmetry
1421 for (int I=0;I<NSPACE;I++)
1422 for (int J=0;J<NSPACE;J++)
1423 {
1424 G[I*NSPACE+J] = 0.0;
1425 for (int K=0;K<NSPACE;K++)
1426 G[I*NSPACE+J] += jacInv[K*NSPACE+I]*jacInv[K*NSPACE+J];
1427 }
1428 G_dd_G = 0.0;
1429 tr_G = 0.0;
1430 for (int I=0;I<NSPACE;I++)
1431 {
1432 tr_G += G[I*NSPACE+I];
1433 for (int J=0;J<NSPACE;J++)
1434 {
1435 G_dd_G += G[I*NSPACE+J]*G[I*NSPACE+J];
1436 }
1437 }
1438 }
1439 inline void calculateGScale(double* G,double* v,double& h)
1440 {
1441 h = 0.0;
1442 for (int I=0;I<NSPACE;I++)
1443 for (int J=0;J<NSPACE;J++)
1444 h += v[I]*G[I*NSPACE+J]*v[J];
1445 h = 1.0/sqrt(h+1.0e-16);//cek hack
1446 }
1447 inline void valFromDOF(const double* dof,const int* l2g_element,const double* trial_ref,double& val)
1448 {
1449 val=0.0;
1450 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1451 val+=dof[l2g_element[j]]*trial_ref[j];
1452 }
1453
1454 inline void gradFromDOF(const double* dof,const int* l2g_element,const double* grad_trial,double* grad)
1455 {
1456 for(int I=0;I<NSPACE;I++)
1457 grad[I] = 0.0;
1458 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1459 for(int I=0;I<NSPACE;I++)
1460 grad[I] += dof[l2g_element[j]]*grad_trial[j*NSPACE+I];
1461 }
1462
1463 inline void hessFromDOF(const double* dof,const int* l2g_element,const double* hess_trial,double* hess)
1464 {
1465 const int NSPACE2=NSPACE*NSPACE;
1466 for(int I=0;I<NSPACE;I++)
1467 for(int J=0;J<NSPACE;J++)
1468 hess[I*NSPACE+J] = 0.0;
1469 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1470 for(int I=0;I<NSPACE;I++)
1471 for(int J=0;J<NSPACE;J++)
1472 hess[I*NSPACE+J] += dof[l2g_element[j]]*hess_trial[j*NSPACE2+I*NSPACE+J];
1473 }
1474
1475 inline void valFromElementDOF(const double* dof,const double* trial_ref,double& val)
1476 {
1477 val=0.0;
1478 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1479 val+=dof[j]*trial_ref[j];
1480 }
1481
1482 inline void gradFromElementDOF(const double* dof,const double* grad_trial,double* grad)
1483 {
1484 for(int I=0;I<NSPACE;I++)
1485 grad[I] = 0.0;
1486 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1487 for(int I=0;I<NSPACE;I++)
1488 grad[I] += dof[j]*grad_trial[j*NSPACE+I];
1489 }
1490
1491 inline void gradTrialFromRef(const double* grad_trial_ref, const double* jacInv, double* grad_trial)
1492 {
1493 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1494 for(int I=0;I<NSPACE;I++)
1495 grad_trial[j*NSPACE+I] = 0.0;
1496 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1497 for(int I=0;I<NSPACE;I++)
1498 for(int J=0;J<NSPACE;J++)
1499 grad_trial[j*NSPACE+I] += jacInv[J*NSPACE+I]*grad_trial_ref[j*NSPACE+J];
1500 }
1501
1502 /*
1503 * DOFaverage
1504 * ----------
1505 *
1506 * Calculate the average DOF value for at a given mesh element.
1507 *
1508 * @param dof array of finite element DOF values
1509 * @param l2g_element local 2 global mapping for the current mesh element
1510 * @param val return value with the average DOF values
1511 */
1512
1513 inline void DOFaverage (const double* dof, const int* l2g_element, double& val)
1514 {
1515 val = 0.0;
1516
1517 for (int j=0; j<NDOF_MESH_TRIAL_ELEMENT; j++)
1518 val+=dof[l2g_element[j]];
1519
1520 val /= NDOF_MESH_TRIAL_ELEMENT;
1521 }
1522
1523 inline void hessTrialFromRef(const double* hess_trial_ref, const double* jacInv, double* hess_trial)
1524 {
1525 const int NSPACE2=NSPACE*NSPACE;
1526 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1527 for(int I=0;I<NSPACE;I++)
1528 for(int J=0;J<NSPACE;J++)
1529 hess_trial[j*NSPACE2+I*NSPACE+J] = 0.0;
1530 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1531 for(int I=0;I<NSPACE;I++)
1532 for(int J=0;J<NSPACE;J++)
1533 for(int K=0;K<NSPACE;K++)
1534 for(int L=0;L<NSPACE;L++)
1535 hess_trial[j*NSPACE2+I*NSPACE+J] += hess_trial_ref[j*NSPACE2+K*NSPACE+L]*jacInv[L*NSPACE+J]*jacInv[K*NSPACE+I];
1536 }
1537
1538 inline void gradTestFromRef(const double* grad_test_ref, const double* jacInv, double* grad_test)
1539 {
1540 for (int i=0;i<NDOF_TEST_ELEMENT;i++)
1541 for(int I=0;I<NSPACE;I++)
1542 grad_test[i*NSPACE+I] = 0.0;
1543 for (int i=0;i<NDOF_TEST_ELEMENT;i++)
1544 for(int I=0;I<NSPACE;I++)
1545 for(int J=0;J<NSPACE;J++)
1546 grad_test[i*NSPACE+I] += jacInv[J*NSPACE+I]*grad_test_ref[i*NSPACE+J];
1547 }
1548
1549 inline void backwardEuler(const double& dt, const double& m_old, const double& m, const double& dm, double& mt, double& dmt)
1550 {
1551 mt =(m-m_old)/dt;
1552 dmt = dm/dt;
1553 }
1554
1555 inline void bdf(const double& alpha, const double& beta, const double& m, const double& dm, double& mt, double& dmt)
1556 {
1557 mt =alpha*m + beta;
1558 dmt = alpha*dm;
1559 }
1560 inline void bdfC2(const double& alpha, const double& beta, const double& m, const double& dm, const double& dm2, double& mt, double& dmt, double& dm2t)
1561 {
1562 mt =alpha*m + beta;
1563 dmt = alpha*dm;
1564 dm2t = alpha*dm2;
1565 }
1566
1567 inline double Mass_weak(const double& mt, const double& w_dV)
1568 {
1569 return mt*w_dV;
1570 }
1571
1572 inline double MassJacobian_weak(const double& dmt,
1573 const double& v,
1574 const double& w_dV)
1575 {
1576 return dmt*v*w_dV;
1577 }
1578
1579 inline double Mass_strong(const double& mt)
1580 {
1581 return mt;
1582 }
1583
1584 inline double MassJacobian_strong(const double& dmt,
1585 const double& v)
1586 {
1587 return dmt*v;
1588 }
1589
1590 inline double Mass_adjoint(const double& dmt,
1591 const double& w_dV)
1592 {
1593 return dmt*w_dV;
1594 }
1595
1596 /*
1597 * pressureProjection_weak
1598 * -----------------------
1599 *
1600 * Inner product calculation of the pressure projection
1601 * stablization method of Bochev, Dohrmann and
1602 * Gunzburger (2006).
1603 *
1604 * @param viscosity viscosity at point
1605 * @param p the pressure value (either trial function
1606 * or actual value)
1607 * @param p_avg the pressure projection value (either
1608 * 1./3. for test functions of the average
1609 * value of the pressure on the element)
1610 * @param dV this is the integral weight
1611 *
1612 */
1613
1614 inline double pressureProjection_weak(const double& viscosity,
1615 const double& p,
1616 const double& p_avg,
1617 const double& q,
1618 const double& dV)
1619 {
1620 if (viscosity==0.){ return 0.;}
1621 return (1./viscosity)*(p-p_avg)*(q-1./3.)*dV;
1622 }
1623
1624 inline double Advection_weak(const double f[NSPACE],
1625 const double grad_w_dV[NSPACE])
1626 {
1627 double tmp=0.0;
1628 for(int I=0;I<NSPACE;I++)
1629 tmp -= f[I]*grad_w_dV[I];
1630 return tmp;
1631 }
1632
1633 inline double AdvectionJacobian_weak(const double df[NSPACE],
1634 const double& v,
1635 const double grad_w_dV[NSPACE])
1636 {
1637 double tmp=0.0;
1638 for(int I=0;I<NSPACE;I++)
1639 tmp -= df[I]*v*grad_w_dV[I];
1640 return tmp;
1641 }
1642
1643 inline double Advection_strong(const double df[NSPACE],
1644 const double grad_u[NSPACE])
1645 {
1646 double tmp=0.0;
1647 for(int I=0;I<NSPACE;I++)
1648 tmp += df[I]*grad_u[I];
1649 return tmp;
1650 }
1651
1652 inline double AdvectionJacobian_strong(const double df[NSPACE],
1653 const double grad_v[NSPACE])
1654 {
1655 double tmp=0.0;
1656 for(int I=0;I<NSPACE;I++)
1657 tmp += df[I]*grad_v[I];
1658 return tmp;
1659 }
1660
1661 inline double Advection_adjoint(const double df[NSPACE],
1662 const double grad_w_dV[NSPACE])
1663 {
1664 double tmp=0.0;
1665 for(int I=0;I<NSPACE;I++)
1666 tmp -= df[I]*grad_w_dV[I];
1667 return tmp;
1668 }
1669
1670 inline double Hamiltonian_weak(const double& H,
1671 const double& w_dV)
1672 {
1673 return H*w_dV;
1674 }
1675
1676 inline double HamiltonianJacobian_weak(const double dH[NSPACE],
1677 const double grad_v[NSPACE],
1678 const double& w_dV)
1679 {
1680 double tmp=0.0;
1681 for(int I=0;I<NSPACE;I++)
1682 tmp += dH[I]*grad_v[I]*w_dV;
1683 return tmp;
1684 }
1685
1686 inline double Hamiltonian_strong(const double dH[NSPACE],
1687 const double grad_u[NSPACE])
1688 {
1689 double tmp=0.0;
1690 for(int I=0;I<NSPACE;I++)
1691 tmp += dH[I]*grad_u[I];
1692 return tmp;
1693 }
1694
1695 inline double HamiltonianJacobian_strong(const double dH[NSPACE],
1696 const double grad_v[NSPACE])
1697 {
1698 double tmp=0.0;
1699 for(int I=0;I<NSPACE;I++)
1700 tmp += dH[I]*grad_v[I];
1701 return tmp;
1702 }
1703
1704 inline double Hamiltonian_adjoint(const double dH[NSPACE],
1705 const double grad_w_dV[NSPACE])
1706 {
1707 double tmp=0.0;
1708 for(int I=0;I<NSPACE;I++)
1709 tmp -= dH[I]*grad_w_dV[I];
1710 return tmp;
1711 }
1712
1713 inline double Diffusion_weak(int* rowptr,
1714 int* colind,
1715 double* a,
1716 const double grad_phi[NSPACE],
1717 const double grad_w_dV[NSPACE])
1718 {
1719 double tmp=0.0;
1720 for(int I=0;I<NSPACE;I++)
1721 for (int m=rowptr[I];m<rowptr[I+1];m++)
1722 tmp += a[m]*grad_phi[colind[m]]*grad_w_dV[I];
1723 return tmp;
1724 }
1725
1726 inline double DiffusionJacobian_weak(int* rowptr,
1727 int* colind,
1728 double* a,
1729 double* da,
1730 const double grad_phi[NSPACE],
1731 const double grad_w_dV[NSPACE],
1732 const double& dphi,
1733 const double& v,
1734 const double grad_v[NSPACE])
1735 {
1736 double daProduct=0.0,dphiProduct=0.0;
1737 for (int I=0;I<NSPACE;I++)
1738 for (int m=rowptr[I];m<rowptr[I+1];m++)
1739 {
1740 daProduct += da[m]*grad_phi[colind[m]]*grad_w_dV[I];
1741 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
1742 }
1743 return daProduct*v+dphiProduct*dphi;
1744 }
1745
1746 inline double SimpleDiffusionJacobian_weak(int* rowptr,
1747 int* colind,
1748 double* a,
1749 const double grad_v[NSPACE],
1750 const double grad_w_dV[NSPACE])
1751 {
1752 double dphiProduct=0.0;
1753 for (int I=0;I<NSPACE;I++)
1754 for (int m=rowptr[I];m<rowptr[I+1];m++)
1755 {
1756 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
1757 }
1758 return dphiProduct;
1759 }
1760
1761 inline double Reaction_weak(const double& r,
1762 const double& w_dV)
1763 {
1764 return r*w_dV;
1765 }
1766
1767 inline double ReactionJacobian_weak(const double& dr,
1768 const double& v,
1769 const double& w_dV)
1770 {
1771 return dr*v*w_dV;
1772 }
1773
1774 inline double Reaction_strong(const double& r)
1775 {
1776 return r;
1777 }
1778
1779 inline double ReactionJacobian_strong(const double& dr,
1780 const double& v)
1781 {
1782 return dr*v;
1783 }
1784
1785 inline double Reaction_adjoint(const double& dr,
1786 const double& w_dV)
1787 {
1788 return dr*w_dV;
1789 }
1790
1791 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
1792 const double& elementDiameter,
1793 const double& strong_residual,
1794 const double grad_u[NSPACE],
1795 double& numDiff)
1796 {
1797 double h,
1798 num,
1799 den,
1800 n_grad_u;
1801 h = elementDiameter;
1802 n_grad_u = 0.0;
1803 for (int I=0;I<NSPACE;I++)
1804 n_grad_u += grad_u[I]*grad_u[I];
1805 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
1806 den = sqrt(n_grad_u+1.0e-12);
1807 numDiff = num/den;
1808 }
1809
1810
1811 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
1812 const double G[NSPACE*NSPACE],
1813 const double& strong_residual,
1814 const double grad_u[NSPACE],
1815 double& numDiff)
1816 {
1817 double den = 0.0;
1818 for (int I=0;I<NSPACE;I++)
1819 for (int J=0;J<NSPACE;J++)
1820 den += grad_u[I]*G[I*NSPACE+J]*grad_u[J];
1821 numDiff = shockCapturingDiffusion*fabs(strong_residual)/(sqrt(den+1.0e-12));
1822 }
1823
1824 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
1825 const double& uref, const double& beta,
1826 const double G[NSPACE*NSPACE],
1827 const double& G_dd_G,
1828 const double& strong_residual,
1829 const double grad_u[NSPACE],
1830 double& numDiff)
1831 {
1832 double den = 0.0;
1833 for (int I=0;I<NSPACE;I++)
1834 for (int J=0;J<NSPACE;J++)
1835 den += grad_u[I]*G[I*NSPACE+J]*grad_u[J];
1836
1837 double h2_uref_1 = 1.0/(sqrt(den+1.0e-12));
1838 double h2_uref_2 = 1.0/(uref*sqrt(G_dd_G+1.0e-12));
1839 numDiff = shockCapturingDiffusion*fabs(strong_residual)*pow(h2_uref_1, 2.0-beta)*pow(h2_uref_2,beta-1.0);
1840 }
1841
1842
1843
1844 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
1845 const double G[NSPACE*NSPACE],
1846 const double& strong_residual,
1847 const double vel[NSPACE],
1848 const double grad_u[NSPACE],
1849 double& numDiff)
1850 {
1851 double den1 = 0.0,den2=0.0, nom=0.0;
1852 for (int I=0;I<NSPACE;I++)
1853 {
1854 nom += vel[I]*vel[I];
1855 den2+= grad_u[I]*grad_u[I];
1856 for (int J=0;J<NSPACE;J++)
1857 den1 += vel[I]*G[I*NSPACE+J]*vel[J];
1858 }
1859 numDiff = shockCapturingDiffusion*fabs(strong_residual)*(sqrt(nom/(den1*den2 + 1.0e-12)));
1860 }
1861
1862
1863 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
1864 const double& elementDiameter,
1865 const double& strong_residual,
1866 const double grad_u[NSPACE],
1867 double& gradNorm,
1868 double& gradNorm_last,
1869 double& numDiff)
1870 {
1871 double h,
1872 num,
1873 n_grad_u;
1874 h = elementDiameter;
1875 n_grad_u = 0.0;
1876 for (int I=0;I<NSPACE;I++)
1877 n_grad_u += grad_u[I]*grad_u[I];
1878 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
1879 gradNorm = sqrt(n_grad_u+1.0e-12);
1880 //cek hack shockCapturingDiffusion*fabs(strong_residual)*grad_phi_G_grad_phi
1881 numDiff = num/gradNorm_last;
1882 }
1883
1884 inline double SubgridError(const double& error,
1885 const double& Lstar_w_dV)
1886 {
1887 return error*Lstar_w_dV;
1888 }
1889
1890 inline double SubgridErrorJacobian(const double& derror,
1891 const double& Lstar_w_dV)
1892 {
1893 return derror*Lstar_w_dV;
1894 }
1895
1896 inline double NumericalDiffusion(const double& numDiff,
1897 const double grad_u[NSPACE],
1898 const double grad_w_dV[NSPACE])
1899 {
1900 double tmp=0.0;
1901 for (int I=0;I<NSPACE;I++)
1902 tmp += numDiff*grad_u[I]*grad_w_dV[I];
1903 return tmp;
1904 }
1905
1906 inline double NumericalDiffusionJacobian(const double& numDiff,
1907 const double grad_v[NSPACE],
1908 const double grad_w_dV[NSPACE])
1909 {
1910 double tmp=0.0;
1911 for (int I=0;I<NSPACE;I++)
1912 tmp += numDiff*grad_v[I]*grad_w_dV[I];
1913 return tmp;
1914 }
1915
1916
1917
1918 inline double ExteriorElementBoundaryFlux(const double& flux,
1919 const double& w_dS)
1920 {
1921 return flux*w_dS;
1922 }
1923
1924 inline double InteriorElementBoundaryFlux(const double& flux,
1925 const double& w_dS)
1926 {
1927 return flux*w_dS;
1928 }
1929
1930 inline double ExteriorNumericalAdvectiveFluxJacobian(const double& dflux_left,
1931 const double& v)
1932 {
1933 return dflux_left*v;
1934 }
1935
1936 inline double InteriorNumericalAdvectiveFluxJacobian(const double& dflux_left,
1937 const double& v)
1938 {
1939 return dflux_left*v;
1940 }
1941
1942 inline double ExteriorElementBoundaryScalarDiffusionAdjoint(const int& isDOFBoundary,
1943 const int& isFluxBoundary,
1944 const double& sigma,
1945 const double& u,
1946 const double& bc_u,
1947 const double normal[NSPACE],
1948 const double& a,
1949 const double grad_w_dS[NSPACE])
1950 {
1951 double tmp=0.0;
1952 for(int I=0;I<NSPACE;I++)
1953 {
1954 tmp += normal[I]*grad_w_dS[I];
1955 }
1956 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*(u-bc_u)*a;
1957 return tmp;
1958 }
1959
1960 inline double ExteriorElementBoundaryScalarDiffusionAdjointJacobian(const int& isDOFBoundary,
1961 const int& isFluxBoundary,
1962 const double& sigma,
1963 const double& v,
1964 const double normal[NSPACE],
1965 const double& a,
1966 const double grad_w_dS[NSPACE])
1967 {
1968 double tmp=0.0;
1969 for(int I=0;I<NSPACE;I++)
1970 {
1971 tmp += normal[I]*grad_w_dS[I];
1972 }
1973 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*v*a;
1974 return tmp;
1975 }
1976
1977 inline double ExteriorElementBoundaryDiffusionAdjoint(const int& isDOFBoundary,
1978 const int& isFluxBoundary,
1979 const double& sigma,
1980 const double& u,
1981 const double& bc_u,
1982 const double normal[NSPACE],
1983 int* rowptr,
1984 int* colind,
1985 double* a,
1986 const double grad_w_dS[NSPACE])
1987 {
1988 double tmp=0.0;
1989 for(int I=0;I<NSPACE;I++)
1990 for (int m=rowptr[I];m<rowptr[I+1];m++)
1991 tmp += (1.0-isFluxBoundary)*isDOFBoundary*sigma*(u-bc_u)*a[m]*normal[colind[m]]*grad_w_dS[I];
1992 return tmp;
1993 }
1994
1995 inline double ExteriorElementBoundaryDiffusionAdjointJacobian(const int& isDOFBoundary,
1996 const int& isFluxBoundary,
1997 const double& sigma,
1998 const double& v,
1999 const double normal[NSPACE],
2000 int* rowptr,
2001 int* colind,
2002 double* a,
2003 const double grad_w_dS[NSPACE])
2004 {
2005 double tmp=0.0;
2006 for(int I=0;I<NSPACE;I++)
2007 for (int m=rowptr[I];m<rowptr[I+1];m++)
2008 tmp += (1.0-isFluxBoundary)*isDOFBoundary*sigma*v*a[m]*normal[colind[m]]*grad_w_dS[I];
2009 return tmp;
2010 }
2011
2012 inline void calculateMapping_element(const int eN,
2013 const int k,
2014 double* mesh_dof,
2015 int* mesh_l2g,
2016 //xt::pyarray<double>& mesh_trial_ref,
2017 double* mesh_trial_ref,
2018 double* mesh_grad_trial_ref,
2019 double* jac,
2020 double& jacDet,
2021 double* jacInv,
2022 double& x,
2023 double& y,
2024 double& z)
2025 {
2026 mapping.calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y,z);
2027 }
2028
2029 inline void calculateH_element(const int eN,
2030 const int k,
2031 double* h_dof,
2032 int* mesh_l2g,
2033 //xt::pyarray<double>& mesh_trial_ref,
2034 double* mesh_trial_ref,
2035 double& h)
2036 {
2038 k,
2039 h_dof,
2040 mesh_l2g,
2041 mesh_trial_ref,
2042 h);
2043 }
2044
2045 inline void calculateMappingVelocity_element(const int eN,
2046 const int k,
2047 double* meshVelocity_dof,
2048 int* mesh_l2g,
2049 //xt::pyarray<double>& mesh_trial_ref,
2050 double* mesh_trial_ref,
2051 double& xt,
2052 double& yt,
2053 double& zt)
2054 {
2055 mapping.calculateMappingVelocity_element(eN,k,meshVelocity_dof,mesh_l2g,mesh_trial_ref,xt,yt,zt);
2056 }
2057
2058 inline
2060 const int ebN_local,
2061 const int kb,
2062 const int ebN_local_kb,
2063 double* mesh_dof,
2064 int* mesh_l2g,
2065 double* mesh_trial_trace_ref,
2066 double* mesh_grad_trial_trace_ref,
2067 double* boundaryJac_ref,
2068 double* jac,
2069 double& jacDet,
2070 double* jacInv,
2071 double* boundaryJac,
2072 double* metricTensor,
2073 double& metricTensorDetSqrt,
2074 double* normal_ref,
2075 double* normal,
2076 double& x,
2077 double& y,
2078 double& z)
2079 {
2081 ebN_local,
2082 kb,
2083 ebN_local_kb,
2084 mesh_dof,
2085 mesh_l2g,
2086 mesh_trial_trace_ref,
2087 mesh_grad_trial_trace_ref,
2088 boundaryJac_ref,
2089 jac,
2090 jacDet,
2091 jacInv,
2092 boundaryJac,
2093 metricTensor,
2094 metricTensorDetSqrt,
2095 normal_ref,
2096 normal,
2097 x,
2098 y,
2099 z);
2100 }
2101
2102 inline
2104 const int ebN_local,
2105 const int kb,
2106 const int ebN_local_kb,
2107 double* mesh_velocity_dof,
2108 int* mesh_l2g,
2109 double* mesh_trial_trace_ref,
2110 double& xt,
2111 double& yt,
2112 double& zt,
2113 double* normal,
2114 double* boundaryJac,
2115 double* metricTensor,
2116 double& metricTensorDetSqrt)
2117 {
2119 ebN_local,
2120 kb,
2121 ebN_local_kb,
2122 mesh_velocity_dof,
2123 mesh_l2g,
2124 mesh_trial_trace_ref,
2125 xt,
2126 yt,
2127 zt,
2128 normal,
2129 boundaryJac,
2130 metricTensor,
2131 metricTensorDetSqrt);
2132 }
2133 double Stress_u_weak(double* stress, double* grad_test_dV)
2134 {
2135 return stress[sXX]*grad_test_dV[X] + stress[sXY]*grad_test_dV[Y] + stress[sXZ]*grad_test_dV[Z];
2136 }
2137 double StressJacobian_u_u_weak(double* dstress, double* grad_trial, double* grad_test_dV)
2138 {
2139 return
2140 (dstress[sXX*nSymTen+sXX]*grad_trial[X]+dstress[sXX*nSymTen+sXY]*grad_trial[Y]+dstress[sXX*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[X] +
2141 (dstress[sXY*nSymTen+sXX]*grad_trial[X]+dstress[sXY*nSymTen+sXY]*grad_trial[Y]+dstress[sXY*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[Y] +
2142 (dstress[sXZ*nSymTen+sXX]*grad_trial[X]+dstress[sXZ*nSymTen+sXY]*grad_trial[Y]+dstress[sXZ*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[Z];
2143 }
2144 double StressJacobian_u_v_weak(double* dstress, double* grad_trial, double* grad_test_dV)
2145 {
2146 return
2147 (dstress[sXX*nSymTen+sYX]*grad_trial[X]+dstress[sXX*nSymTen+sYY]*grad_trial[Y]+dstress[sXX*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[X] +
2148 (dstress[sXY*nSymTen+sYX]*grad_trial[X]+dstress[sXY*nSymTen+sYY]*grad_trial[Y]+dstress[sXY*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[Y] +
2149 (dstress[sXZ*nSymTen+sYX]*grad_trial[X]+dstress[sXZ*nSymTen+sYY]*grad_trial[Y]+dstress[sXZ*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[Z];
2150 }
2151 double StressJacobian_u_w_weak(double* dstress, double* grad_trial, double* grad_test_dV)
2152 {
2153 return
2154 (dstress[sXX*nSymTen+sZX]*grad_trial[X]+dstress[sXX*nSymTen+sZY]*grad_trial[Y]+dstress[sXX*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[X] +
2155 (dstress[sXY*nSymTen+sZX]*grad_trial[X]+dstress[sXY*nSymTen+sZY]*grad_trial[Y]+dstress[sXY*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[Y] +
2156 (dstress[sXZ*nSymTen+sZX]*grad_trial[X]+dstress[sXZ*nSymTen+sZY]*grad_trial[Y]+dstress[sXZ*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[Z];
2157 }
2158 double Stress_v_weak(double* stress, double* grad_test_dV)
2159 {
2160 return stress[sYX]*grad_test_dV[X] + stress[sYY]*grad_test_dV[Y] + stress[sYZ]*grad_test_dV[Z];
2161 }
2162 double StressJacobian_v_u_weak(double* dstress, double* grad_trial,double* grad_test_dV)
2163 {
2164 return
2165 (dstress[sYX*nSymTen+sXX]*grad_trial[X]+dstress[sYX*nSymTen+sXY]*grad_trial[Y]+dstress[sYX*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[X] +
2166 (dstress[sYY*nSymTen+sXX]*grad_trial[X]+dstress[sYY*nSymTen+sXY]*grad_trial[Y]+dstress[sYY*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[Y] +
2167 (dstress[sYZ*nSymTen+sXX]*grad_trial[X]+dstress[sYZ*nSymTen+sXY]*grad_trial[Y]+dstress[sYZ*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[Z];
2168 }
2169 double StressJacobian_v_v_weak(double* dstress, double* grad_trial,double* grad_test_dV)
2170 {
2171 return
2172 (dstress[sYX*nSymTen+sYX]*grad_trial[X]+dstress[sYX*nSymTen+sYY]*grad_trial[Y]+dstress[sYX*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[X] +
2173 (dstress[sYY*nSymTen+sYX]*grad_trial[X]+dstress[sYY*nSymTen+sYY]*grad_trial[Y]+dstress[sYY*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[Y] +
2174 (dstress[sYZ*nSymTen+sYX]*grad_trial[X]+dstress[sYZ*nSymTen+sYY]*grad_trial[Y]+dstress[sYZ*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[Z];
2175 }
2176 double StressJacobian_v_w_weak(double* dstress, double* grad_trial,double* grad_test_dV)
2177 {
2178 return
2179 (dstress[sYX*nSymTen+sZX]*grad_trial[X]+dstress[sYX*nSymTen+sZY]*grad_trial[Y]+dstress[sYX*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[X] +
2180 (dstress[sYY*nSymTen+sZX]*grad_trial[X]+dstress[sYY*nSymTen+sZY]*grad_trial[Y]+dstress[sYY*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[Y] +
2181 (dstress[sYZ*nSymTen+sZX]*grad_trial[X]+dstress[sYZ*nSymTen+sZY]*grad_trial[Y]+dstress[sYZ*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[Z];
2182 }
2183 double Stress_w_weak(double* stress, double* grad_test_dV)
2184 {
2185 return stress[sZX]*grad_test_dV[X] + stress[sZY]*grad_test_dV[Y] + stress[sZZ]*grad_test_dV[Z];
2186 }
2187 double StressJacobian_w_u_weak(double* dstress, double* grad_trial, double* grad_test_dV)
2188 {
2189 return
2190 (dstress[sZX*nSymTen+sXX]*grad_trial[X]+dstress[sZX*nSymTen+sXY]*grad_trial[Y]+dstress[sZX*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[X] +
2191 (dstress[sZY*nSymTen+sXX]*grad_trial[X]+dstress[sZY*nSymTen+sXY]*grad_trial[Y]+dstress[sZY*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[Y] +
2192 (dstress[sZZ*nSymTen+sXX]*grad_trial[X]+dstress[sZZ*nSymTen+sXY]*grad_trial[Y]+dstress[sZZ*nSymTen+sXZ]*grad_trial[Z])*grad_test_dV[Z];
2193 }
2194 double StressJacobian_w_v_weak(double* dstress, double* grad_trial, double* grad_test_dV)
2195 {
2196 return
2197 (dstress[sZX*nSymTen+sYX]*grad_trial[X]+dstress[sZX*nSymTen+sYY]*grad_trial[Y]+dstress[sZX*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[X] +
2198 (dstress[sZY*nSymTen+sYX]*grad_trial[X]+dstress[sZY*nSymTen+sYY]*grad_trial[Y]+dstress[sZY*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[Y] +
2199 (dstress[sZZ*nSymTen+sYX]*grad_trial[X]+dstress[sZZ*nSymTen+sYY]*grad_trial[Y]+dstress[sZZ*nSymTen+sYZ]*grad_trial[Z])*grad_test_dV[Z];
2200 }
2201 double StressJacobian_w_w_weak(double* dstress, double* grad_trial, double* grad_test_dV)
2202 {
2203 return
2204 (dstress[sZX*nSymTen+sZX]*grad_trial[X]+dstress[sZX*nSymTen+sZY]*grad_trial[Y]+dstress[sZX*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[X] +
2205 (dstress[sZY*nSymTen+sZX]*grad_trial[X]+dstress[sZY*nSymTen+sZY]*grad_trial[Y]+dstress[sZY*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[Y] +
2206 (dstress[sZZ*nSymTen+sZX]*grad_trial[X]+dstress[sZZ*nSymTen+sZY]*grad_trial[Y]+dstress[sZZ*nSymTen+sZZ]*grad_trial[Z])*grad_test_dV[Z];
2207 }
2208 double ExteriorElementBoundaryStressFlux(const double& stressFlux,const double& disp_test_dS)
2209 {
2210 return stressFlux*disp_test_dS;
2211 }
2212 double ExteriorElementBoundaryStressFluxJacobian(const double& dstressFlux,const double& disp_test_dS)
2213 {
2214 return dstressFlux*disp_test_dS;
2215 }
2216};
2217
2218//specialization for 2D
2219template<int NDOF_MESH_TRIAL_ELEMENT, int NDOF_TRIAL_ELEMENT, int NDOF_TEST_ELEMENT>
2220class CompKernel<2,NDOF_MESH_TRIAL_ELEMENT,NDOF_TRIAL_ELEMENT,NDOF_TEST_ELEMENT>
2221{
2222public:
2224 const int X,
2235 X(mapping.X),
2236 Y(mapping.Y),
2242 XHX(mapping.XHX),
2243 YHX(mapping.YHX),
2245 {}
2246 inline void calculateG(double* jacInv,double* G,double& G_dd_G, double& tr_G)
2247 {
2248 //get the metric tensor
2249 //cek todo use symmetry
2250 for (int I=0;I<2;I++)
2251 for (int J=0;J<2;J++)
2252 {
2253 G[I*2+J] = 0.0;
2254 for (int K=0;K<2;K++)
2255 G[I*2+J] += jacInv[K*2+I]*jacInv[K*2+J];
2256 }
2257 G_dd_G = 0.0;
2258 tr_G = 0.0;
2259 for (int I=0;I<2;I++)
2260 {
2261 tr_G += G[I*2+I];
2262 for (int J=0;J<2;J++)
2263 {
2264 G_dd_G += G[I*2+J]*G[I*2+J];
2265 }
2266 }
2267 }
2268 inline void calculateGScale(double* G,double* v,double& h)
2269 {
2270 h = 0.0;
2271 for (int I=0;I<2;I++)
2272 for (int J=0;J<2;J++)
2273 h += v[I]*G[I*2+J]*v[J];
2274 h = 1.0/sqrt(h+1.0e-12);
2275 }
2276 inline void valFromDOF(const double* dof,const int* l2g_element,const double* trial_ref,double& val)
2277 {
2278 val=0.0;
2279 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2280 val+=dof[l2g_element[j]]*trial_ref[j];
2281 }
2282
2283 inline void gradFromDOF(const double* dof,const int* l2g_element,const double* grad_trial,double* grad)
2284 {
2285 for(int I=0;I<2;I++)
2286 grad[I] = 0.0;
2287 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2288 for(int I=0;I<2;I++)
2289 grad[I] += dof[l2g_element[j]]*grad_trial[j*2+I];
2290 }
2291
2292 inline void hessFromDOF(const double* dof,const int* l2g_element,const double* hess_trial,double* hess)
2293 {
2294 for(int I=0;I<2;I++)
2295 for(int J=0;J<2;J++)
2296 hess[I*2+J] = 0.0;
2297 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2298 for(int I=0;I<2;I++)
2299 for(int J=0;J<2;J++)
2300 hess[I*2+J] += dof[l2g_element[j]]*hess_trial[j*4+I*2+J];
2301 }
2302
2303 inline void valFromElementDOF(const double* dof,const double* trial_ref,double& val)
2304 {
2305 val=0.0;
2306 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2307 val+=dof[j]*trial_ref[j];
2308 }
2309
2310 inline void gradFromElementDOF(const double* dof,const double* grad_trial,double* grad)
2311 {
2312 for(int I=0;I<2;I++)
2313 grad[I] = 0.0;
2314 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2315 for(int I=0;I<2;I++)
2316 grad[I] += dof[j]*grad_trial[j*2+I];
2317 }
2318
2319 inline void gradTrialFromRef(const double* grad_trial_ref, const double* jacInv, double* grad_trial)
2320 {
2321 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2322 for(int I=0;I<2;I++)
2323 grad_trial[j*2+I] = 0.0;
2324 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2325 for(int I=0;I<2;I++)
2326 for(int J=0;J<2;J++)
2327 grad_trial[j*2+I] += jacInv[J*2+I]*grad_trial_ref[j*2+J];
2328 }
2329
2330 /*
2331 * DOFaverage
2332 * ----------
2333 *
2334 * Calculate the average DOF value for at a given mesh element.
2335 *
2336 * @param dof array of finite element DOF values
2337 * @param l2g_element local 2 global mapping for the current mesh element
2338 * @param val return value with the average DOF values
2339 */
2340
2341 inline void DOFaverage (const double* dof, const int* l2g_element, double& val)
2342 {
2343 val = 0.0;
2344
2345 for (int j=0; j<NDOF_MESH_TRIAL_ELEMENT; j++)
2346 val+=dof[l2g_element[j]];
2347
2348 val /= NDOF_MESH_TRIAL_ELEMENT;
2349 }
2350
2351
2352 inline void hessTrialFromRef(const double* hess_trial_ref, const double* jacInv, double* hess_trial)
2353 {
2354 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2355 for(int I=0;I<2;I++)
2356 for(int J=0;J<2;J++)
2357 hess_trial[j*4+I*2+J] = 0.0;
2358 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2359 for(int I=0;I<2;I++)
2360 for(int J=0;J<2;J++)
2361 for(int K=0;K<2;K++)
2362 for(int L=0;L<2;L++)
2363 hess_trial[j*4+I*2+J] += hess_trial_ref[j*4+K*2+L]*jacInv[L*2+J]*jacInv[K*2+I];
2364 }
2365
2366 inline void gradTestFromRef(const double* grad_test_ref, const double* jacInv, double* grad_test)
2367 {
2368 for (int i=0;i<NDOF_TEST_ELEMENT;i++)
2369 for(int I=0;I<2;I++)
2370 grad_test[i*2+I] = 0.0;
2371 for (int i=0;i<NDOF_TEST_ELEMENT;i++)
2372 for(int I=0;I<2;I++)
2373 for(int J=0;J<2;J++)
2374 grad_test[i*2+I] += jacInv[J*2+I]*grad_test_ref[i*2+J];
2375 }
2376
2377 inline void backwardEuler(const double& dt, const double& m_old, const double& m, const double& dm, double& mt, double& dmt)
2378 {
2379 mt =(m-m_old)/dt;
2380 dmt = dm/dt;
2381 }
2382
2383 inline void bdf(const double& alpha, const double& beta, const double& m, const double& dm, double& mt, double& dmt)
2384 {
2385 mt =alpha*m + beta;
2386 dmt = alpha*dm;
2387 }
2388
2389 inline void bdfC2(const double& alpha, const double& beta, const double& m, const double& dm, const double& dm2, double& mt, double& dmt, double& dm2t)
2390 {
2391 mt =alpha*m + beta;
2392 dmt = alpha*dm;
2393 dm2t = alpha*dm2;
2394 }
2395
2396 inline double Mass_weak(const double& mt, const double& w_dV)
2397 {
2398 return mt*w_dV;
2399 }
2400
2401 inline double MassJacobian_weak(const double& dmt,
2402 const double& v,
2403 const double& w_dV)
2404 {
2405 return dmt*v*w_dV;
2406 }
2407
2408 inline double Mass_strong(const double& mt)
2409 {
2410 return mt;
2411 }
2412
2413 inline double MassJacobian_strong(const double& dmt,
2414 const double& v)
2415 {
2416 return dmt*v;
2417 }
2418
2419 inline double Mass_adjoint(const double& dmt,
2420 const double& w_dV)
2421 {
2422 return dmt*w_dV;
2423 }
2424
2425 /*
2426 * pressureProjection_weak
2427 * -----------------------
2428 *
2429 * Inner product calculation of the pressure projection
2430 * stablization method of Bochev, Dohrmann and
2431 * Gunzburger (2006).
2432 *
2433 * @param viscosity viscosity at point
2434 * @param p the pressure value (either trial function
2435 * or actual value)
2436 * @param p_avg the pressure projection value (either
2437 * 1./3. for test functions of the average
2438 * value of the pressure on the element)
2439 * @param dV this is the integral weight
2440 *
2441 */
2442
2443 inline double pressureProjection_weak(const double& viscosity,
2444 const double& p,
2445 const double& p_avg,
2446 const double& q,
2447 const double& dV)
2448 {
2449 if (viscosity==0.){ return 0.;}
2450 return (1./viscosity)*(p-p_avg)*(q-1./3.)*dV;
2451 }
2452
2453 inline double Advection_weak(const double f[2],
2454 const double grad_w_dV[2])
2455 {
2456 double tmp=0.0;
2457 for(int I=0;I<2;I++)
2458 tmp -= f[I]*grad_w_dV[I];
2459 return tmp;
2460 }
2461
2462 inline double AdvectionJacobian_weak(const double df[2],
2463 const double& v,
2464 const double grad_w_dV[2])
2465 {
2466 double tmp=0.0;
2467 for(int I=0;I<2;I++)
2468 tmp -= df[I]*v*grad_w_dV[I];
2469 return tmp;
2470 }
2471
2472 inline double Advection_strong(const double df[2],
2473 const double grad_u[2])
2474 {
2475 double tmp=0.0;
2476 for(int I=0;I<2;I++)
2477 tmp += df[I]*grad_u[I];
2478 return tmp;
2479 }
2480
2481 inline double AdvectionJacobian_strong(const double df[2],
2482 const double grad_v[2])
2483 {
2484 double tmp=0.0;
2485 for(int I=0;I<2;I++)
2486 tmp += df[I]*grad_v[I];
2487 return tmp;
2488 }
2489
2490 inline double Advection_adjoint(const double df[2],
2491 const double grad_w_dV[2])
2492 {
2493 double tmp=0.0;
2494 for(int I=0;I<2;I++)
2495 tmp -= df[I]*grad_w_dV[I];
2496 return tmp;
2497 }
2498
2499 inline double Hamiltonian_weak(const double& H,
2500 const double& w_dV)
2501 {
2502 return H*w_dV;
2503 }
2504
2505 inline double HamiltonianJacobian_weak(const double dH[2],
2506 const double grad_v[2],
2507 const double& w_dV)
2508 {
2509 double tmp=0.0;
2510 for(int I=0;I<2;I++)
2511 tmp += dH[I]*grad_v[I]*w_dV;
2512 return tmp;
2513 }
2514
2515 inline double Hamiltonian_strong(const double dH[2],
2516 const double grad_u[2])
2517 {
2518 double tmp=0.0;
2519 for(int I=0;I<2;I++)
2520 tmp += dH[I]*grad_u[I];
2521 return tmp;
2522 }
2523
2524 inline double HamiltonianJacobian_strong(const double dH[2],
2525 const double grad_v[2])
2526 {
2527 double tmp=0.0;
2528 for(int I=0;I<2;I++)
2529 tmp += dH[I]*grad_v[I];
2530 return tmp;
2531 }
2532
2533 inline double Hamiltonian_adjoint(const double dH[2],
2534 const double grad_w_dV[2])
2535 {
2536 double tmp=0.0;
2537 for(int I=0;I<2;I++)
2538 tmp -= dH[I]*grad_w_dV[I];
2539 return tmp;
2540 }
2541
2542 inline double Diffusion_weak(int* rowptr,
2543 int* colind,
2544 double* a,
2545 const double grad_phi[2],
2546 const double grad_w_dV[2])
2547 {
2548 double tmp=0.0;
2549 for(int I=0;I<2;I++)
2550 for (int m=rowptr[I];m<rowptr[I+1];m++)
2551 tmp += a[m]*grad_phi[colind[m]]*grad_w_dV[I];
2552 return tmp;
2553 }
2554
2555 inline double DiffusionJacobian_weak(int* rowptr,
2556 int* colind,
2557 double* a,
2558 double* da,
2559 const double grad_phi[2],
2560 const double grad_w_dV[2],
2561 const double& dphi,
2562 const double& v,
2563 const double grad_v[2])
2564 {
2565 double daProduct=0.0,dphiProduct=0.0;
2566 for (int I=0;I<2;I++)
2567 for (int m=rowptr[I];m<rowptr[I+1];m++)
2568 {
2569 daProduct += da[m]*grad_phi[colind[m]]*grad_w_dV[I];
2570 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
2571 }
2572 return daProduct*v+dphiProduct*dphi;
2573 }
2574
2575 inline double SimpleDiffusionJacobian_weak(int* rowptr,
2576 int* colind,
2577 double* a,
2578 const double grad_v[2],
2579 const double grad_w_dV[2])
2580 {
2581 double dphiProduct=0.0;
2582 for (int I=0;I<2;I++)
2583 for (int m=rowptr[I];m<rowptr[I+1];m++)
2584 {
2585 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
2586 }
2587 return dphiProduct;
2588 }
2589
2590 inline double Reaction_weak(const double& r,
2591 const double& w_dV)
2592 {
2593 return r*w_dV;
2594 }
2595
2596 inline double ReactionJacobian_weak(const double& dr,
2597 const double& v,
2598 const double& w_dV)
2599 {
2600 return dr*v*w_dV;
2601 }
2602
2603 inline double Reaction_strong(const double& r)
2604 {
2605 return r;
2606 }
2607
2608 inline double ReactionJacobian_strong(const double& dr,
2609 const double& v)
2610 {
2611 return dr*v;
2612 }
2613
2614 inline double Reaction_adjoint(const double& dr,
2615 const double& w_dV)
2616 {
2617 return dr*w_dV;
2618 }
2619
2620 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
2621 const double& elementDiameter,
2622 const double& strong_residual,
2623 const double grad_u[2],
2624 double& numDiff)
2625 {
2626 double h,
2627 num,
2628 den,
2629 n_grad_u;
2630 h = elementDiameter;
2631 n_grad_u = 0.0;
2632 for (int I=0;I<2;I++)
2633 n_grad_u += grad_u[I]*grad_u[I];
2634 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
2635 den = sqrt(n_grad_u+1.0e-12);
2636 //cek hack shockCapturingDiffusion*fabs(strong_residual)*grad_phi_G_grad_phi
2637 numDiff = num/den;
2638 }
2639
2640
2641 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
2642 const double G[2*2],
2643 const double& strong_residual,
2644 const double grad_u[2],
2645 double& numDiff)
2646 {
2647 double den = 0.0;
2648 for (int I=0;I<2;I++)
2649 for (int J=0;J<2;J++)
2650 den += grad_u[I]*G[I*2+J]*grad_u[J];
2651
2652 numDiff = shockCapturingDiffusion*fabs(strong_residual)/(sqrt(den+1.0e-12));
2653 }
2654
2655 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
2656 const double& uref, const double& beta,
2657 const double G[2*2],
2658 const double& G_dd_G,
2659 const double& strong_residual,
2660 const double grad_u[2],
2661 double& numDiff)
2662 {
2663 double den = 0.0;
2664 for (int I=0;I<2;I++)
2665 for (int J=0;J<2;J++)
2666 den += grad_u[I]*G[I*2+J]*grad_u[J];
2667
2668 double h2_uref_1 = 1.0/(sqrt(den+1.0e-12));
2669 double h2_uref_2 = 1.0/(uref*sqrt(G_dd_G+1.0e-12));
2670 numDiff = shockCapturingDiffusion*fabs(strong_residual)*pow(h2_uref_1, 2.0-beta)*pow(h2_uref_2,beta-1.0);
2671 }
2672
2673
2674
2675 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
2676 const double G[2*2],
2677 const double& strong_residual,
2678 const double vel[2],
2679 const double grad_u[2],
2680 double& numDiff)
2681 {
2682 double den1 = 0.0,den2=0.0, nom=0.0;
2683 for (int I=0;I<2;I++)
2684 {
2685 nom += vel[I]*vel[I];
2686 den2+= grad_u[I]*grad_u[I];
2687 for (int J=0;J<2;J++)
2688 den1 += vel[I]*G[I*2+J]*vel[J];
2689 }
2690 numDiff = shockCapturingDiffusion*fabs(strong_residual)*(sqrt(nom/(den1*den2 + 1.0e-12)));
2691 }
2692
2693
2694 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
2695 const double& elementDiameter,
2696 const double& strong_residual,
2697 const double grad_u[2],
2698 double& gradNorm,
2699 double& gradNorm_last,
2700 double& numDiff)
2701 {
2702 double h,
2703 num,
2704 n_grad_u;
2705 h = elementDiameter;
2706 n_grad_u = 0.0;
2707 for (int I=0;I<2;I++)
2708 n_grad_u += grad_u[I]*grad_u[I];
2709 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
2710 gradNorm = sqrt(n_grad_u+1.0e-12);
2711 //cek hack shockCapturingDiffusion*fabs(strong_residual)*grad_phi_G_grad_phi
2712 numDiff = num/gradNorm_last;
2713 }
2714
2715 inline double SubgridError(const double& error,
2716 const double& Lstar_w_dV)
2717 {
2718 return error*Lstar_w_dV;
2719 }
2720
2721 inline double SubgridErrorJacobian(const double& derror,
2722 const double& Lstar_w_dV)
2723 {
2724 return derror*Lstar_w_dV;
2725 }
2726
2727 inline double NumericalDiffusion(const double& numDiff,
2728 const double grad_u[2],
2729 const double grad_w_dV[2])
2730 {
2731 double tmp=0.0;
2732 for (int I=0;I<2;I++)
2733 tmp += numDiff*grad_u[I]*grad_w_dV[I];
2734 return tmp;
2735 }
2736
2737 inline double NumericalDiffusionJacobian(const double& numDiff,
2738 const double grad_v[2],
2739 const double grad_w_dV[2])
2740 {
2741 double tmp=0.0;
2742 for (int I=0;I<2;I++)
2743 tmp += numDiff*grad_v[I]*grad_w_dV[I];
2744 return tmp;
2745 }
2746
2747
2748
2749 inline double ExteriorElementBoundaryFlux(const double& flux,
2750 const double& w_dS)
2751 {
2752 return flux*w_dS;
2753 }
2754
2755 inline double InteriorElementBoundaryFlux(const double& flux,
2756 const double& w_dS)
2757 {
2758 return flux*w_dS;
2759 }
2760
2761 inline double ExteriorNumericalAdvectiveFluxJacobian(const double& dflux_left,
2762 const double& v)
2763 {
2764 return dflux_left*v;
2765 }
2766
2767 inline double InteriorNumericalAdvectiveFluxJacobian(const double& dflux_left,
2768 const double& v)
2769 {
2770 return dflux_left*v;
2771 }
2772
2773 inline double ExteriorElementBoundaryScalarDiffusionAdjoint(const int& isDOFBoundary,
2774 const int& isFluxBoundary,
2775 const double& sigma,
2776 const double& u,
2777 const double& bc_u,
2778 const double normal[2],
2779 const double& a,
2780 const double grad_w_dS[2])
2781 {
2782 double tmp=0.0;
2783 for(int I=0;I<2;I++)
2784 {
2785 tmp += normal[I]*grad_w_dS[I];
2786 }
2787 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*(u-bc_u)*a;
2788 return tmp;
2789 }
2790
2791 inline double ExteriorElementBoundaryScalarDiffusionAdjointJacobian(const int& isDOFBoundary,
2792 const int& isFluxBoundary,
2793 const double& sigma,
2794 const double& v,
2795 const double normal[2],
2796 const double& a,
2797 const double grad_w_dS[2])
2798 {
2799 double tmp=0.0;
2800 for(int I=0;I<2;I++)
2801 {
2802 tmp += normal[I]*grad_w_dS[I];
2803 }
2804 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*v*a;
2805 return tmp;
2806 }
2807
2808 inline double ExteriorElementBoundaryDiffusionAdjoint(const int& isDOFBoundary,
2809 const int& isFluxBoundary,
2810 const double& sigma,
2811 const double& u,
2812 const double& bc_u,
2813 const double normal[2],
2814 int* rowptr,
2815 int* colind,
2816 double* a,
2817 const double grad_w_dS[2])
2818 {
2819 double tmp=0.0;
2820 for(int I=0;I<2;I++)
2821 for (int m=rowptr[I];m<rowptr[I+1];m++)
2822 tmp += (1.0-isFluxBoundary)*isDOFBoundary*sigma*(u-bc_u)*a[m]*normal[colind[m]]*grad_w_dS[I];
2823 return tmp;
2824 }
2825
2826 inline double ExteriorElementBoundaryDiffusionAdjointJacobian(const int& isDOFBoundary,
2827 const int& isFluxBoundary,
2828 const double& sigma,
2829 const double& v,
2830 const double normal[2],
2831 int* rowptr,
2832 int* colind,
2833 double* a,
2834 const double grad_w_dS[2])
2835 {
2836 double tmp=0.0;
2837 for(int I=0;I<2;I++)
2838 for (int m=rowptr[I];m<rowptr[I+1];m++)
2839 tmp += (1.0-isFluxBoundary)*isDOFBoundary*sigma*v*a[m]*normal[colind[m]]*grad_w_dS[I];
2840 return tmp;
2841 }
2842
2843 inline void calculateMapping_element(const int eN,
2844 const int k,
2845 double* mesh_dof,
2846 int* mesh_l2g,
2847 //xt::pyarray<double>& mesh_trial_ref,
2848 double* mesh_trial_ref,
2849 double* mesh_grad_trial_ref,
2850 double* jac,
2851 double& jacDet,
2852 double* jacInv,
2853 double& x,
2854 double& y)
2855 {
2856 mapping.calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y);
2857 }
2858
2859 inline void calculateH_element(const int eN,
2860 const int k,
2861 double* h_dof,
2862 int* mesh_l2g,
2863 //xt::pyarray<double>& mesh_trial_ref,
2864 double* mesh_trial_ref,
2865 double& h)
2866 {
2868 k,
2869 h_dof,
2870 mesh_l2g,
2871 mesh_trial_ref,
2872 h);
2873 }
2874
2875 inline void calculateMapping_element(const int eN,
2876 const int k,
2877 double* mesh_dof,
2878 int* mesh_l2g,
2879 //xt::pyarray<double>& mesh_trial_ref,
2880 double* mesh_trial_ref,
2881 double* mesh_grad_trial_ref,
2882 double* jac,
2883 double& jacDet,
2884 double* jacInv,
2885 double& x,
2886 double& y,
2887 double& z)
2888 {
2889 mapping.calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y,z);
2890 }
2891
2892 inline void calculateMappingVelocity_element(const int eN,
2893 const int k,
2894 double* meshVelocity_dof,
2895 int* mesh_l2g,
2896 //xt::pyarray<double>& mesh_trial_ref,
2897 double* mesh_trial_ref,
2898 double& xt,
2899 double& yt)
2900 {
2901 mapping.calculateMappingVelocity_element(eN,k,meshVelocity_dof,mesh_l2g,mesh_trial_ref,xt,yt);
2902 }
2903
2904 inline void calculateMappingVelocity_element(const int eN,
2905 const int k,
2906 double* meshVelocity_dof,
2907 int* mesh_l2g,
2908 //xt::pyarray<double>& mesh_trial_ref,
2909 double* mesh_trial_ref,
2910 double& xt,
2911 double& yt,
2912 double& zt)
2913 {
2914 mapping.calculateMappingVelocity_element(eN,k,meshVelocity_dof,mesh_l2g,mesh_trial_ref,xt,yt,zt);
2915 }
2916
2917 inline
2919 const int ebN_local,
2920 const int kb,
2921 const int ebN_local_kb,
2922 double* mesh_dof,
2923 int* mesh_l2g,
2924 double* mesh_trial_trace_ref,
2925 double* mesh_grad_trial_trace_ref,
2926 double* boundaryJac_ref,
2927 double* jac,
2928 double& jacDet,
2929 double* jacInv,
2930 double* boundaryJac,
2931 double* metricTensor,
2932 double& metricTensorDetSqrt,
2933 double* normal_ref,
2934 double* normal,
2935 double& x,
2936 double& y)
2937 {
2939 ebN_local,
2940 kb,
2941 ebN_local_kb,
2942 mesh_dof,
2943 mesh_l2g,
2944 mesh_trial_trace_ref,
2945 mesh_grad_trial_trace_ref,
2946 boundaryJac_ref,
2947 jac,
2948 jacDet,
2949 jacInv,
2950 boundaryJac,
2951 metricTensor,
2952 metricTensorDetSqrt,
2953 normal_ref,
2954 normal,
2955 x,
2956 y);
2957 }
2958
2959 inline
2961 const int ebN_local,
2962 const int kb,
2963 const int ebN_local_kb,
2964 double* mesh_dof,
2965 int* mesh_l2g,
2966 double* mesh_trial_trace_ref,
2967 double* mesh_grad_trial_trace_ref,
2968 double* boundaryJac_ref,
2969 double* jac,
2970 double& jacDet,
2971 double* jacInv,
2972 double* boundaryJac,
2973 double* metricTensor,
2974 double& metricTensorDetSqrt,
2975 double* normal_ref,
2976 double* normal,
2977 double& x,
2978 double& y,
2979 double& z)
2980 {
2982 ebN_local,
2983 kb,
2984 ebN_local_kb,
2985 mesh_dof,
2986 mesh_l2g,
2987 mesh_trial_trace_ref,
2988 mesh_grad_trial_trace_ref,
2989 boundaryJac_ref,
2990 jac,
2991 jacDet,
2992 jacInv,
2993 boundaryJac,
2994 metricTensor,
2995 metricTensorDetSqrt,
2996 normal_ref,
2997 normal,
2998 x,
2999 y,
3000 z);
3001 }
3002
3003 inline
3005 const int ebN_local,
3006 const int kb,
3007 const int ebN_local_kb,
3008 double* mesh_velocity_dof,
3009 int* mesh_l2g,
3010 double* mesh_trial_trace_ref,
3011 double& xt,
3012 double& yt,
3013 double* normal,
3014 double* boundaryJac,
3015 double* metricTensor,
3016 double& metricTensorDetSqrt)
3017 {
3019 ebN_local,
3020 kb,
3021 ebN_local_kb,
3022 mesh_velocity_dof,
3023 mesh_l2g,
3024 mesh_trial_trace_ref,
3025 xt,
3026 yt,
3027 normal,
3028 boundaryJac,
3029 metricTensor,
3030 metricTensorDetSqrt);
3031 }
3032 inline
3034 const int ebN_local,
3035 const int kb,
3036 const int ebN_local_kb,
3037 double* mesh_velocity_dof,
3038 int* mesh_l2g,
3039 double* mesh_trial_trace_ref,
3040 double& xt,
3041 double& yt,
3042 double& zt,
3043 double* normal,
3044 double* boundaryJac,
3045 double* metricTensor,
3046 double& metricTensorDetSqrt)
3047 {
3049 ebN_local,
3050 kb,
3051 ebN_local_kb,
3052 mesh_velocity_dof,
3053 mesh_l2g,
3054 mesh_trial_trace_ref,
3055 xt,
3056 yt,
3057 zt,
3058 normal,
3059 boundaryJac,
3060 metricTensor,
3061 metricTensorDetSqrt);
3062 }
3063 double Stress_u_weak(double* stress, double* grad_test_dV)
3064 {
3065 return stress[sXX]*grad_test_dV[X] + stress[sXY]*grad_test_dV[Y];
3066 }
3067 double StressJacobian_u_u_weak(double* dstress, double* grad_trial, double* grad_test_dV)
3068 {
3069 return
3070 (dstress[sXX*nSymTen+sXX]*grad_trial[X]+dstress[sXX*nSymTen+sXY]*grad_trial[Y])*grad_test_dV[X] +
3071 (dstress[sXY*nSymTen+sXX]*grad_trial[X]+dstress[sXY*nSymTen+sXY]*grad_trial[Y])*grad_test_dV[Y];
3072 }
3073 double StressJacobian_u_v_weak(double* dstress, double* grad_trial, double* grad_test_dV)
3074 {
3075 return
3076 (dstress[sXX*nSymTen+sYX]*grad_trial[X]+dstress[sXX*nSymTen+sYY]*grad_trial[Y])*grad_test_dV[X] +
3077 (dstress[sXY*nSymTen+sYX]*grad_trial[X]+dstress[sXY*nSymTen+sYY]*grad_trial[Y])*grad_test_dV[Y];
3078 }
3079 double Stress_v_weak(double* stress, double* grad_test_dV)
3080 {
3081 return stress[sYX]*grad_test_dV[X] + stress[sYY]*grad_test_dV[Y];
3082 }
3083 double StressJacobian_v_u_weak(double* dstress, double* grad_trial,double* grad_test_dV)
3084 {
3085 return
3086 (dstress[sYX*nSymTen+sXX]*grad_trial[X]+dstress[sYX*nSymTen+sXY]*grad_trial[Y])*grad_test_dV[X] +
3087 (dstress[sYY*nSymTen+sXX]*grad_trial[X]+dstress[sYY*nSymTen+sXY]*grad_trial[Y])*grad_test_dV[Y];
3088 }
3089 double StressJacobian_v_v_weak(double* dstress, double* grad_trial,double* grad_test_dV)
3090 {
3091 return
3092 (dstress[sYX*nSymTen+sYX]*grad_trial[X]+dstress[sYX*nSymTen+sYY]*grad_trial[Y])*grad_test_dV[X] +
3093 (dstress[sYY*nSymTen+sYX]*grad_trial[X]+dstress[sYY*nSymTen+sYY]*grad_trial[Y])*grad_test_dV[Y];
3094 }
3095 double ExteriorElementBoundaryStressFlux(const double& stressFlux,const double& disp_test_dS)
3096 {
3097 return stressFlux*disp_test_dS;
3098 }
3099 double ExteriorElementBoundaryStressFluxJacobian(const double& dstressFlux,const double& disp_test_dS)
3100 {
3101 return dstressFlux*disp_test_dS;
3102 }
3103};
3104
3105//specialization for 1D
3106template<int NDOF_MESH_TRIAL_ELEMENT, int NDOF_TRIAL_ELEMENT, int NDOF_TEST_ELEMENT>
3107class CompKernel<1,NDOF_MESH_TRIAL_ELEMENT,NDOF_TRIAL_ELEMENT,NDOF_TEST_ELEMENT>
3108{
3109public:
3111 const int X,
3118 X(mapping.X),
3119 XX(mapping.XX),
3120 sXX(mapping.sXX),
3122 XHX(mapping.XHX),
3124 {}
3125 inline void calculateG(double* jacInv,double* G,double& G_dd_G, double& tr_G)
3126 {
3127 for (int I=0;I<1;I++)
3128 for (int J=0;J<1;J++)
3129 {
3130 G[I+J] = 0.0;
3131 for (int K=0;K<1;K++)
3132 G[I+J] += jacInv[K+I]*jacInv[K+J];
3133 }
3134 G_dd_G = 0.0;
3135 tr_G = 0.0;
3136 for (int I=0;I<1;I++)
3137 {
3138 tr_G += G[I+I];
3139 for (int J=0;J<1;J++)
3140 {
3141 G_dd_G += G[I+J]*G[I+J];
3142 }
3143 }
3144 }
3145 inline void calculateGScale(double* G,double* v,double& h)
3146 {
3147 h = 0.0;
3148 for (int I=0;I<1;I++)
3149 for (int J=0;J<1;J++)
3150 h += v[I]*G[I+J]*v[J];
3151 h = 1.0/sqrt(h+1.0e-12);
3152 }
3153 inline void valFromDOF(const double* dof,const int* l2g_element,const double* trial_ref,double& val)
3154 {
3155 val=0.0;
3156 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3157 val+=dof[l2g_element[j]]*trial_ref[j];
3158 }
3159
3160 inline void gradFromDOF(const double* dof,const int* l2g_element,const double* grad_trial,double* grad)
3161 {
3162 for(int I=0;I<1;I++)
3163 grad[I] = 0.0;
3164 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3165 for(int I=0;I<1;I++)
3166 grad[I] += dof[l2g_element[j]]*grad_trial[j+I];
3167 }
3168
3169 inline void hessFromDOF(const double* dof,const int* l2g_element,const double* hess_trial,double* hess)
3170 {
3171 for(int I=0;I<1;I++)
3172 for(int J=0;J<1;J++)
3173 hess[I+J] = 0.0;
3174 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3175 for(int I=0;I<1;I++)
3176 for(int J=0;J<1;J++)
3177 hess[I*2+J] += dof[l2g_element[j]]*hess_trial[j+I+J];
3178 }
3179
3180 inline void valFromElementDOF(const double* dof,const double* trial_ref,double& val)
3181 {
3182 val=0.0;
3183 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3184 val+=dof[j]*trial_ref[j];
3185 }
3186
3187 inline void gradFromElementDOF(const double* dof,const double* grad_trial,double* grad)
3188 {
3189 for(int I=0;I<1;I++)
3190 grad[I] = 0.0;
3191 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3192 for(int I=0;I<1;I++)
3193 grad[I] += dof[j]*grad_trial[j+I];
3194 }
3195
3196 inline void gradTrialFromRef(const double* grad_trial_ref, const double* jacInv, double* grad_trial)
3197 {
3198 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3199 for(int I=0;I<1;I++)
3200 grad_trial[j+I] = 0.0;
3201 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3202 for(int I=0;I<1;I++)
3203 for(int J=0;J<1;J++)
3204 grad_trial[j+I] += jacInv[J+I]*grad_trial_ref[j+J];
3205 }
3206
3207 /*
3208 * DOFaverage
3209 * ----------
3210 *
3211 * Calculate the average DOF value for at a given mesh element.
3212 *
3213 * @param dof array of finite element DOF values
3214 * @param l2g_element local 2 global mapping for the current mesh element
3215 * @param val return value with the average DOF values
3216 */
3217
3218 inline void DOFaverage (const double* dof, const int* l2g_element, double& val)
3219 {
3220 val = 0.0;
3221
3222 for (int j=0; j<NDOF_MESH_TRIAL_ELEMENT; j++)
3223 val+=dof[l2g_element[j]];
3224
3225 val /= NDOF_MESH_TRIAL_ELEMENT;
3226 }
3227
3228
3229 inline void hessTrialFromRef(const double* hess_trial_ref, const double* jacInv, double* hess_trial)
3230 {
3231 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3232 for(int I=0;I<1;I++)
3233 for(int J=0;J<1;J++)
3234 hess_trial[j+I+J] = 0.0;
3235 for (int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3236 for(int I=0;I<1;I++)
3237 for(int J=0;J<1;J++)
3238 for(int K=0;K<1;K++)
3239 for(int L=0;L<1;L++)
3240 hess_trial[j+I+J] += hess_trial_ref[j+K+L]*jacInv[L+J]*jacInv[K+I];
3241 }
3242
3243 inline void gradTestFromRef(const double* grad_test_ref, const double* jacInv, double* grad_test)
3244 {
3245 for (int i=0;i<NDOF_TEST_ELEMENT;i++)
3246 for(int I=0;I<1;I++)
3247 grad_test[i+I] = 0.0;
3248 for (int i=0;i<NDOF_TEST_ELEMENT;i++)
3249 for(int I=0;I<1;I++)
3250 for(int J=0;J<1;J++)
3251 grad_test[i+I] += jacInv[J+I]*grad_test_ref[i+J];
3252 }
3253
3254 inline void backwardEuler(const double& dt, const double& m_old, const double& m, const double& dm, double& mt, double& dmt)
3255 {
3256 mt =(m-m_old)/dt;
3257 dmt = dm/dt;
3258 }
3259
3260 inline void bdf(const double& alpha, const double& beta, const double& m, const double& dm, double& mt, double& dmt)
3261 {
3262 mt =alpha*m + beta;
3263 dmt = alpha*dm;
3264 }
3265
3266 inline void bdfC2(const double& alpha, const double& beta, const double& m, const double& dm, const double& dm2, double& mt, double& dmt, double& dm2t)
3267 {
3268 mt =alpha*m + beta;
3269 dmt = alpha*dm;
3270 dm2t = alpha*dm2;
3271 }
3272
3273 inline double Mass_weak(const double& mt, const double& w_dV)
3274 {
3275 return mt*w_dV;
3276 }
3277
3278 inline double MassJacobian_weak(const double& dmt,
3279 const double& v,
3280 const double& w_dV)
3281 {
3282 return dmt*v*w_dV;
3283 }
3284
3285 inline double Mass_strong(const double& mt)
3286 {
3287 return mt;
3288 }
3289
3290 inline double MassJacobian_strong(const double& dmt,
3291 const double& v)
3292 {
3293 return dmt*v;
3294 }
3295
3296 inline double Mass_adjoint(const double& dmt,
3297 const double& w_dV)
3298 {
3299 return dmt*w_dV;
3300 }
3301
3302 /*
3303 * pressureProjection_weak
3304 * -----------------------
3305 *
3306 * Inner product calculation of the pressure projection
3307 * stablization method of Bochev, Dohrmann and
3308 * Gunzburger (2006).
3309 *
3310 * @param viscosity viscosity at point
3311 * @param p the pressure value (either trial function
3312 * or actual value)
3313 * @param p_avg the pressure projection value (either
3314 * 1./3. for test functions of the average
3315 * value of the pressure on the element)
3316 * @param dV this is the integral weight
3317 *
3318 */
3319
3320 inline double pressureProjection_weak(const double& viscosity,
3321 const double& p,
3322 const double& p_avg,
3323 const double& q,
3324 const double& dV)
3325 {
3326 if (viscosity==0.){ return 0.;}
3327 return (1./viscosity)*(p-p_avg)*(q-1./3.)*dV;
3328 }
3329
3330 inline double Advection_weak(const double f[1],
3331 const double grad_w_dV[1])
3332 {
3333 double tmp=0.0;
3334 for(int I=0;I<1;I++)
3335 tmp -= f[I]*grad_w_dV[I];
3336 return tmp;
3337 }
3338
3339 inline double AdvectionJacobian_weak(const double df[1],
3340 const double& v,
3341 const double grad_w_dV[1])
3342 {
3343 double tmp=0.0;
3344 for(int I=0;I<1;I++)
3345 tmp -= df[I]*v*grad_w_dV[I];
3346 return tmp;
3347 }
3348
3349 inline double Advection_strong(const double df[1],
3350 const double grad_u[1])
3351 {
3352 double tmp=0.0;
3353 for(int I=0;I<1;I++)
3354 tmp += df[I]*grad_u[I];
3355 return tmp;
3356 }
3357
3358 inline double AdvectionJacobian_strong(const double df[1],
3359 const double grad_v[1])
3360 {
3361 double tmp=0.0;
3362 for(int I=0;I<1;I++)
3363 tmp += df[I]*grad_v[I];
3364 return tmp;
3365 }
3366
3367 inline double Advection_adjoint(const double df[1],
3368 const double grad_w_dV[1])
3369 {
3370 double tmp=0.0;
3371 for(int I=0;I<1;I++)
3372 tmp -= df[I]*grad_w_dV[I];
3373 return tmp;
3374 }
3375
3376 inline double Hamiltonian_weak(const double& H,
3377 const double& w_dV)
3378 {
3379 return H*w_dV;
3380 }
3381
3382 inline double HamiltonianJacobian_weak(const double dH[1],
3383 const double grad_v[1],
3384 const double& w_dV)
3385 {
3386 double tmp=0.0;
3387 for(int I=0;I<1;I++)
3388 tmp += dH[I]*grad_v[I]*w_dV;
3389 return tmp;
3390 }
3391
3392 inline double Hamiltonian_strong(const double dH[1],
3393 const double grad_u[1])
3394 {
3395 double tmp=0.0;
3396 for(int I=0;I<1;I++)
3397 tmp += dH[I]*grad_u[I];
3398 return tmp;
3399 }
3400
3401 inline double HamiltonianJacobian_strong(const double dH[1],
3402 const double grad_v[1])
3403 {
3404 double tmp=0.0;
3405 for(int I=0;I<1;I++)
3406 tmp += dH[I]*grad_v[I];
3407 return tmp;
3408 }
3409
3410 inline double Hamiltonian_adjoint(const double dH[1],
3411 const double grad_w_dV[1])
3412 {
3413 double tmp=0.0;
3414 for(int I=0;I<1;I++)
3415 tmp -= dH[I]*grad_w_dV[I];
3416 return tmp;
3417 }
3418
3419 inline double Diffusion_weak(int* rowptr,
3420 int* colind,
3421 double* a,
3422 const double grad_phi[1],
3423 const double grad_w_dV[1])
3424 {
3425 double tmp=0.0;
3426 for(int I=0;I<1;I++)
3427 for (int m=rowptr[I];m<rowptr[I+1];m++)
3428 tmp += a[m]*grad_phi[colind[m]]*grad_w_dV[I];
3429 return tmp;
3430 }
3431
3432 inline double DiffusionJacobian_weak(int* rowptr,
3433 int* colind,
3434 double* a,
3435 double* da,
3436 const double grad_phi[1],
3437 const double grad_w_dV[1],
3438 const double& dphi,
3439 const double& v,
3440 const double grad_v[1])
3441 {
3442 double daProduct=0.0,dphiProduct=0.0;
3443 for (int I=0;I<1;I++)
3444 for (int m=rowptr[I];m<rowptr[I+1];m++)
3445 {
3446 daProduct += da[m]*grad_phi[colind[m]]*grad_w_dV[I];
3447 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
3448 }
3449 return daProduct*v+dphiProduct*dphi;
3450 }
3451
3452 inline double SimpleDiffusionJacobian_weak(int* rowptr,
3453 int* colind,
3454 double* a,
3455 const double grad_v[1],
3456 const double grad_w_dV[1])
3457 {
3458 double dphiProduct=0.0;
3459 for (int I=0;I<1;I++)
3460 for (int m=rowptr[I];m<rowptr[I+1];m++)
3461 {
3462 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
3463 }
3464 return dphiProduct;
3465 }
3466
3467 inline double Reaction_weak(const double& r,
3468 const double& w_dV)
3469 {
3470 return r*w_dV;
3471 }
3472
3473 inline double ReactionJacobian_weak(const double& dr,
3474 const double& v,
3475 const double& w_dV)
3476 {
3477 return dr*v*w_dV;
3478 }
3479
3480 inline double Reaction_strong(const double& r)
3481 {
3482 return r;
3483 }
3484
3485 inline double ReactionJacobian_strong(const double& dr,
3486 const double& v)
3487 {
3488 return dr*v;
3489 }
3490
3491 inline double Reaction_adjoint(const double& dr,
3492 const double& w_dV)
3493 {
3494 return dr*w_dV;
3495 }
3496
3497 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
3498 const double& elementDiameter,
3499 const double& strong_residual,
3500 const double grad_u[1],
3501 double& numDiff)
3502 {
3503 double h,
3504 num,
3505 den,
3506 n_grad_u;
3507 h = elementDiameter;
3508 n_grad_u = 0.0;
3509 for (int I=0;I<1;I++)
3510 n_grad_u += grad_u[I]*grad_u[I];
3511 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
3512 den = sqrt(n_grad_u+1.0e-12);
3513 //cek hack shockCapturingDiffusion*fabs(strong_residual)*grad_phi_G_grad_phi
3514 numDiff = num/den;
3515 }
3516
3517
3518 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
3519 const double G[1*1],
3520 const double& strong_residual,
3521 const double grad_u[1],
3522 double& numDiff)
3523 {
3524 double den = 0.0;
3525 for (int I=0;I<1;I++)
3526 for (int J=0;J<1;J++)
3527 den += grad_u[I]*G[I+J]*grad_u[J];
3528
3529 numDiff = shockCapturingDiffusion*fabs(strong_residual)/(sqrt(den+1.0e-12));
3530 }
3531
3532 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
3533 const double& uref, const double& beta,
3534 const double G[1*1],
3535 const double& G_dd_G,
3536 const double& strong_residual,
3537 const double grad_u[1],
3538 double& numDiff)
3539 {
3540 double den = 0.0;
3541 for (int I=0;I<1;I++)
3542 for (int J=0;J<1;J++)
3543 den += grad_u[I]*G[I+J]*grad_u[J];
3544
3545 double h2_uref_1 = 1.0/(sqrt(den+1.0e-12));
3546 double h2_uref_2 = 1.0/(uref*sqrt(G_dd_G+1.0e-12));
3547 numDiff = shockCapturingDiffusion*fabs(strong_residual)*pow(h2_uref_1, 2.0-beta)*pow(h2_uref_2,beta-1.0);
3548 }
3549
3550
3551
3552 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
3553 const double G[1*1],
3554 const double& strong_residual,
3555 const double vel[1],
3556 const double grad_u[1],
3557 double& numDiff)
3558 {
3559 double den1 = 0.0,den2=0.0, nom=0.0;
3560 for (int I=0;I<1;I++)
3561 {
3562 nom += vel[I]*vel[I];
3563 den2+= grad_u[I]*grad_u[I];
3564 for (int J=0;J<1;J++)
3565 den1 += vel[I]*G[I+J]*vel[J];
3566 }
3567 numDiff = shockCapturingDiffusion*fabs(strong_residual)*(sqrt(nom/(den1*den2 + 1.0e-12)));
3568 }
3569
3570
3571 inline void calculateNumericalDiffusion(const double& shockCapturingDiffusion,
3572 const double& elementDiameter,
3573 const double& strong_residual,
3574 const double grad_u[1],
3575 double& gradNorm,
3576 double& gradNorm_last,
3577 double& numDiff)
3578 {
3579 double h,
3580 num,
3581 n_grad_u;
3582 h = elementDiameter;
3583 n_grad_u = 0.0;
3584 for (int I=0;I<1;I++)
3585 n_grad_u += grad_u[I]*grad_u[I];
3586 num = shockCapturingDiffusion*0.5*h*fabs(strong_residual);
3587 gradNorm = sqrt(n_grad_u+1.0e-12);
3588 //cek hack shockCapturingDiffusion*fabs(strong_residual)*grad_phi_G_grad_phi
3589 numDiff = num/gradNorm_last;
3590 }
3591
3592 inline double SubgridError(const double& error,
3593 const double& Lstar_w_dV)
3594 {
3595 return error*Lstar_w_dV;
3596 }
3597
3598 inline double SubgridErrorJacobian(const double& derror,
3599 const double& Lstar_w_dV)
3600 {
3601 return derror*Lstar_w_dV;
3602 }
3603
3604 inline double NumericalDiffusion(const double& numDiff,
3605 const double grad_u[1],
3606 const double grad_w_dV[1])
3607 {
3608 double tmp=0.0;
3609 for (int I=0;I<1;I++)
3610 tmp += numDiff*grad_u[I]*grad_w_dV[I];
3611 return tmp;
3612 }
3613
3614 inline double NumericalDiffusionJacobian(const double& numDiff,
3615 const double grad_v[1],
3616 const double grad_w_dV[1])
3617 {
3618 double tmp=0.0;
3619 for (int I=0;I<1;I++)
3620 tmp += numDiff*grad_v[I]*grad_w_dV[I];
3621 return tmp;
3622 }
3623
3624
3625
3626 inline double ExteriorElementBoundaryFlux(const double& flux,
3627 const double& w_dS)
3628 {
3629 return flux*w_dS;
3630 }
3631
3632 inline double InteriorElementBoundaryFlux(const double& flux,
3633 const double& w_dS)
3634 {
3635 return flux*w_dS;
3636 }
3637
3638 inline double ExteriorNumericalAdvectiveFluxJacobian(const double& dflux_left,
3639 const double& v)
3640 {
3641 return dflux_left*v;
3642 }
3643
3644 inline double InteriorNumericalAdvectiveFluxJacobian(const double& dflux_left,
3645 const double& v)
3646 {
3647 return dflux_left*v;
3648 }
3649
3650 inline double ExteriorElementBoundaryScalarDiffusionAdjoint(const int& isDOFBoundary,
3651 const int& isFluxBoundary,
3652 const double& sigma,
3653 const double& u,
3654 const double& bc_u,
3655 const double normal[1],
3656 const double& a,
3657 const double grad_w_dS[1])
3658 {
3659 double tmp=0.0;
3660 for(int I=0;I<1;I++)
3661 {
3662 tmp += normal[I]*grad_w_dS[I];
3663 }
3664 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*(u-bc_u)*a;
3665 return tmp;
3666 }
3667
3668 inline double ExteriorElementBoundaryScalarDiffusionAdjointJacobian(const int& isDOFBoundary,
3669 const int& isFluxBoundary,
3670 const double& sigma,
3671 const double& v,
3672 const double normal[1],
3673 const double& a,
3674 const double grad_w_dS[1])
3675 {
3676 double tmp=0.0;
3677 for(int I=0;I<1;I++)
3678 {
3679 tmp += normal[I]*grad_w_dS[I];
3680 }
3681 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*v*a;
3682 return tmp;
3683 }
3684
3685 inline double ExteriorElementBoundaryDiffusionAdjoint(const int& isDOFBoundary,
3686 const int& isFluxBoundary,
3687 const double& sigma,
3688 const double& u,
3689 const double& bc_u,
3690 const double normal[1],
3691 int* rowptr,
3692 int* colind,
3693 double* a,
3694 const double grad_w_dS[1])
3695 {
3696 double tmp=0.0;
3697 for(int I=0;I<1;I++)
3698 for (int m=rowptr[I];m<rowptr[I+1];m++)
3699 tmp += (1.0-isFluxBoundary)*isDOFBoundary*sigma*(u-bc_u)*a[m]*normal[colind[m]]*grad_w_dS[I];
3700 return tmp;
3701 }
3702
3703 inline double ExteriorElementBoundaryDiffusionAdjointJacobian(const int& isDOFBoundary,
3704 const int& isFluxBoundary,
3705 const double& sigma,
3706 const double& v,
3707 const double normal[1],
3708 int* rowptr,
3709 int* colind,
3710 double* a,
3711 const double grad_w_dS[1])
3712 {
3713 double tmp=0.0;
3714 for(int I=0;I<1;I++)
3715 for (int m=rowptr[I];m<rowptr[I+1];m++)
3716 tmp += (1.0-isFluxBoundary)*isDOFBoundary*sigma*v*a[m]*normal[colind[m]]*grad_w_dS[I];
3717 return tmp;
3718 }
3719
3720 inline void calculateMapping_element(const int eN,
3721 const int k,
3722 double* mesh_dof,
3723 int* mesh_l2g,
3724 //xt::pyarray<double>& mesh_trial_ref,
3725 double* mesh_trial_ref,
3726 double* mesh_grad_trial_ref,
3727 double* jac,
3728 double& jacDet,
3729 double* jacInv,
3730 double& x,
3731 double& y)
3732 {
3733 mapping.calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y);
3734 }
3735
3736 inline void calculateH_element(const int eN,
3737 const int k,
3738 double* h_dof,
3739 int* mesh_l2g,
3740 //xt::pyarray<double>& mesh_trial_ref,
3741 double* mesh_trial_ref,
3742 double& h)
3743 {
3745 k,
3746 h_dof,
3747 mesh_l2g,
3748 mesh_trial_ref,
3749 h);
3750 }
3751
3752 inline void calculateMapping_element(const int eN,
3753 const int k,
3754 double* mesh_dof,
3755 int* mesh_l2g,
3756 //xt::pyarray<double>& mesh_trial_ref,
3757 double* mesh_trial_ref,
3758 double* mesh_grad_trial_ref,
3759 double* jac,
3760 double& jacDet,
3761 double* jacInv,
3762 double& x,
3763 double& y,
3764 double& z)
3765 {
3766 mapping.calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y,z);
3767 }
3768
3769 inline void calculateMappingVelocity_element(const int eN,
3770 const int k,
3771 double* meshVelocity_dof,
3772 int* mesh_l2g,
3773 //xt::pyarray<double>& mesh_trial_ref,
3774 double* mesh_trial_ref,
3775 double& xt,
3776 double& yt)
3777 {
3778 mapping.calculateMappingVelocity_element(eN,k,meshVelocity_dof,mesh_l2g,mesh_trial_ref,xt,yt);
3779 }
3780
3781 inline void calculateMappingVelocity_element(const int eN,
3782 const int k,
3783 double* meshVelocity_dof,
3784 int* mesh_l2g,
3785 //xt::pyarray<double>& mesh_trial_ref,
3786 double* mesh_trial_ref,
3787 double& xt,
3788 double& yt,
3789 double& zt)
3790 {
3791 mapping.calculateMappingVelocity_element(eN,k,meshVelocity_dof,mesh_l2g,mesh_trial_ref,xt,yt,zt);
3792 }
3793
3794 inline
3796 const int ebN_local,
3797 const int kb,
3798 const int ebN_local_kb,
3799 double* mesh_dof,
3800 int* mesh_l2g,
3801 double* mesh_trial_trace_ref,
3802 double* mesh_grad_trial_trace_ref,
3803 double* boundaryJac_ref,
3804 double* jac,
3805 double& jacDet,
3806 double* jacInv,
3807 double* boundaryJac,
3808 double* metricTensor,
3809 double& metricTensorDetSqrt,
3810 double* normal_ref,
3811 double* normal,
3812 double& x,
3813 double& y)
3814 {
3816 ebN_local,
3817 kb,
3818 ebN_local_kb,
3819 mesh_dof,
3820 mesh_l2g,
3821 mesh_trial_trace_ref,
3822 mesh_grad_trial_trace_ref,
3823 boundaryJac_ref,
3824 jac,
3825 jacDet,
3826 jacInv,
3827 boundaryJac,
3828 metricTensor,
3829 metricTensorDetSqrt,
3830 normal_ref,
3831 normal,
3832 x,
3833 y);
3834 }
3835
3836 inline
3838 const int ebN_local,
3839 const int kb,
3840 const int ebN_local_kb,
3841 double* mesh_dof,
3842 int* mesh_l2g,
3843 double* mesh_trial_trace_ref,
3844 double* mesh_grad_trial_trace_ref,
3845 double* boundaryJac_ref,
3846 double* jac,
3847 double& jacDet,
3848 double* jacInv,
3849 double* boundaryJac,
3850 double* metricTensor,
3851 double& metricTensorDetSqrt,
3852 double* normal_ref,
3853 double* normal,
3854 double& x,
3855 double& y,
3856 double& z)
3857 {
3859 ebN_local,
3860 kb,
3861 ebN_local_kb,
3862 mesh_dof,
3863 mesh_l2g,
3864 mesh_trial_trace_ref,
3865 mesh_grad_trial_trace_ref,
3866 boundaryJac_ref,
3867 jac,
3868 jacDet,
3869 jacInv,
3870 boundaryJac,
3871 metricTensor,
3872 metricTensorDetSqrt,
3873 normal_ref,
3874 normal,
3875 x,
3876 y,
3877 z);
3878 }
3879
3880 inline
3882 const int ebN_local,
3883 const int kb,
3884 const int ebN_local_kb,
3885 double* mesh_velocity_dof,
3886 int* mesh_l2g,
3887 double* mesh_trial_trace_ref,
3888 double& xt,
3889 double& yt,
3890 double* normal,
3891 double* boundaryJac,
3892 double* metricTensor,
3893 double& metricTensorDetSqrt)
3894 {
3896 ebN_local,
3897 kb,
3898 ebN_local_kb,
3899 mesh_velocity_dof,
3900 mesh_l2g,
3901 mesh_trial_trace_ref,
3902 xt,
3903 yt,
3904 normal,
3905 boundaryJac,
3906 metricTensor,
3907 metricTensorDetSqrt);
3908 }
3909 inline
3911 const int ebN_local,
3912 const int kb,
3913 const int ebN_local_kb,
3914 double* mesh_velocity_dof,
3915 int* mesh_l2g,
3916 double* mesh_trial_trace_ref,
3917 double& xt,
3918 double& yt,
3919 double& zt,
3920 double* normal,
3921 double* boundaryJac,
3922 double* metricTensor,
3923 double& metricTensorDetSqrt)
3924 {
3926 ebN_local,
3927 kb,
3928 ebN_local_kb,
3929 mesh_velocity_dof,
3930 mesh_l2g,
3931 mesh_trial_trace_ref,
3932 xt,
3933 yt,
3934 zt,
3935 normal,
3936 boundaryJac,
3937 metricTensor,
3938 metricTensorDetSqrt);
3939 }
3940 double Stress_u_weak(double* stress, double* grad_test_dV)
3941 {
3942 return stress[sXX]*grad_test_dV[X];
3943 }
3944 double StressJacobian_u_u_weak(double* dstress, double* grad_trial, double* grad_test_dV)
3945 {
3946 return
3947 dstress[sXX*nSymTen]*grad_trial[X];
3948 }
3949
3950 double ExteriorElementBoundaryStressFlux(const double& stressFlux,const double& disp_test_dS)
3951 {
3952 return stressFlux*disp_test_dS;
3953 }
3954 double ExteriorElementBoundaryStressFluxJacobian(const double& dstressFlux,const double& disp_test_dS)
3955 {
3956 return dstressFlux*disp_test_dS;
3957 }
3958};
3959#endif
Double q
Definition Headers.h:81
Double L
Definition Headers.h:72
Double r
Definition Headers.h:83
Double H
Definition Headers.h:65
Double f
Definition Headers.h:64
Double u
Definition Headers.h:89
Int num
Definition Headers.h:32
Double * z
Definition Headers.h:49
Double v
Definition Headers.h:95
double InteriorNumericalAdvectiveFluxJacobian(const double &dflux_left, const double &v)
double Advection_weak(const double f[1], const double grad_w_dV[1])
void DOFaverage(const double *dof, const int *l2g_element, double &val)
double SimpleDiffusionJacobian_weak(int *rowptr, int *colind, double *a, const double grad_v[1], const double grad_w_dV[1])
double SubgridError(const double &error, const double &Lstar_w_dV)
double ExteriorElementBoundaryScalarDiffusionAdjoint(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &u, const double &bc_u, const double normal[1], const double &a, const double grad_w_dS[1])
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y, double &z)
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
double Advection_adjoint(const double df[1], const double grad_w_dV[1])
void bdf(const double &alpha, const double &beta, const double &m, const double &dm, double &mt, double &dmt)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[1], double &gradNorm, double &gradNorm_last, double &numDiff)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double G[1 *1], const double &strong_residual, const double vel[1], const double grad_u[1], double &numDiff)
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double &zt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &uref, const double &beta, const double G[1 *1], const double &G_dd_G, const double &strong_residual, const double grad_u[1], double &numDiff)
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y)
double Diffusion_weak(int *rowptr, int *colind, double *a, const double grad_phi[1], const double grad_w_dV[1])
void backwardEuler(const double &dt, const double &m_old, const double &m, const double &dm, double &mt, double &dmt)
double SubgridErrorJacobian(const double &derror, const double &Lstar_w_dV)
double NumericalDiffusion(const double &numDiff, const double grad_u[1], const double grad_w_dV[1])
double ExteriorElementBoundaryStressFluxJacobian(const double &dstressFlux, const double &disp_test_dS)
double ExteriorNumericalAdvectiveFluxJacobian(const double &dflux_left, const double &v)
void calculateMappingVelocity_element(const int eN, const int k, double *meshVelocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt, double &zt)
double DiffusionJacobian_weak(int *rowptr, int *colind, double *a, double *da, const double grad_phi[1], const double grad_w_dV[1], const double &dphi, const double &v, const double grad_v[1])
void calculateG(double *jacInv, double *G, double &G_dd_G, double &tr_G)
double Advection_strong(const double df[1], const double grad_u[1])
double HamiltonianJacobian_weak(const double dH[1], const double grad_v[1], const double &w_dV)
double ExteriorElementBoundaryDiffusionAdjoint(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &u, const double &bc_u, const double normal[1], int *rowptr, int *colind, double *a, const double grad_w_dS[1])
double MassJacobian_weak(const double &dmt, const double &v, const double &w_dV)
void hessTrialFromRef(const double *hess_trial_ref, const double *jacInv, double *hess_trial)
double HamiltonianJacobian_strong(const double dH[1], const double grad_v[1])
void gradFromElementDOF(const double *dof, const double *grad_trial, double *grad)
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
CompKernelSpaceMapping< 1, NDOF_MESH_TRIAL_ELEMENT > mapping
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
double NumericalDiffusionJacobian(const double &numDiff, const double grad_v[1], const double grad_w_dV[1])
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double G[1 *1], const double &strong_residual, const double grad_u[1], double &numDiff)
double AdvectionJacobian_weak(const double df[1], const double &v, const double grad_w_dV[1])
void calculateMappingVelocity_element(const int eN, const int k, double *meshVelocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt)
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
double pressureProjection_weak(const double &viscosity, const double &p, const double &p_avg, const double &q, const double &dV)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[1], double &numDiff)
double InteriorElementBoundaryFlux(const double &flux, const double &w_dS)
double AdvectionJacobian_strong(const double df[1], const double grad_v[1])
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y, double &z)
double ExteriorElementBoundaryFlux(const double &flux, const double &w_dS)
double ExteriorElementBoundaryDiffusionAdjointJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &v, const double normal[1], int *rowptr, int *colind, double *a, const double grad_w_dS[1])
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y)
void bdfC2(const double &alpha, const double &beta, const double &m, const double &dm, const double &dm2, double &mt, double &dmt, double &dm2t)
void hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
double Hamiltonian_strong(const double dH[1], const double grad_u[1])
void gradTrialFromRef(const double *grad_trial_ref, const double *jacInv, double *grad_trial)
double ReactionJacobian_weak(const double &dr, const double &v, const double &w_dV)
double ExteriorElementBoundaryScalarDiffusionAdjointJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &v, const double normal[1], const double &a, const double grad_w_dS[1])
void gradTestFromRef(const double *grad_test_ref, const double *jacInv, double *grad_test)
double ExteriorElementBoundaryStressFlux(const double &stressFlux, const double &disp_test_dS)
void valFromElementDOF(const double *dof, const double *trial_ref, double &val)
double Hamiltonian_adjoint(const double dH[1], const double grad_w_dV[1])
double StressJacobian_u_u_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double NumericalDiffusion(const double &numDiff, const double grad_u[2], const double grad_w_dV[2])
void gradTestFromRef(const double *grad_test_ref, const double *jacInv, double *grad_test)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[2], double &gradNorm, double &gradNorm_last, double &numDiff)
CompKernelSpaceMapping< 2, NDOF_MESH_TRIAL_ELEMENT > mapping
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &uref, const double &beta, const double G[2 *2], const double &G_dd_G, const double &strong_residual, const double grad_u[2], double &numDiff)
double NumericalDiffusionJacobian(const double &numDiff, const double grad_v[2], const double grad_w_dV[2])
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
double ExteriorElementBoundaryScalarDiffusionAdjoint(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &u, const double &bc_u, const double normal[2], const double &a, const double grad_w_dS[2])
double ExteriorElementBoundaryDiffusionAdjoint(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &u, const double &bc_u, const double normal[2], int *rowptr, int *colind, double *a, const double grad_w_dS[2])
void calculateMappingVelocity_element(const int eN, const int k, double *meshVelocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt)
void calculateG(double *jacInv, double *G, double &G_dd_G, double &tr_G)
void hessTrialFromRef(const double *hess_trial_ref, const double *jacInv, double *hess_trial)
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
double InteriorElementBoundaryFlux(const double &flux, const double &w_dS)
void bdf(const double &alpha, const double &beta, const double &m, const double &dm, double &mt, double &dmt)
double pressureProjection_weak(const double &viscosity, const double &p, const double &p_avg, const double &q, const double &dV)
double SimpleDiffusionJacobian_weak(int *rowptr, int *colind, double *a, const double grad_v[2], const double grad_w_dV[2])
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double &zt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
double StressJacobian_v_v_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double Advection_weak(const double f[2], const double grad_w_dV[2])
double DiffusionJacobian_weak(int *rowptr, int *colind, double *a, double *da, const double grad_phi[2], const double grad_w_dV[2], const double &dphi, const double &v, const double grad_v[2])
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
double StressJacobian_u_u_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double Hamiltonian_strong(const double dH[2], const double grad_u[2])
void hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y)
double ExteriorElementBoundaryFlux(const double &flux, const double &w_dS)
double ExteriorElementBoundaryStressFluxJacobian(const double &dstressFlux, const double &disp_test_dS)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[2], double &numDiff)
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y, double &z)
double Advection_strong(const double df[2], const double grad_u[2])
double StressJacobian_u_v_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double HamiltonianJacobian_strong(const double dH[2], const double grad_v[2])
double ReactionJacobian_weak(const double &dr, const double &v, const double &w_dV)
double AdvectionJacobian_weak(const double df[2], const double &v, const double grad_w_dV[2])
double SubgridErrorJacobian(const double &derror, const double &Lstar_w_dV)
double Advection_adjoint(const double df[2], const double grad_w_dV[2])
double ExteriorElementBoundaryDiffusionAdjointJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &v, const double normal[2], int *rowptr, int *colind, double *a, const double grad_w_dS[2])
void gradFromElementDOF(const double *dof, const double *grad_trial, double *grad)
void DOFaverage(const double *dof, const int *l2g_element, double &val)
double ExteriorNumericalAdvectiveFluxJacobian(const double &dflux_left, const double &v)
void bdfC2(const double &alpha, const double &beta, const double &m, const double &dm, const double &dm2, double &mt, double &dmt, double &dm2t)
double InteriorNumericalAdvectiveFluxJacobian(const double &dflux_left, const double &v)
void valFromElementDOF(const double *dof, const double *trial_ref, double &val)
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y)
double Diffusion_weak(int *rowptr, int *colind, double *a, const double grad_phi[2], const double grad_w_dV[2])
double MassJacobian_weak(const double &dmt, const double &v, const double &w_dV)
double StressJacobian_v_u_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double ExteriorElementBoundaryStressFlux(const double &stressFlux, const double &disp_test_dS)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double G[2 *2], const double &strong_residual, const double grad_u[2], double &numDiff)
void calculateMappingVelocity_element(const int eN, const int k, double *meshVelocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt, double &zt)
void gradTrialFromRef(const double *grad_trial_ref, const double *jacInv, double *grad_trial)
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y, double &z)
void backwardEuler(const double &dt, const double &m_old, const double &m, const double &dm, double &mt, double &dmt)
double AdvectionJacobian_strong(const double df[2], const double grad_v[2])
double Hamiltonian_adjoint(const double dH[2], const double grad_w_dV[2])
double SubgridError(const double &error, const double &Lstar_w_dV)
double ExteriorElementBoundaryScalarDiffusionAdjointJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &v, const double normal[2], const double &a, const double grad_w_dS[2])
double HamiltonianJacobian_weak(const double dH[2], const double grad_v[2], const double &w_dV)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double G[2 *2], const double &strong_residual, const double vel[2], const double grad_u[2], double &numDiff)
double ExteriorElementBoundaryStressFlux(const double &stressFlux, const double &disp_test_dS)
double Advection_strong(const double df[NSPACE], const double grad_u[NSPACE])
double pressureProjection_weak(const double &viscosity, const double &p, const double &p_avg, const double &q, const double &dV)
void hessTrialFromRef(const double *hess_trial_ref, const double *jacInv, double *hess_trial)
double AdvectionJacobian_strong(const double df[NSPACE], const double grad_v[NSPACE])
double Reaction_adjoint(const double &dr, const double &w_dV)
const int ZY
double Stress_v_weak(double *stress, double *grad_test_dV)
double ExteriorElementBoundaryScalarDiffusionAdjoint(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &u, const double &bc_u, const double normal[NSPACE], const double &a, const double grad_w_dS[NSPACE])
CompKernelSpaceMapping< NSPACE, NDOF_MESH_TRIAL_ELEMENT > mapping
double ExteriorElementBoundaryFlux(const double &flux, const double &w_dS)
const int YZ
const int sYY
const int sYZ
double ReactionJacobian_weak(const double &dr, const double &v, const double &w_dV)
void calculateG(double *jacInv, double *G, double &G_dd_G, double &tr_G)
double StressJacobian_v_v_weak(double *dstress, double *grad_trial, double *grad_test_dV)
void valFromElementDOF(const double *dof, const double *trial_ref, double &val)
const int nSymTen
const int HXHY
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[NSPACE], double &numDiff)
double StressJacobian_u_u_weak(double *dstress, double *grad_trial, double *grad_test_dV)
void bdf(const double &alpha, const double &beta, const double &m, const double &dm, double &mt, double &dmt)
double Stress_w_weak(double *stress, double *grad_test_dV)
double Mass_adjoint(const double &dmt, const double &w_dV)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double G[NSPACE *NSPACE], const double &strong_residual, const double grad_u[NSPACE], double &numDiff)
const int YHX
void gradTestFromRef(const double *grad_test_ref, const double *jacInv, double *grad_test)
const int sYX
void calculateMappingVelocity_element(const int eN, const int k, double *meshVelocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt, double &zt)
double InteriorElementBoundaryFlux(const double &flux, const double &w_dS)
double Advection_weak(const double f[NSPACE], const double grad_w_dV[NSPACE])
double AdvectionJacobian_weak(const double df[NSPACE], const double &v, const double grad_w_dV[NSPACE])
void hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
const int ZHX
double HamiltonianJacobian_weak(const double dH[NSPACE], const double grad_v[NSPACE], const double &w_dV)
double Advection_adjoint(const double df[NSPACE], const double grad_w_dV[NSPACE])
const int HXHX
double ExteriorNumericalAdvectiveFluxJacobian(const double &dflux_left, const double &v)
double StressJacobian_w_u_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double MassJacobian_weak(const double &dmt, const double &v, const double &w_dV)
const int ZHY
const int sZY
double Hamiltonian_strong(const double dH[NSPACE], const double grad_u[NSPACE])
double DiffusionJacobian_weak(int *rowptr, int *colind, double *a, double *da, const double grad_phi[NSPACE], const double grad_w_dV[NSPACE], const double &dphi, const double &v, const double grad_v[NSPACE])
double Mass_weak(const double &mt, const double &w_dV)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &uref, const double &beta, const double G[NSPACE *NSPACE], const double &G_dd_G, const double &strong_residual, const double grad_u[NSPACE], double &numDiff)
const int ZZ
double ReactionJacobian_strong(const double &dr, const double &v)
const int YY
void gradTrialFromRef(const double *grad_trial_ref, const double *jacInv, double *grad_trial)
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
const int HYHX
double ExteriorElementBoundaryDiffusionAdjointJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &v, const double normal[NSPACE], int *rowptr, int *colind, double *a, const double grad_w_dS[NSPACE])
double ExteriorElementBoundaryScalarDiffusionAdjointJacobian(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &v, const double normal[NSPACE], const double &a, const double grad_w_dS[NSPACE])
double NumericalDiffusion(const double &numDiff, const double grad_u[NSPACE], const double grad_w_dV[NSPACE])
const int HYHY
const int sZX
double HamiltonianJacobian_strong(const double dH[NSPACE], const double grad_v[NSPACE])
double NumericalDiffusionJacobian(const double &numDiff, const double grad_v[NSPACE], const double grad_w_dV[NSPACE])
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double &zt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
double StressJacobian_u_v_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double Hamiltonian_adjoint(const double dH[NSPACE], const double grad_w_dV[NSPACE])
double Reaction_weak(const double &r, const double &w_dV)
double Stress_u_weak(double *stress, double *grad_test_dV)
double SubgridError(const double &error, const double &Lstar_w_dV)
double StressJacobian_v_u_weak(double *dstress, double *grad_trial, double *grad_test_dV)
void backwardEuler(const double &dt, const double &m_old, const double &m, const double &dm, double &mt, double &dmt)
double ExteriorElementBoundaryDiffusionAdjoint(const int &isDOFBoundary, const int &isFluxBoundary, const double &sigma, const double &u, const double &bc_u, const double normal[NSPACE], int *rowptr, int *colind, double *a, const double grad_w_dS[NSPACE])
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
double StressJacobian_v_w_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double Diffusion_weak(int *rowptr, int *colind, double *a, const double grad_phi[NSPACE], const double grad_w_dV[NSPACE])
const int sZZ
void gradFromElementDOF(const double *dof, const double *grad_trial, double *grad)
const int XHX
const int YX
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y, double &z)
void DOFaverage(const double *dof, const int *l2g_element, double &val)
double StressJacobian_w_w_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double Mass_strong(const double &mt)
const int sXX
double Hamiltonian_weak(const double &H, const double &w_dV)
const int Y
double Reaction_strong(const double &r)
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double &elementDiameter, const double &strong_residual, const double grad_u[NSPACE], double &gradNorm, double &gradNorm_last, double &numDiff)
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y, double &z)
double StressJacobian_w_v_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double SimpleDiffusionJacobian_weak(int *rowptr, int *colind, double *a, const double grad_v[NSPACE], const double grad_w_dV[NSPACE])
const int XY
const int sXY
const int XHY
void bdfC2(const double &alpha, const double &beta, const double &m, const double &dm, const double &dm2, double &mt, double &dmt, double &dm2t)
const int XX
void calculateGScale(double *G, double *v, double &h)
double InteriorNumericalAdvectiveFluxJacobian(const double &dflux_left, const double &v)
double SubgridErrorJacobian(const double &derror, const double &Lstar_w_dV)
double MassJacobian_strong(const double &dmt, const double &v)
const int sXZ
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
double StressJacobian_u_w_weak(double *dstress, double *grad_trial, double *grad_test_dV)
double ExteriorElementBoundaryStressFluxJacobian(const double &dstressFlux, const double &disp_test_dS)
const int YHY
const int X
const int ZX
const int Z
void calculateNumericalDiffusion(const double &shockCapturingDiffusion, const double G[NSPACE *NSPACE], const double &strong_residual, const double vel[NSPACE], const double grad_u[NSPACE], double &numDiff)
const int XZ
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y, double &z)
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y, double &z)
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x)
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x)
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt, double &zt)
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double &zt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
void hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt)
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y, double &z)
Definition CompKernel.h:964
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y)
Definition CompKernel.h:666
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
Definition CompKernel.h:721
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
Definition CompKernel.h:853
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y)
Definition CompKernel.h:763
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
Definition CompKernel.h:900
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt)
Definition CompKernel.h:738
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double &zt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y, double &z)
Definition CompKernel.h:920
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
Definition CompKernel.h:893
void hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
Definition CompKernel.h:909
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt, double &zt)
Definition CompKernel.h:946
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
Definition CompKernel.h:608
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt)
Definition CompKernel.h:549
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt, double &zt)
Definition CompKernel.h:312
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y, double &z)
Definition CompKernel.h:338
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y, double &z)
Definition CompKernel.h:230
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
Definition CompKernel.h:495
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
Definition CompKernel.h:502
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
Definition CompKernel.h:295
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y)
Definition CompKernel.h:522
void hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
Definition CompKernel.h:511
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double &zt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
Definition CompKernel.h:445
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y)
Definition CompKernel.h:567
void calculateMapping_element(const int eN, const int k, double *mesh_dof, int *mesh_l2g, double *mesh_trial_ref, double *mesh_grad_trial_ref, double *jac, double &jacDet, double *jacInv, double &x, double &y, double &z)
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
Definition CompKernel.h:165
void calculateMappingVelocity_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double &xt, double &yt, double &zt, double *normal, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt)
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt, double &zt)
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
void hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
Definition CompKernel.h:181
void calculateMapping_elementBoundary(const int eN, const int ebN_local, const int kb, const int ebN_local_kb, double *mesh_dof, int *mesh_l2g, double *mesh_trial_trace_ref, double *mesh_grad_trial_trace_ref, double *boundaryJac_ref, double *jac, double &jacDet, double *jacInv, double *boundaryJac, double *metricTensor, double &metricTensorDetSqrt, double *normal_ref, double *normal, double &x, double &y, double &z)
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
Definition CompKernel.h:172
const int X
Definition CompKernel.h:81
const int nSymTen
Definition CompKernel.h:84
const int HXHX
Definition CompKernel.h:86
const int XHX
Definition CompKernel.h:85
const int sXX
Definition CompKernel.h:83
const int XX
Definition CompKernel.h:82
const int YHX
Definition CompKernel.h:61
const int sYY
Definition CompKernel.h:58
const int sXY
Definition CompKernel.h:57
const int HXHX
Definition CompKernel.h:62
const int sXX
Definition CompKernel.h:57
const int YX
Definition CompKernel.h:56
const int YY
Definition CompKernel.h:56
const int XHX
Definition CompKernel.h:60
const int XHY
Definition CompKernel.h:60
const int sYX
Definition CompKernel.h:58
const int XX
Definition CompKernel.h:55
const int Y
Definition CompKernel.h:54
const int nSymTen
Definition CompKernel.h:59
const int X
Definition CompKernel.h:54
const int XY
Definition CompKernel.h:55
const int YHY
Definition CompKernel.h:61
const int XY
Definition CompKernel.h:19
const int XHX
Definition CompKernel.h:26
const int ZHY
Definition CompKernel.h:28
const int ZX
Definition CompKernel.h:21
const int YHX
Definition CompKernel.h:27
const int sXZ
Definition CompKernel.h:22
const int Z
Definition CompKernel.h:18
const int sZX
Definition CompKernel.h:24
const int ZZ
Definition CompKernel.h:21
const int YX
Definition CompKernel.h:20
const int sXX
Definition CompKernel.h:22
const int HYHY
Definition CompKernel.h:30
const int HXHX
Definition CompKernel.h:29
const int YY
Definition CompKernel.h:20
const int sYX
Definition CompKernel.h:23
const int XZ
Definition CompKernel.h:19
const int YZ
Definition CompKernel.h:20
const int sXY
Definition CompKernel.h:22
const int ZY
Definition CompKernel.h:21
const int Y
Definition CompKernel.h:18
const int X
Definition CompKernel.h:18
const int XHY
Definition CompKernel.h:26
const int sYY
Definition CompKernel.h:23
const int nSymTen
Definition CompKernel.h:25
const int sZY
Definition CompKernel.h:24
const int HYHX
Definition CompKernel.h:30
const int ZHX
Definition CompKernel.h:28
const int XX
Definition CompKernel.h:19
const int sYZ
Definition CompKernel.h:23
const int HXHY
Definition CompKernel.h:29
const int sZZ
Definition CompKernel.h:24
const int YHY
Definition CompKernel.h:27
double df(double C, double b, double a, int q, int r)
void vel(double rS, double norm_v, double r, double theta, double *vR, double *vTHETA)