98template<
int NSPACE,
int NDOF_MESH_TRIAL_ELEMENT>
107 double* mesh_trial_ref,
108 double* mesh_grad_trial_ref,
120 double* mesh_trial_ref,
124 double* mesh_velocity_dof,
127 double* mesh_trial_ref,
134 const int ebN_local_kb,
137 double* mesh_trial_trace_ref,
138 double* mesh_grad_trial_trace_ref,
139 double* boundaryJac_ref,
144 double* metricTensor,
145 double& metricTensorDetSqrt,
154 const int ebN_local_kb,
155 double* mesh_velocity_dof,
157 double* mesh_trial_trace_ref,
163 double* metricTensor,
164 double& metricTensorDetSqrt);
165 inline void valFromDOF(
const double* dof,
const int* l2g_element,
const double* trial_ref,
double& val)
168 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
169 val+=dof[l2g_element[j]]*trial_ref[j];
172 inline void gradFromDOF(
const double* dof,
const int* l2g_element,
const double* grad_trial,
double* grad)
174 for(
int I=0;I<NSPACE;I++)
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];
181 inline void hessFromDOF(
const double* dof,
const int* l2g_element,
const double* hess_trial,
double* hess)
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];
195template<
int NDOF_MESH_TRIAL_ELEMENT>
235 double* mesh_trial_ref,
236 double* mesh_grad_trial_ref,
244 double Grad_x[3],Grad_y[3],Grad_z[3],oneOverJacDet;
250 for (
int I=0;I<3;I++)
252 Grad_x[I]=0.0;Grad_y[I]=0.0;Grad_z[I]=0.0;
254 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
256 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+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++)
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];
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]);
300 double* mesh_trial_ref,
304 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
306 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
308 h += h_dof[mesh_l2g[eN_j]]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
314 double* mesh_velocity_dof,
317 double* mesh_trial_ref,
325 xt=0.0;yt=0.0;zt=0.0;
326 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
328 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+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];
341 const int ebN_local_kb,
344 double* mesh_trial_trace_ref,
345 double* mesh_grad_trial_trace_ref,
346 double* boundaryJac_ref,
351 double* metricTensor,
352 double& metricTensorDetSqrt,
359 const int ebN_local_kb_nSpace = ebN_local_kb*3,
360 ebN_local_kb_nSpace_nSpacem1 = ebN_local_kb*3*2;
362 double Grad_x_ext[3],Grad_y_ext[3],Grad_z_ext[3],oneOverJacDet,norm_normal=0.0;
367 for (
int I=0;I<3;I++)
373 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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++)
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];
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];
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]);
414 for (
int I=0;I<3;I++)
416 for (
int I=0;I<3;I++)
418 for (
int J=0;J<3;J++)
420 normal[I] += jacInv[J*3+I]*normal_ref[ebN_local_kb_nSpace+J];
422 norm_normal+=normal[I]*normal[I];
424 norm_normal = sqrt(norm_normal);
425 for (
int I=0;I<3;I++)
427 normal[I] /= norm_normal;
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];
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];
442 metricTensorDetSqrt=sqrt(metricTensor[
HXHX]*metricTensor[
HYHY]- metricTensor[
HXHY]*metricTensor[
HYHX]);
448 const int ebN_local_kb,
449 double* mesh_velocity_dof,
451 double* mesh_trial_trace_ref,
457 double* metricTensor,
458 double& metricTensorDetSqrt)
465 xt=0.0;yt=0.0;zt=0.0;
466 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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];
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));
495 inline void valFromDOF(
const double* dof,
const int* l2g_element,
const double* trial_ref,
double& val)
498 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
499 val+=dof[l2g_element[j]]*trial_ref[j];
502 inline void gradFromDOF(
const double* dof,
const int* l2g_element,
const double* grad_trial,
double* grad)
506 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
508 grad[I] += dof[l2g_element[j]]*grad_trial[j*3+I];
511 inline void hessFromDOF(
const double* dof,
const int* l2g_element,
const double* hess_trial,
double* hess)
516 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
519 hess[I*3+J] += dof[l2g_element[j]]*hess_trial[j*9+I*3+J];
527 double* mesh_trial_ref,
528 double* mesh_grad_trial_ref,
551 double* mesh_velocity_dof,
554 double* mesh_trial_ref,
570 const int ebN_local_kb,
573 double* mesh_trial_trace_ref,
574 double* mesh_grad_trial_trace_ref,
575 double* boundaryJac_ref,
580 double* metricTensor,
581 double& metricTensorDetSqrt,
593 mesh_trial_trace_ref,
594 mesh_grad_trial_trace_ref,
611 const int ebN_local_kb,
612 double* mesh_velocity_dof,
614 double* mesh_trial_trace_ref,
619 double* metricTensor,
620 double& metricTensorDetSqrt)
628 mesh_trial_trace_ref,
635 metricTensorDetSqrt);
641template<
int NDOF_MESH_TRIAL_ELEMENT>
671 double* mesh_trial_ref,
672 double* mesh_grad_trial_ref,
679 double Grad_x[2],Grad_y[2],oneOverJacDet;
685 for (
int I=0;I<2;I++)
687 Grad_x[I]=0.0;Grad_y[I]=0.0;
689 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
691 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+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++)
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];
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];
726 double* mesh_trial_ref,
730 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
732 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+j;
734 h += h_dof[mesh_l2g[eN_j]]*mesh_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j];
740 double* mesh_velocity_dof,
743 double* mesh_trial_ref,
751 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
753 int eN_j=eN*NDOF_MESH_TRIAL_ELEMENT+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];
766 const int ebN_local_kb,
769 double* mesh_trial_trace_ref,
770 double* mesh_grad_trial_trace_ref,
771 double* boundaryJac_ref,
776 double* metricTensor,
777 double& metricTensorDetSqrt,
783 const int ebN_local_kb_nSpace = ebN_local_kb*2,
784 ebN_local_kb_nSpace_nSpacem1 = ebN_local_kb*2*1;
786 double Grad_x_ext[2],Grad_y_ext[2],oneOverJacDet,norm_normal=0.0;
791 for (
int I=0;I<2;I++)
796 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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;
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++)
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];
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];
829 for (
int I=0;I<2;I++)
831 for (
int I=0;I<2;I++)
833 for (
int J=0;J<2;J++)
835 normal[I] += jacInv[J*2+I]*normal_ref[ebN_local_kb_nSpace+J];
837 norm_normal+=normal[I]*normal[I];
839 norm_normal = sqrt(norm_normal);
840 for (
int I=0;I<2;I++)
842 normal[I] /= norm_normal;
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];
848 metricTensor[
HXHX] = boundaryJac[
XHX]*boundaryJac[
XHX]+boundaryJac[
YHX]*boundaryJac[
YHX];
850 metricTensorDetSqrt=sqrt(metricTensor[
HXHX]);
856 const int ebN_local_kb,
857 double* mesh_velocity_dof,
859 double* mesh_trial_trace_ref,
864 double* metricTensor,
865 double& metricTensorDetSqrt)
873 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
875 int eN_j = eN*NDOF_MESH_TRIAL_ELEMENT+j;
876 int ebN_local_kb_j = ebN_local_kb*NDOF_MESH_TRIAL_ELEMENT+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];
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));
893 inline void valFromDOF(
const double* dof,
const int* l2g_element,
const double* trial_ref,
double& val)
896 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
897 val+=dof[l2g_element[j]]*trial_ref[j];
900 inline void gradFromDOF(
const double* dof,
const int* l2g_element,
const double* grad_trial,
double* grad)
904 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
906 grad[I] += dof[l2g_element[j]]*grad_trial[j*2+I];
909 inline void hessFromDOF(
const double* dof,
const int* l2g_element,
const double* hess_trial,
double* hess)
914 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
917 hess[I*2+J] += dof[l2g_element[j]]*hess_trial[j*4+I*2+J];
925 double* mesh_trial_ref,
926 double* mesh_grad_trial_ref,
948 double* mesh_velocity_dof,
951 double* mesh_trial_ref,
967 const int ebN_local_kb,
970 double* mesh_trial_trace_ref,
971 double* mesh_grad_trial_trace_ref,
972 double* boundaryJac_ref,
977 double* metricTensor,
978 double& metricTensorDetSqrt,
991 mesh_trial_trace_ref,
992 mesh_grad_trial_trace_ref,
1010 const int ebN_local,
1012 const int ebN_local_kb,
1013 double* mesh_velocity_dof,
1015 double* mesh_trial_trace_ref,
1020 double* boundaryJac,
1021 double* metricTensor,
1022 double& metricTensorDetSqrt)
1030 mesh_trial_trace_ref,
1036 metricTensorDetSqrt);
1041template<
int NDOF_MESH_TRIAL_ELEMENT>
1064 double* mesh_trial_ref,
1065 double* mesh_grad_trial_ref,
1071 double Grad_x[1],Grad_y[1],oneOverJacDet;
1077 for (
int I=0;I<1;I++)
1081 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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++)
1087 Grad_x[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_ref[k*NDOF_MESH_TRIAL_ELEMENT+j+I];
1090 jac[
XX] = Grad_x[
X];
1092 oneOverJacDet = 1.0/jacDet;
1093 jacInv[
XX] = oneOverJacDet;
1101 double* mesh_trial_ref,
1105 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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];
1114 double* mesh_velocity_dof,
1117 double* mesh_trial_ref,
1124 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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];
1132 const int ebN_local,
1134 const int ebN_local_kb,
1137 double* mesh_trial_trace_ref,
1138 double* mesh_grad_trial_trace_ref,
1139 double* boundaryJac_ref,
1143 double* boundaryJac,
1144 double* metricTensor,
1145 double& metricTensorDetSqrt,
1150 const int ebN_local_kb_nSpace = ebN_local_kb*1;
1152 double Grad_x_ext[1],oneOverJacDet,norm_normal=0.0;
1157 for (
int I=0;I<1;I++)
1159 Grad_x_ext[I] = 0.0;
1161 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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++)
1169 Grad_x_ext[I] += mesh_dof[mesh_l2g[eN_j]*3+0]*mesh_grad_trial_trace_ref[ebN_local_kb_j_nSpace+I];
1173 jac[
XX] = Grad_x_ext[
X];
1175 oneOverJacDet = 1.0/jacDet;
1176 jacInv[
XX] = oneOverJacDet;
1179 for (
int I=0;I<1;I++)
1181 for (
int I=0;I<1;I++)
1183 for (
int J=0;J<1;J++)
1185 normal[I] += jacInv[J+I]*normal_ref[ebN_local_kb_nSpace+J];
1187 norm_normal+=normal[I]*normal[I];
1189 norm_normal = sqrt(norm_normal);
1190 for (
int I=0;I<1;I++)
1192 normal[I] /= norm_normal;
1195 boundaryJac[
XHX] = 1.0;
1197 metricTensor[
HXHX] = 1.0;
1199 metricTensorDetSqrt=1.0;
1203 const int ebN_local,
1205 const int ebN_local_kb,
1206 double* mesh_velocity_dof,
1208 double* mesh_trial_trace_ref,
1211 double* boundaryJac,
1212 double* metricTensor,
1213 double& metricTensorDetSqrt)
1219 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
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];
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));
1236 inline void valFromDOF(
const double* dof,
const int* l2g_element,
const double* trial_ref,
double& val)
1239 for (
int j=0;j<NDOF_MESH_TRIAL_ELEMENT;j++)
1240 val+=dof[l2g_element[j]]*trial_ref[j];
1243 inline void gradFromDOF(
const double* dof,
const int* l2g_element,
const double* grad_trial,
double* grad)
1245 for(
int I=0;I<1;I++)
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];
1252 inline void hessFromDOF(
const double* dof,
const int* l2g_element,
const double* hess_trial,
double* hess)
1254 for(
int I=0;I<1;I++)
1255 for(
int J=0;J<1;J++)
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];
1268 double* mesh_trial_ref,
1269 double* mesh_grad_trial_ref,
1282 mesh_grad_trial_ref,
1290 double* mesh_velocity_dof,
1293 double* mesh_trial_ref,
1306 const int ebN_local,
1308 const int ebN_local_kb,
1311 double* mesh_trial_trace_ref,
1312 double* mesh_grad_trial_trace_ref,
1313 double* boundaryJac_ref,
1317 double* boundaryJac,
1318 double* metricTensor,
1319 double& metricTensorDetSqrt,
1332 mesh_trial_trace_ref,
1333 mesh_grad_trial_trace_ref,
1340 metricTensorDetSqrt,
1350 const int ebN_local,
1352 const int ebN_local_kb,
1353 double* mesh_velocity_dof,
1355 double* mesh_trial_trace_ref,
1360 double* boundaryJac,
1361 double* metricTensor,
1362 double& metricTensorDetSqrt)
1370 mesh_trial_trace_ref,
1375 metricTensorDetSqrt);
1380template<
int NSPACE,
int NDOF_MESH_TRIAL_ELEMENT,
int NDOF_TRIAL_ELEMENT,
int NDOF_TEST_ELEMENT>
1417 inline void calculateG(
double* jacInv,
double* G,
double& G_dd_G,
double& tr_G)
1421 for (
int I=0;I<NSPACE;I++)
1422 for (
int J=0;J<NSPACE;J++)
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];
1430 for (
int I=0;I<NSPACE;I++)
1432 tr_G += G[I*NSPACE+I];
1433 for (
int J=0;J<NSPACE;J++)
1435 G_dd_G += G[I*NSPACE+J]*G[I*NSPACE+J];
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);
1447 inline void valFromDOF(
const double* dof,
const int* l2g_element,
const double* trial_ref,
double& val)
1450 for (
int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1451 val+=dof[l2g_element[j]]*trial_ref[j];
1454 inline void gradFromDOF(
const double* dof,
const int* l2g_element,
const double* grad_trial,
double* grad)
1456 for(
int I=0;I<NSPACE;I++)
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];
1463 inline void hessFromDOF(
const double* dof,
const int* l2g_element,
const double* hess_trial,
double* hess)
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];
1478 for (
int j=0;j<NDOF_TRIAL_ELEMENT;j++)
1479 val+=dof[j]*trial_ref[j];
1484 for(
int I=0;I<NSPACE;I++)
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];
1491 inline void gradTrialFromRef(
const double* grad_trial_ref,
const double* jacInv,
double* grad_trial)
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];
1513 inline void DOFaverage (
const double* dof,
const int* l2g_element,
double& val)
1517 for (
int j=0; j<NDOF_MESH_TRIAL_ELEMENT; j++)
1518 val+=dof[l2g_element[j]];
1520 val /= NDOF_MESH_TRIAL_ELEMENT;
1523 inline void hessTrialFromRef(
const double* hess_trial_ref,
const double* jacInv,
double* hess_trial)
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];
1538 inline void gradTestFromRef(
const double* grad_test_ref,
const double* jacInv,
double* grad_test)
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];
1549 inline void backwardEuler(
const double& dt,
const double& m_old,
const double& m,
const double& dm,
double& mt,
double& dmt)
1555 inline void bdf(
const double& alpha,
const double& beta,
const double& m,
const double& dm,
double& mt,
double& dmt)
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)
1567 inline double Mass_weak(
const double& mt,
const double& w_dV)
1616 const double& p_avg,
1620 if (viscosity==0.){
return 0.;}
1621 return (1./viscosity)*(p-p_avg)*(
q-1./3.)*dV;
1625 const double grad_w_dV[NSPACE])
1628 for(
int I=0;I<NSPACE;I++)
1629 tmp -=
f[I]*grad_w_dV[I];
1635 const double grad_w_dV[NSPACE])
1638 for(
int I=0;I<NSPACE;I++)
1639 tmp -=
df[I]*
v*grad_w_dV[I];
1644 const double grad_u[NSPACE])
1647 for(
int I=0;I<NSPACE;I++)
1648 tmp +=
df[I]*grad_u[I];
1653 const double grad_v[NSPACE])
1656 for(
int I=0;I<NSPACE;I++)
1657 tmp +=
df[I]*grad_v[I];
1662 const double grad_w_dV[NSPACE])
1665 for(
int I=0;I<NSPACE;I++)
1666 tmp -=
df[I]*grad_w_dV[I];
1677 const double grad_v[NSPACE],
1681 for(
int I=0;I<NSPACE;I++)
1682 tmp += dH[I]*grad_v[I]*w_dV;
1687 const double grad_u[NSPACE])
1690 for(
int I=0;I<NSPACE;I++)
1691 tmp += dH[I]*grad_u[I];
1696 const double grad_v[NSPACE])
1699 for(
int I=0;I<NSPACE;I++)
1700 tmp += dH[I]*grad_v[I];
1705 const double grad_w_dV[NSPACE])
1708 for(
int I=0;I<NSPACE;I++)
1709 tmp -= dH[I]*grad_w_dV[I];
1716 const double grad_phi[NSPACE],
1717 const double grad_w_dV[NSPACE])
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];
1730 const double grad_phi[NSPACE],
1731 const double grad_w_dV[NSPACE],
1734 const double grad_v[NSPACE])
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++)
1740 daProduct += da[m]*grad_phi[colind[m]]*grad_w_dV[I];
1741 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
1743 return daProduct*
v+dphiProduct*dphi;
1749 const double grad_v[NSPACE],
1750 const double grad_w_dV[NSPACE])
1752 double dphiProduct=0.0;
1753 for (
int I=0;I<NSPACE;I++)
1754 for (
int m=rowptr[I];m<rowptr[I+1];m++)
1756 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
1792 const double& elementDiameter,
1793 const double& strong_residual,
1794 const double grad_u[NSPACE],
1801 h = elementDiameter;
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);
1812 const double G[NSPACE*NSPACE],
1813 const double& strong_residual,
1814 const double grad_u[NSPACE],
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));
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],
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];
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);
1845 const double G[NSPACE*NSPACE],
1846 const double& strong_residual,
1847 const double vel[NSPACE],
1848 const double grad_u[NSPACE],
1851 double den1 = 0.0,den2=0.0, nom=0.0;
1852 for (
int I=0;I<NSPACE;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];
1859 numDiff = shockCapturingDiffusion*fabs(strong_residual)*(sqrt(nom/(den1*den2 + 1.0e-12)));
1864 const double& elementDiameter,
1865 const double& strong_residual,
1866 const double grad_u[NSPACE],
1868 double& gradNorm_last,
1874 h = elementDiameter;
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);
1881 numDiff =
num/gradNorm_last;
1885 const double& Lstar_w_dV)
1887 return error*Lstar_w_dV;
1891 const double& Lstar_w_dV)
1893 return derror*Lstar_w_dV;
1897 const double grad_u[NSPACE],
1898 const double grad_w_dV[NSPACE])
1901 for (
int I=0;I<NSPACE;I++)
1902 tmp += numDiff*grad_u[I]*grad_w_dV[I];
1907 const double grad_v[NSPACE],
1908 const double grad_w_dV[NSPACE])
1911 for (
int I=0;I<NSPACE;I++)
1912 tmp += numDiff*grad_v[I]*grad_w_dV[I];
1933 return dflux_left*
v;
1939 return dflux_left*
v;
1943 const int& isFluxBoundary,
1944 const double& sigma,
1947 const double normal[NSPACE],
1949 const double grad_w_dS[NSPACE])
1952 for(
int I=0;I<NSPACE;I++)
1954 tmp += normal[I]*grad_w_dS[I];
1956 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*(
u-bc_u)*a;
1961 const int& isFluxBoundary,
1962 const double& sigma,
1964 const double normal[NSPACE],
1966 const double grad_w_dS[NSPACE])
1969 for(
int I=0;I<NSPACE;I++)
1971 tmp += normal[I]*grad_w_dS[I];
1973 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*
v*a;
1978 const int& isFluxBoundary,
1979 const double& sigma,
1982 const double normal[NSPACE],
1986 const double grad_w_dS[NSPACE])
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];
1996 const int& isFluxBoundary,
1997 const double& sigma,
1999 const double normal[NSPACE],
2003 const double grad_w_dS[NSPACE])
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];
2017 double* mesh_trial_ref,
2018 double* mesh_grad_trial_ref,
2026 mapping.
calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y,
z);
2034 double* mesh_trial_ref,
2047 double* meshVelocity_dof,
2050 double* mesh_trial_ref,
2060 const int ebN_local,
2062 const int ebN_local_kb,
2065 double* mesh_trial_trace_ref,
2066 double* mesh_grad_trial_trace_ref,
2067 double* boundaryJac_ref,
2071 double* boundaryJac,
2072 double* metricTensor,
2073 double& metricTensorDetSqrt,
2086 mesh_trial_trace_ref,
2087 mesh_grad_trial_trace_ref,
2094 metricTensorDetSqrt,
2104 const int ebN_local,
2106 const int ebN_local_kb,
2107 double* mesh_velocity_dof,
2109 double* mesh_trial_trace_ref,
2114 double* boundaryJac,
2115 double* metricTensor,
2116 double& metricTensorDetSqrt)
2124 mesh_trial_trace_ref,
2131 metricTensorDetSqrt);
2135 return stress[
sXX]*grad_test_dV[
X] + stress[
sXY]*grad_test_dV[
Y] + stress[
sXZ]*grad_test_dV[
Z];
2160 return stress[
sYX]*grad_test_dV[
X] + stress[
sYY]*grad_test_dV[
Y] + stress[
sYZ]*grad_test_dV[
Z];
2185 return stress[
sZX]*grad_test_dV[
X] + stress[
sZY]*grad_test_dV[
Y] + stress[
sZZ]*grad_test_dV[
Z];
2210 return stressFlux*disp_test_dS;
2214 return dstressFlux*disp_test_dS;
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>
2246 inline void calculateG(
double* jacInv,
double* G,
double& G_dd_G,
double& tr_G)
2250 for (
int I=0;I<2;I++)
2251 for (
int J=0;J<2;J++)
2254 for (
int K=0;K<2;K++)
2255 G[I*2+J] += jacInv[K*2+I]*jacInv[K*2+J];
2259 for (
int I=0;I<2;I++)
2262 for (
int J=0;J<2;J++)
2264 G_dd_G += G[I*2+J]*G[I*2+J];
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);
2276 inline void valFromDOF(
const double* dof,
const int* l2g_element,
const double* trial_ref,
double& val)
2279 for (
int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2280 val+=dof[l2g_element[j]]*trial_ref[j];
2283 inline void gradFromDOF(
const double* dof,
const int* l2g_element,
const double* grad_trial,
double* grad)
2285 for(
int I=0;I<2;I++)
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];
2292 inline void hessFromDOF(
const double* dof,
const int* l2g_element,
const double* hess_trial,
double* hess)
2294 for(
int I=0;I<2;I++)
2295 for(
int J=0;J<2;J++)
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];
2306 for (
int j=0;j<NDOF_TRIAL_ELEMENT;j++)
2307 val+=dof[j]*trial_ref[j];
2312 for(
int I=0;I<2;I++)
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];
2319 inline void gradTrialFromRef(
const double* grad_trial_ref,
const double* jacInv,
double* grad_trial)
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];
2341 inline void DOFaverage (
const double* dof,
const int* l2g_element,
double& val)
2345 for (
int j=0; j<NDOF_MESH_TRIAL_ELEMENT; j++)
2346 val+=dof[l2g_element[j]];
2348 val /= NDOF_MESH_TRIAL_ELEMENT;
2352 inline void hessTrialFromRef(
const double* hess_trial_ref,
const double* jacInv,
double* hess_trial)
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];
2366 inline void gradTestFromRef(
const double* grad_test_ref,
const double* jacInv,
double* grad_test)
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];
2377 inline void backwardEuler(
const double& dt,
const double& m_old,
const double& m,
const double& dm,
double& mt,
double& dmt)
2383 inline void bdf(
const double& alpha,
const double& beta,
const double& m,
const double& dm,
double& mt,
double& dmt)
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)
2396 inline double Mass_weak(
const double& mt,
const double& w_dV)
2445 const double& p_avg,
2449 if (viscosity==0.){
return 0.;}
2450 return (1./viscosity)*(p-p_avg)*(
q-1./3.)*dV;
2454 const double grad_w_dV[2])
2457 for(
int I=0;I<2;I++)
2458 tmp -=
f[I]*grad_w_dV[I];
2464 const double grad_w_dV[2])
2467 for(
int I=0;I<2;I++)
2468 tmp -=
df[I]*
v*grad_w_dV[I];
2473 const double grad_u[2])
2476 for(
int I=0;I<2;I++)
2477 tmp +=
df[I]*grad_u[I];
2482 const double grad_v[2])
2485 for(
int I=0;I<2;I++)
2486 tmp +=
df[I]*grad_v[I];
2491 const double grad_w_dV[2])
2494 for(
int I=0;I<2;I++)
2495 tmp -=
df[I]*grad_w_dV[I];
2506 const double grad_v[2],
2510 for(
int I=0;I<2;I++)
2511 tmp += dH[I]*grad_v[I]*w_dV;
2516 const double grad_u[2])
2519 for(
int I=0;I<2;I++)
2520 tmp += dH[I]*grad_u[I];
2525 const double grad_v[2])
2528 for(
int I=0;I<2;I++)
2529 tmp += dH[I]*grad_v[I];
2534 const double grad_w_dV[2])
2537 for(
int I=0;I<2;I++)
2538 tmp -= dH[I]*grad_w_dV[I];
2545 const double grad_phi[2],
2546 const double grad_w_dV[2])
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];
2559 const double grad_phi[2],
2560 const double grad_w_dV[2],
2563 const double grad_v[2])
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++)
2569 daProduct += da[m]*grad_phi[colind[m]]*grad_w_dV[I];
2570 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
2572 return daProduct*
v+dphiProduct*dphi;
2578 const double grad_v[2],
2579 const double grad_w_dV[2])
2581 double dphiProduct=0.0;
2582 for (
int I=0;I<2;I++)
2583 for (
int m=rowptr[I];m<rowptr[I+1];m++)
2585 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
2621 const double& elementDiameter,
2622 const double& strong_residual,
2623 const double grad_u[2],
2630 h = elementDiameter;
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);
2642 const double G[2*2],
2643 const double& strong_residual,
2644 const double grad_u[2],
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];
2652 numDiff = shockCapturingDiffusion*fabs(strong_residual)/(sqrt(den+1.0e-12));
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],
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];
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);
2676 const double G[2*2],
2677 const double& strong_residual,
2678 const double vel[2],
2679 const double grad_u[2],
2682 double den1 = 0.0,den2=0.0, nom=0.0;
2683 for (
int I=0;I<2;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];
2690 numDiff = shockCapturingDiffusion*fabs(strong_residual)*(sqrt(nom/(den1*den2 + 1.0e-12)));
2695 const double& elementDiameter,
2696 const double& strong_residual,
2697 const double grad_u[2],
2699 double& gradNorm_last,
2705 h = elementDiameter;
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);
2712 numDiff =
num/gradNorm_last;
2716 const double& Lstar_w_dV)
2718 return error*Lstar_w_dV;
2722 const double& Lstar_w_dV)
2724 return derror*Lstar_w_dV;
2728 const double grad_u[2],
2729 const double grad_w_dV[2])
2732 for (
int I=0;I<2;I++)
2733 tmp += numDiff*grad_u[I]*grad_w_dV[I];
2738 const double grad_v[2],
2739 const double grad_w_dV[2])
2742 for (
int I=0;I<2;I++)
2743 tmp += numDiff*grad_v[I]*grad_w_dV[I];
2764 return dflux_left*
v;
2770 return dflux_left*
v;
2774 const int& isFluxBoundary,
2775 const double& sigma,
2778 const double normal[2],
2780 const double grad_w_dS[2])
2783 for(
int I=0;I<2;I++)
2785 tmp += normal[I]*grad_w_dS[I];
2787 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*(
u-bc_u)*a;
2792 const int& isFluxBoundary,
2793 const double& sigma,
2795 const double normal[2],
2797 const double grad_w_dS[2])
2800 for(
int I=0;I<2;I++)
2802 tmp += normal[I]*grad_w_dS[I];
2804 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*
v*a;
2809 const int& isFluxBoundary,
2810 const double& sigma,
2813 const double normal[2],
2817 const double grad_w_dS[2])
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];
2827 const int& isFluxBoundary,
2828 const double& sigma,
2830 const double normal[2],
2834 const double grad_w_dS[2])
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];
2848 double* mesh_trial_ref,
2849 double* mesh_grad_trial_ref,
2856 mapping.
calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y);
2864 double* mesh_trial_ref,
2880 double* mesh_trial_ref,
2881 double* mesh_grad_trial_ref,
2889 mapping.
calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y,
z);
2894 double* meshVelocity_dof,
2897 double* mesh_trial_ref,
2906 double* meshVelocity_dof,
2909 double* mesh_trial_ref,
2919 const int ebN_local,
2921 const int ebN_local_kb,
2924 double* mesh_trial_trace_ref,
2925 double* mesh_grad_trial_trace_ref,
2926 double* boundaryJac_ref,
2930 double* boundaryJac,
2931 double* metricTensor,
2932 double& metricTensorDetSqrt,
2944 mesh_trial_trace_ref,
2945 mesh_grad_trial_trace_ref,
2952 metricTensorDetSqrt,
2961 const int ebN_local,
2963 const int ebN_local_kb,
2966 double* mesh_trial_trace_ref,
2967 double* mesh_grad_trial_trace_ref,
2968 double* boundaryJac_ref,
2972 double* boundaryJac,
2973 double* metricTensor,
2974 double& metricTensorDetSqrt,
2987 mesh_trial_trace_ref,
2988 mesh_grad_trial_trace_ref,
2995 metricTensorDetSqrt,
3005 const int ebN_local,
3007 const int ebN_local_kb,
3008 double* mesh_velocity_dof,
3010 double* mesh_trial_trace_ref,
3014 double* boundaryJac,
3015 double* metricTensor,
3016 double& metricTensorDetSqrt)
3024 mesh_trial_trace_ref,
3030 metricTensorDetSqrt);
3034 const int ebN_local,
3036 const int ebN_local_kb,
3037 double* mesh_velocity_dof,
3039 double* mesh_trial_trace_ref,
3044 double* boundaryJac,
3045 double* metricTensor,
3046 double& metricTensorDetSqrt)
3054 mesh_trial_trace_ref,
3061 metricTensorDetSqrt);
3065 return stress[
sXX]*grad_test_dV[
X] + stress[
sXY]*grad_test_dV[
Y];
3081 return stress[
sYX]*grad_test_dV[
X] + stress[
sYY]*grad_test_dV[
Y];
3097 return stressFlux*disp_test_dS;
3101 return dstressFlux*disp_test_dS;
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>
3125 inline void calculateG(
double* jacInv,
double* G,
double& G_dd_G,
double& tr_G)
3127 for (
int I=0;I<1;I++)
3128 for (
int J=0;J<1;J++)
3131 for (
int K=0;K<1;K++)
3132 G[I+J] += jacInv[K+I]*jacInv[K+J];
3136 for (
int I=0;I<1;I++)
3139 for (
int J=0;J<1;J++)
3141 G_dd_G += G[I+J]*G[I+J];
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);
3153 inline void valFromDOF(
const double* dof,
const int* l2g_element,
const double* trial_ref,
double& val)
3156 for (
int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3157 val+=dof[l2g_element[j]]*trial_ref[j];
3160 inline void gradFromDOF(
const double* dof,
const int* l2g_element,
const double* grad_trial,
double* grad)
3162 for(
int I=0;I<1;I++)
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];
3169 inline void hessFromDOF(
const double* dof,
const int* l2g_element,
const double* hess_trial,
double* hess)
3171 for(
int I=0;I<1;I++)
3172 for(
int J=0;J<1;J++)
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];
3183 for (
int j=0;j<NDOF_TRIAL_ELEMENT;j++)
3184 val+=dof[j]*trial_ref[j];
3189 for(
int I=0;I<1;I++)
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];
3196 inline void gradTrialFromRef(
const double* grad_trial_ref,
const double* jacInv,
double* grad_trial)
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];
3218 inline void DOFaverage (
const double* dof,
const int* l2g_element,
double& val)
3222 for (
int j=0; j<NDOF_MESH_TRIAL_ELEMENT; j++)
3223 val+=dof[l2g_element[j]];
3225 val /= NDOF_MESH_TRIAL_ELEMENT;
3229 inline void hessTrialFromRef(
const double* hess_trial_ref,
const double* jacInv,
double* hess_trial)
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];
3243 inline void gradTestFromRef(
const double* grad_test_ref,
const double* jacInv,
double* grad_test)
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];
3254 inline void backwardEuler(
const double& dt,
const double& m_old,
const double& m,
const double& dm,
double& mt,
double& dmt)
3260 inline void bdf(
const double& alpha,
const double& beta,
const double& m,
const double& dm,
double& mt,
double& dmt)
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)
3273 inline double Mass_weak(
const double& mt,
const double& w_dV)
3322 const double& p_avg,
3326 if (viscosity==0.){
return 0.;}
3327 return (1./viscosity)*(p-p_avg)*(
q-1./3.)*dV;
3331 const double grad_w_dV[1])
3334 for(
int I=0;I<1;I++)
3335 tmp -=
f[I]*grad_w_dV[I];
3341 const double grad_w_dV[1])
3344 for(
int I=0;I<1;I++)
3345 tmp -=
df[I]*
v*grad_w_dV[I];
3350 const double grad_u[1])
3353 for(
int I=0;I<1;I++)
3354 tmp +=
df[I]*grad_u[I];
3359 const double grad_v[1])
3362 for(
int I=0;I<1;I++)
3363 tmp +=
df[I]*grad_v[I];
3368 const double grad_w_dV[1])
3371 for(
int I=0;I<1;I++)
3372 tmp -=
df[I]*grad_w_dV[I];
3383 const double grad_v[1],
3387 for(
int I=0;I<1;I++)
3388 tmp += dH[I]*grad_v[I]*w_dV;
3393 const double grad_u[1])
3396 for(
int I=0;I<1;I++)
3397 tmp += dH[I]*grad_u[I];
3402 const double grad_v[1])
3405 for(
int I=0;I<1;I++)
3406 tmp += dH[I]*grad_v[I];
3411 const double grad_w_dV[1])
3414 for(
int I=0;I<1;I++)
3415 tmp -= dH[I]*grad_w_dV[I];
3422 const double grad_phi[1],
3423 const double grad_w_dV[1])
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];
3436 const double grad_phi[1],
3437 const double grad_w_dV[1],
3440 const double grad_v[1])
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++)
3446 daProduct += da[m]*grad_phi[colind[m]]*grad_w_dV[I];
3447 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
3449 return daProduct*
v+dphiProduct*dphi;
3455 const double grad_v[1],
3456 const double grad_w_dV[1])
3458 double dphiProduct=0.0;
3459 for (
int I=0;I<1;I++)
3460 for (
int m=rowptr[I];m<rowptr[I+1];m++)
3462 dphiProduct += a[m]*grad_v[colind[m]]*grad_w_dV[I];
3498 const double& elementDiameter,
3499 const double& strong_residual,
3500 const double grad_u[1],
3507 h = elementDiameter;
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);
3519 const double G[1*1],
3520 const double& strong_residual,
3521 const double grad_u[1],
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];
3529 numDiff = shockCapturingDiffusion*fabs(strong_residual)/(sqrt(den+1.0e-12));
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],
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];
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);
3553 const double G[1*1],
3554 const double& strong_residual,
3555 const double vel[1],
3556 const double grad_u[1],
3559 double den1 = 0.0,den2=0.0, nom=0.0;
3560 for (
int I=0;I<1;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];
3567 numDiff = shockCapturingDiffusion*fabs(strong_residual)*(sqrt(nom/(den1*den2 + 1.0e-12)));
3572 const double& elementDiameter,
3573 const double& strong_residual,
3574 const double grad_u[1],
3576 double& gradNorm_last,
3582 h = elementDiameter;
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);
3589 numDiff =
num/gradNorm_last;
3593 const double& Lstar_w_dV)
3595 return error*Lstar_w_dV;
3599 const double& Lstar_w_dV)
3601 return derror*Lstar_w_dV;
3605 const double grad_u[1],
3606 const double grad_w_dV[1])
3609 for (
int I=0;I<1;I++)
3610 tmp += numDiff*grad_u[I]*grad_w_dV[I];
3615 const double grad_v[1],
3616 const double grad_w_dV[1])
3619 for (
int I=0;I<1;I++)
3620 tmp += numDiff*grad_v[I]*grad_w_dV[I];
3641 return dflux_left*
v;
3647 return dflux_left*
v;
3651 const int& isFluxBoundary,
3652 const double& sigma,
3655 const double normal[1],
3657 const double grad_w_dS[1])
3660 for(
int I=0;I<1;I++)
3662 tmp += normal[I]*grad_w_dS[I];
3664 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*(
u-bc_u)*a;
3669 const int& isFluxBoundary,
3670 const double& sigma,
3672 const double normal[1],
3674 const double grad_w_dS[1])
3677 for(
int I=0;I<1;I++)
3679 tmp += normal[I]*grad_w_dS[I];
3681 tmp *= (1.0-isFluxBoundary)*isDOFBoundary*sigma*
v*a;
3686 const int& isFluxBoundary,
3687 const double& sigma,
3690 const double normal[1],
3694 const double grad_w_dS[1])
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];
3704 const int& isFluxBoundary,
3705 const double& sigma,
3707 const double normal[1],
3711 const double grad_w_dS[1])
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];
3725 double* mesh_trial_ref,
3726 double* mesh_grad_trial_ref,
3733 mapping.
calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y);
3741 double* mesh_trial_ref,
3757 double* mesh_trial_ref,
3758 double* mesh_grad_trial_ref,
3766 mapping.
calculateMapping_element(eN,k,mesh_dof,mesh_l2g,mesh_trial_ref,mesh_grad_trial_ref,jac,jacDet,jacInv,x,y,
z);
3771 double* meshVelocity_dof,
3774 double* mesh_trial_ref,
3783 double* meshVelocity_dof,
3786 double* mesh_trial_ref,
3796 const int ebN_local,
3798 const int ebN_local_kb,
3801 double* mesh_trial_trace_ref,
3802 double* mesh_grad_trial_trace_ref,
3803 double* boundaryJac_ref,
3807 double* boundaryJac,
3808 double* metricTensor,
3809 double& metricTensorDetSqrt,
3821 mesh_trial_trace_ref,
3822 mesh_grad_trial_trace_ref,
3829 metricTensorDetSqrt,
3838 const int ebN_local,
3840 const int ebN_local_kb,
3843 double* mesh_trial_trace_ref,
3844 double* mesh_grad_trial_trace_ref,
3845 double* boundaryJac_ref,
3849 double* boundaryJac,
3850 double* metricTensor,
3851 double& metricTensorDetSqrt,
3864 mesh_trial_trace_ref,
3865 mesh_grad_trial_trace_ref,
3872 metricTensorDetSqrt,
3882 const int ebN_local,
3884 const int ebN_local_kb,
3885 double* mesh_velocity_dof,
3887 double* mesh_trial_trace_ref,
3891 double* boundaryJac,
3892 double* metricTensor,
3893 double& metricTensorDetSqrt)
3901 mesh_trial_trace_ref,
3907 metricTensorDetSqrt);
3911 const int ebN_local,
3913 const int ebN_local_kb,
3914 double* mesh_velocity_dof,
3916 double* mesh_trial_trace_ref,
3921 double* boundaryJac,
3922 double* metricTensor,
3923 double& metricTensorDetSqrt)
3931 mesh_trial_trace_ref,
3938 metricTensorDetSqrt);
3942 return stress[
sXX]*grad_test_dV[
X];
3952 return stressFlux*disp_test_dS;
3956 return dstressFlux*disp_test_dS;
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 Mass_strong(const double &mt)
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])
double Stress_u_weak(double *stress, double *grad_test_dV)
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 ReactionJacobian_strong(const double &dr, const double &v)
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)
double Hamiltonian_weak(const double &H, const double &w_dV)
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 Mass_adjoint(const double &dmt, const double &w_dV)
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)
void calculateGScale(double *G, double *v, double &h)
double pressureProjection_weak(const double &viscosity, const double &p, const double &p_avg, const double &q, const double &dV)
double MassJacobian_strong(const double &dmt, const double &v)
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])
double Reaction_adjoint(const double &dr, const double &w_dV)
double Reaction_weak(const double &r, const double &w_dV)
double Mass_weak(const double &mt, const double &w_dV)
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 Reaction_strong(const double &r)
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 Mass_weak(const double &mt, const double &w_dV)
double NumericalDiffusion(const double &numDiff, const double grad_u[2], const double grad_w_dV[2])
double Mass_strong(const double &mt)
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)
double Stress_v_weak(double *stress, double *grad_test_dV)
void bdf(const double &alpha, const double &beta, const double &m, const double &dm, double &mt, double &dmt)
double ReactionJacobian_strong(const double &dr, const double &v)
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 Reaction_strong(const double &r)
void calculateGScale(double *G, double *v, double &h)
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 Stress_u_weak(double *stress, double *grad_test_dV)
double ReactionJacobian_weak(const double &dr, const double &v, const double &w_dV)
double Reaction_adjoint(const double &dr, const double &w_dV)
double AdvectionJacobian_weak(const double df[2], const double &v, const double grad_w_dV[2])
double Hamiltonian_weak(const double &H, const double &w_dV)
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 Reaction_weak(const double &r, const double &w_dV)
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 MassJacobian_strong(const double &dmt, const double &v)
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 Mass_adjoint(const double &dmt, const double &w_dV)
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)
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)
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)
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)
void gradTestFromRef(const double *grad_test_ref, const double *jacInv, double *grad_test)
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)
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])
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)
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)
double ReactionJacobian_strong(const double &dr, const double &v)
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)
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])
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])
void gradFromElementDOF(const double *dof, 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)
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)
double Hamiltonian_weak(const double &H, const double &w_dV)
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])
void bdfC2(const double &alpha, const double &beta, const double &m, const double &dm, const double &dm2, double &mt, double &dmt, double &dm2t)
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)
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)
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)
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)
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 calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
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)
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)
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt)
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)
void valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
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, double &yt, double &zt)
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)
void calculateMappingVelocity_element(const int eN, const int k, double *mesh_velocity_dof, int *mesh_l2g, double *mesh_trial_ref, double &xt, double &yt)
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 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 valFromDOF(const double *dof, const int *l2g_element, const double *trial_ref, double &val)
void gradFromDOF(const double *dof, const int *l2g_element, const double *grad_trial, double *grad)
void calculateH_element(const int eN, const int k, double *h_dof, int *mesh_l2g, double *mesh_trial_ref, double &h)
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 hessFromDOF(const double *dof, const int *l2g_element, const double *hess_trial, double *hess)
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_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)
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)
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)
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)
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)