4 #define TR_ALPHA_EXT 1.0
16 int nInteriorElementBoundaries_global,
17 int nElementBoundaries_element,
18 int nQuadraturePoints_elementBoundary,
20 int* interiorElementBoundaries,
21 int* elementBoundaryElements,
22 int* elementBoundaryLocalElementBoundaries,
31 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,J;
32 double left_flux,right_flux,u_left,u_right,left_speed,right_speed,sonicSpeed;
34 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
36 ebN = interiorElementBoundaries[ebNI];
37 left_eN_global = elementBoundaryElements[ebN*2+0];
38 right_eN_global = elementBoundaryElements[ebN*2+1];
39 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
40 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
41 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
51 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
52 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
56 df[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
57 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
62 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
63 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
67 df[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
68 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
73 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
74 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
78 f[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
79 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
84 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
85 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
89 f[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
90 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
95 u_left =
u[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
96 left_ebN_element*nQuadraturePoints_elementBoundary+
98 u_right=
u[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
99 right_ebN_element*nQuadraturePoints_elementBoundary+
108 if (u_left >= u_right)
110 if (left_flux >= right_flux)
112 flux[ebN*nQuadraturePoints_elementBoundary+
114 dflux_left[ebN*nQuadraturePoints_elementBoundary+
116 dflux_right[ebN*nQuadraturePoints_elementBoundary+
121 flux[ebN*nQuadraturePoints_elementBoundary+
123 dflux_left[ebN*nQuadraturePoints_elementBoundary+
125 dflux_right[ebN*nQuadraturePoints_elementBoundary+
132 flux[ebN*nQuadraturePoints_elementBoundary+
134 dflux_left[ebN*nQuadraturePoints_elementBoundary+
136 dflux_right[ebN*nQuadraturePoints_elementBoundary+
139 if (right_flux <= flux[ebN*nQuadraturePoints_elementBoundary+k])
141 flux[ebN*nQuadraturePoints_elementBoundary+
143 dflux_left[ebN*nQuadraturePoints_elementBoundary+
145 dflux_right[ebN*nQuadraturePoints_elementBoundary+
148 if (u_left <= sonicPoint && sonicPoint <= u_right &&
149 sonicFlux < flux[ebN*nQuadraturePoints_elementBoundary+k])
151 flux[ebN*nQuadraturePoints_elementBoundary+
154 dflux_left[ebN*nQuadraturePoints_elementBoundary+
156 dflux_right[ebN*nQuadraturePoints_elementBoundary+
180 int nInteriorElementBoundaries_global,
181 int nElementBoundaries_element,
182 int nQuadraturePoints_elementBoundary,
183 int nQuadraturePoints_element,
185 int* interiorElementBoundaries,
186 int* elementBoundaryElements,
187 int* elementBoundaryLocalElementBoundaries,
197 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,J;
198 double left_flux,right_flux,u_left,u_right,left_speed,right_speed,
199 maxSpeed_left,maxSpeed_right,maxSpeed_element,maxSpeed,tmp_left,tmp_right;
202 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
204 ebN = interiorElementBoundaries[ebNI];
205 left_eN_global = elementBoundaryElements[ebN*2+0];
206 right_eN_global = elementBoundaryElements[ebN*2+1];
207 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
208 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
213 maxSpeed_left =0.0; maxSpeed_right=0.0; maxSpeed=0.0;
214 for (k=0; k < nQuadraturePoints_element; k++)
216 tmp_left = 0.0; tmp_right=0.0;
218 for (J=0; J < nSpace; J++)
221 df_element[left_eN_global*nQuadraturePoints_element*nSpace+ k*nSpace + J]
223 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
224 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
229 df_element[right_eN_global*nQuadraturePoints_element*nSpace+ k*nSpace + J]
231 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
232 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
236 maxSpeed_left = fabs(tmp_left) > maxSpeed_left ? fabs(tmp_left) : maxSpeed_left;
237 maxSpeed_right= fabs(tmp_right) > maxSpeed_right ? fabs(tmp_right) : maxSpeed_right;
239 maxSpeed_element = maxSpeed_right > maxSpeed_left ? maxSpeed_right : maxSpeed_left;
240 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
246 for(J=0;J<nSpace;J++)
250 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
251 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
255 df[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
256 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
261 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
262 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
266 df[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
267 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
272 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
273 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
277 f[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
278 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
283 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
284 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
288 f[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
289 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
294 u_left =
u[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
295 left_ebN_element*nQuadraturePoints_elementBoundary+
297 u_right=
u[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
298 right_ebN_element*nQuadraturePoints_elementBoundary+
300 maxSpeed = fabs(left_speed) > maxSpeed_element ? fabs(left_speed) : maxSpeed_element;
301 maxSpeed = fabs(right_speed) > maxSpeed ? fabs(right_speed) : maxSpeed;
302 maxSpeed*= safetyFactor;
303 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.5*(left_flux + right_flux) - 0.5*maxSpeed*(u_right-u_left);
304 dflux_left[ebN*nQuadraturePoints_elementBoundary+k] = 0.5*left_speed + 0.5*maxSpeed;
305 dflux_right[ebN*nQuadraturePoints_elementBoundary+k]= 0.5*right_speed - 0.5*maxSpeed;
327 int nExteriorElementBoundaries_global,
328 int nElementBoundaries_element,
329 int nQuadraturePoints_elementBoundary,
330 int nQuadraturePoints_element,
332 int* exteriorElementBoundaries,
333 int* elementBoundaryElements,
334 int* elementBoundaryLocalElementBoundaries,
348 int ebNE,ebN,eN_global,ebN_element,k,J;
349 double left_flux,right_flux,u_left,u_right,left_speed,right_speed,
350 maxSpeed_element,maxSpeed,tmp_left;
353 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
355 ebN = exteriorElementBoundaries[ebNE];
356 eN_global = elementBoundaryElements[ebN*2+0];
357 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
362 maxSpeed_element =0.0; maxSpeed=0.0;
363 for (k=0; k < nQuadraturePoints_element; k++)
367 for (J=0; J < nSpace; J++)
370 df_element[eN_global*nQuadraturePoints_element*nSpace+ k*nSpace + J]
372 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
373 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
377 maxSpeed_element = fabs(tmp_left) > maxSpeed_element ? fabs(tmp_left) : maxSpeed_element;
379 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
385 for(J=0;J<nSpace;J++)
389 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
390 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
394 df[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
395 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
400 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
401 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
405 bc_df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
410 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
411 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
415 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
416 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
421 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
422 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
426 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
431 u_left =
u[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
432 ebN_element*nQuadraturePoints_elementBoundary+
434 u_right= bc_u[ebNE*nQuadraturePoints_elementBoundary+
436 maxSpeed = fabs(left_speed) > maxSpeed_element ? fabs(left_speed) : maxSpeed_element;
437 maxSpeed = fabs(right_speed) > maxSpeed ? fabs(right_speed) : maxSpeed;
438 maxSpeed*= safetyFactor;
439 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.5*(left_flux + right_flux) - 0.5*maxSpeed*(u_right-u_left);
440 dflux[ebN*nQuadraturePoints_elementBoundary+k] = 0.5*left_speed + 0.5*maxSpeed;
464 int nExteriorElementBoundaries_global,
465 int nQuadraturePoints_elementBoundary,
466 int nQuadraturePoints_element,
468 int* exteriorElementBoundaries,
469 int* elementBoundaryElements,
470 int* elementBoundaryLocalElementBoundaries,
484 int ebNE,ebN,eN_global,k,J;
485 double left_flux,right_flux,u_left,u_right,left_speed,right_speed,
486 maxSpeed_element,maxSpeed,tmp_left;
489 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
491 ebN = exteriorElementBoundaries[ebNE];
492 eN_global = elementBoundaryElements[ebN*2+0];
497 maxSpeed_element =0.0; maxSpeed=0.0;
498 for (k=0; k < nQuadraturePoints_element; k++)
502 for (J=0; J < nSpace; J++)
505 df_element[eN_global*nQuadraturePoints_element*nSpace+ k*nSpace + J]
507 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
511 maxSpeed_element = fabs(tmp_left) > maxSpeed_element ? fabs(tmp_left) : maxSpeed_element;
513 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
519 for(J=0;J<nSpace;J++)
523 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
527 df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
532 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
536 bc_df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
541 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
545 f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
550 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
554 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
559 u_left =
u[ebNE*nQuadraturePoints_elementBoundary+
561 u_right= bc_u[ebNE*nQuadraturePoints_elementBoundary+
563 maxSpeed = fabs(left_speed) > maxSpeed_element ? fabs(left_speed) : maxSpeed_element;
564 maxSpeed = fabs(right_speed) > maxSpeed ? fabs(right_speed) : maxSpeed;
565 maxSpeed*= safetyFactor;
566 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.5*(left_flux + right_flux) - 0.5*maxSpeed*(u_right-u_left);
567 dflux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.5*left_speed + 0.5*maxSpeed;
592 int nInteriorElementBoundaries_global,
593 int nElementBoundaries_element,
594 int nQuadraturePoints_elementBoundary,
595 int nQuadraturePoints_element,
597 int* interiorElementBoundaries,
598 int* elementBoundaryElements,
599 int* elementBoundaryLocalElementBoundaries,
603 double* lambda_bar_element,
606 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,J;
607 double left_flux,right_flux,u_left,u_right,
608 maxSpeed_left,maxSpeed_right,maxSpeed_element,maxSpeed,tmp_left,tmp_right;
611 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
613 ebN = interiorElementBoundaries[ebNI];
614 left_eN_global = elementBoundaryElements[ebN*2+0];
615 right_eN_global = elementBoundaryElements[ebN*2+1];
616 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
617 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
622 maxSpeed_left =0.0; maxSpeed_right=0.0; maxSpeed=0.0;
623 for (k=0; k < nQuadraturePoints_element; k++)
625 tmp_left = lambda_bar_element[left_eN_global*nQuadraturePoints_element + k];
626 tmp_right= lambda_bar_element[right_eN_global*nQuadraturePoints_element + k];
627 maxSpeed_left = fabs(tmp_left) > maxSpeed_left ? fabs(tmp_left) : maxSpeed_left;
628 maxSpeed_right= fabs(tmp_right) > maxSpeed_right ? fabs(tmp_right) : maxSpeed_right;
630 maxSpeed_element = maxSpeed_right > maxSpeed_left ? maxSpeed_right : maxSpeed_left;
631 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
635 for(J=0;J<nSpace;J++)
639 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
640 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
644 f[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
645 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
650 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
651 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
655 f[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
656 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
661 u_left =
u[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
662 left_ebN_element*nQuadraturePoints_elementBoundary+
664 u_right=
u[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
665 right_ebN_element*nQuadraturePoints_elementBoundary+
669 maxSpeed = maxSpeed_element;
670 maxSpeed*= safetyFactor;
671 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.5*(left_flux + right_flux) - 0.5*maxSpeed*(u_right-u_left);
693 int nExteriorElementBoundaries_global,
694 int nQuadraturePoints_elementBoundary,
695 int nQuadraturePoints_element,
697 int* exteriorElementBoundaries,
698 int* elementBoundaryElements,
699 int* elementBoundaryLocalElementBoundaries,
710 int ebNE,ebN,eN_global,k,J;
711 double left_flux,right_flux,u_left,u_right,
712 maxSpeed_element,maxSpeed,tmp_left;
715 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
717 ebN = exteriorElementBoundaries[ebNE];
718 eN_global = elementBoundaryElements[ebN*2+0];
723 maxSpeed_element =0.0; maxSpeed=0.0;
724 for (k=0; k < nQuadraturePoints_element; k++)
726 tmp_left = lambda_bar[eN_global*nQuadraturePoints_element + k];
727 maxSpeed_element = fabs(tmp_left) > maxSpeed_element ? fabs(tmp_left) : maxSpeed_element;
729 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
733 for(J=0;J<nSpace;J++)
737 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
741 f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
746 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
750 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
755 u_left =
u[ebNE*nQuadraturePoints_elementBoundary+
757 u_right= bc_u[ebNE*nQuadraturePoints_elementBoundary+
761 maxSpeed =maxSpeed_element;
762 maxSpeed*= safetyFactor;
763 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.5*(left_flux + right_flux) - 0.5*maxSpeed*(u_right-u_left);
785 int nExteriorElementBoundaries_global,
786 int nElementBoundaries_element,
787 int nQuadraturePoints_elementBoundary,
788 int nQuadraturePoints_element,
790 int* exteriorElementBoundaries,
791 int* elementBoundaryElements,
792 int* elementBoundaryLocalElementBoundaries,
803 int ebNE,ebN,eN_global,ebN_element,k,J;
804 double left_flux,right_flux,u_left,u_right,
805 maxSpeed_element,maxSpeed,tmp_left;
808 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
810 ebN = exteriorElementBoundaries[ebNE];
811 eN_global = elementBoundaryElements[ebN*2+0];
812 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
817 maxSpeed_element =0.0; maxSpeed=0.0;
818 for (k=0; k < nQuadraturePoints_element; k++)
820 tmp_left = lambda_bar[eN_global*nQuadraturePoints_element+k];
821 maxSpeed_element = fabs(tmp_left) > maxSpeed_element ? fabs(tmp_left) : maxSpeed_element;
823 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
827 for(J=0;J<nSpace;J++)
831 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
832 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
836 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
837 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
842 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
843 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
847 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
852 u_left =
u[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
853 ebN_element*nQuadraturePoints_elementBoundary+
855 u_right= bc_u[ebNE*nQuadraturePoints_elementBoundary+
859 maxSpeed = maxSpeed_element;
860 maxSpeed*= safetyFactor;
861 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.5*(left_flux + right_flux) - 0.5*maxSpeed*(u_right-u_left);
892 double penalty_floor,
893 int nInteriorElementBoundaries_global,
894 int nElementBoundaries_element,
895 int nQuadraturePoints_elementBoundary,
897 int* interiorElementBoundaries,
898 int* elementBoundaryElements,
899 int* elementBoundaryLocalElementBoundaries,
907 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,J,I,nSpace2=nSpace*nSpace;
908 double diffusiveVelocityComponent_I,max_a=0.0;
909 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
911 ebN = interiorElementBoundaries[ebNI];
912 left_eN_global = elementBoundaryElements[ebN*2+0];
913 right_eN_global = elementBoundaryElements[ebN*2+1];
914 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
915 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
916 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
919 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
920 for(I=0;I<nSpace;I++)
922 diffusiveVelocityComponent_I=0.0;
923 for(J=0;J<nSpace;J++)
925 diffusiveVelocityComponent_I
927 (a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
928 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
933 grad_phi[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
934 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
937 a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
938 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
943 grad_phi[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
944 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
946 max_a = fmax(max_a,a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
947 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
951 max_a = fmax(max_a,a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
952 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
957 flux[ebN*nQuadraturePoints_elementBoundary+
960 diffusiveVelocityComponent_I
962 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
963 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
967 flux[ebN*nQuadraturePoints_elementBoundary+
970 max_a = fmax(max_a,penalty_floor);
971 double penalty_flux = penalty[ebN*nQuadraturePoints_elementBoundary+
974 (
u[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
975 left_ebN_element*nQuadraturePoints_elementBoundary+
977 u[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
978 right_ebN_element*nQuadraturePoints_elementBoundary+
980 if (scale_penalty) penalty_flux *= max_a;
981 flux[ebN*nQuadraturePoints_elementBoundary+
988 double penalty_floor,
989 int nInteriorElementBoundaries_global,
990 int nElementBoundaries_element,
991 int nQuadraturePoints_elementBoundary,
995 int* interiorElementBoundaries,
996 int* elementBoundaryElements,
997 int* elementBoundaryLocalElementBoundaries,
1005 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,I,m,
nnz=rowptr[nSpace];
1006 double diffusiveVelocityComponent_I,max_a=0.0;
1007 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
1009 ebN = interiorElementBoundaries[ebNI];
1010 left_eN_global = elementBoundaryElements[ebN*2+0];
1011 right_eN_global = elementBoundaryElements[ebN*2+1];
1012 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
1013 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
1014 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
1017 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
1018 for(I=0;I<nSpace;I++)
1020 diffusiveVelocityComponent_I=0.0;
1021 for(m=rowptr[I];m<rowptr[I+1];m++)
1023 diffusiveVelocityComponent_I
1025 (a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1026 left_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1030 grad_phi[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1031 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1034 a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1035 right_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1039 grad_phi[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1040 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1041 k*nSpace+colind[m]]);
1042 max_a = fmax(max_a,a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1043 left_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1046 max_a = fmax(max_a,a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1047 right_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1051 flux[ebN*nQuadraturePoints_elementBoundary+
1054 diffusiveVelocityComponent_I
1056 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1057 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1061 flux[ebN*nQuadraturePoints_elementBoundary+
1064 max_a = fmax(max_a,penalty_floor);
1065 double penalty_flux = penalty[ebN*nQuadraturePoints_elementBoundary+
1068 (
u[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
1069 left_ebN_element*nQuadraturePoints_elementBoundary+
1071 u[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
1072 right_ebN_element*nQuadraturePoints_elementBoundary+
1074 if (scale_penalty) penalty_flux *= max_a;
1078 flux[ebN*nQuadraturePoints_elementBoundary+
1256 double penalty_floor,
1257 int nInteriorElementBoundaries_global,
1258 int nElementBoundaries_element,
1259 int nQuadraturePoints_elementBoundary,
1260 int nDOF_trial_element,
1263 int* interiorElementBoundaries,
1264 int* elementBoundaryElements,
1265 int* elementBoundaryLocalElementBoundaries,
1274 double* fluxJacobian)
1276 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,j,left_j_global,right_j_global,I,J,nSpace2=nSpace*nSpace;
1277 double leftJacobian,rightJacobian,diffusiveVelocityComponent_I_leftJacobian,diffusiveVelocityComponent_I_rightJacobian,diffusiveVelocityComponent_I_leftJacobian2,diffusiveVelocityComponent_I_rightJacobian2,max_a=0.0;
1278 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
1280 ebN = interiorElementBoundaries[ebNI];
1281 left_eN_global = elementBoundaryElements[ebN*2+0];
1282 right_eN_global = elementBoundaryElements[ebN*2+1];
1283 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
1284 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
1285 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
1288 for(j=0;j<nDOF_trial_element;j++)
1292 left_j_global = l2g[left_eN_global*nDOF_trial_element+j];
1293 right_j_global= l2g[right_eN_global*nDOF_trial_element+j];
1294 for(I=0;I<nSpace;I++)
1296 diffusiveVelocityComponent_I_leftJacobian=0.0;
1297 diffusiveVelocityComponent_I_leftJacobian2=0.0;
1298 diffusiveVelocityComponent_I_rightJacobian=0.0;
1299 diffusiveVelocityComponent_I_rightJacobian2=0.0;
1300 for(J=0;J<nSpace;J++)
1302 diffusiveVelocityComponent_I_leftJacobian
1304 da[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1305 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1310 grad_phi[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1311 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1314 diffusiveVelocityComponent_I_rightJacobian
1316 da[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1317 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1322 grad_phi[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1323 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1326 diffusiveVelocityComponent_I_leftJacobian2
1328 a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1329 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1334 grad_v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1335 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1336 k*nDOF_trial_element*nSpace+
1339 diffusiveVelocityComponent_I_rightJacobian2
1341 a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1342 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1347 grad_v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1348 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1349 k*nDOF_trial_element*nSpace+
1352 max_a = fmax(max_a,a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1353 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1357 max_a = fmax(max_a,a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1358 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1366 (diffusiveVelocityComponent_I_leftJacobian
1368 v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1369 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1370 k*nDOF_trial_element+
1373 diffusiveVelocityComponent_I_leftJacobian2*
1374 dphi[left_j_global])
1376 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1377 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1382 (diffusiveVelocityComponent_I_rightJacobian
1384 v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1385 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1386 k*nDOF_trial_element+
1389 diffusiveVelocityComponent_I_rightJacobian2*dphi[right_j_global])
1391 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1392 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1396 leftJacobian *= 0.5;
1397 rightJacobian *= 0.5;
1398 max_a = fmax(max_a,penalty_floor);
1399 double penaltyJacobian_term = penalty[ebN*nQuadraturePoints_elementBoundary+
1401 if (scale_penalty) penaltyJacobian_term *= max_a;
1404 penaltyJacobian_term
1406 v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1407 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1408 k*nDOF_trial_element+
1412 penaltyJacobian_term
1414 v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1415 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1416 k*nDOF_trial_element+
1418 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1419 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1420 k*nDOF_trial_element+
1422 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1423 1*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1424 k*nDOF_trial_element+
1425 j] += rightJacobian;
1432 double penalty_floor,
1433 int nInteriorElementBoundaries_global,
1434 int nElementBoundaries_element,
1435 int nQuadraturePoints_elementBoundary,
1436 int nDOF_trial_element,
1441 int* interiorElementBoundaries,
1442 int* elementBoundaryElements,
1443 int* elementBoundaryLocalElementBoundaries,
1452 double* fluxJacobian)
1454 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,j,left_j_global,right_j_global,I,m,
nnz=rowptr[nSpace];
1455 double leftJacobian,rightJacobian,diffusiveVelocityComponent_I_leftJacobian,diffusiveVelocityComponent_I_rightJacobian,diffusiveVelocityComponent_I_leftJacobian2,diffusiveVelocityComponent_I_rightJacobian2,max_a;
1456 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
1458 ebN = interiorElementBoundaries[ebNI];
1459 left_eN_global = elementBoundaryElements[ebN*2+0];
1460 right_eN_global = elementBoundaryElements[ebN*2+1];
1461 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
1462 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
1463 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
1466 for(j=0;j<nDOF_trial_element;j++)
1470 left_j_global = l2g[left_eN_global*nDOF_trial_element+j];
1471 right_j_global= l2g[right_eN_global*nDOF_trial_element+j];
1472 for(I=0;I<nSpace;I++)
1474 diffusiveVelocityComponent_I_leftJacobian=0.0;
1475 diffusiveVelocityComponent_I_leftJacobian2=0.0;
1476 diffusiveVelocityComponent_I_rightJacobian=0.0;
1477 diffusiveVelocityComponent_I_rightJacobian2=0.0;
1478 for(m=rowptr[I];m<rowptr[I+1];m++)
1480 diffusiveVelocityComponent_I_leftJacobian
1482 da[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1483 left_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1487 grad_phi[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1488 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1491 diffusiveVelocityComponent_I_rightJacobian
1493 da[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1494 right_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1498 grad_phi[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1499 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1502 diffusiveVelocityComponent_I_leftJacobian2
1504 a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1505 left_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1509 grad_v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1510 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1511 k*nDOF_trial_element*nSpace+
1514 diffusiveVelocityComponent_I_rightJacobian2
1516 a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1517 right_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1521 grad_v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1522 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
1523 k*nDOF_trial_element*nSpace+
1526 max_a = fmax(max_a,a[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1527 left_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1530 max_a = fmax(max_a,a[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
1531 right_ebN_element*nQuadraturePoints_elementBoundary*
nnz+
1538 (diffusiveVelocityComponent_I_leftJacobian
1540 v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1541 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1542 k*nDOF_trial_element+
1545 diffusiveVelocityComponent_I_leftJacobian2*
1546 dphi[left_j_global])
1548 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1549 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1554 (diffusiveVelocityComponent_I_rightJacobian
1556 v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1557 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1558 k*nDOF_trial_element+
1561 diffusiveVelocityComponent_I_rightJacobian2*dphi[right_j_global])
1563 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1564 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1568 leftJacobian *= 0.5;
1569 rightJacobian *= 0.5;
1570 max_a = fmax(max_a,penalty_floor);
1571 double penaltyJacobian_term = penalty[ebN*nQuadraturePoints_elementBoundary+
1573 if (scale_penalty) penaltyJacobian_term *= max_a;
1576 penaltyJacobian_term
1578 v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1579 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1580 k*nDOF_trial_element+
1584 penaltyJacobian_term
1586 v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1587 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1588 k*nDOF_trial_element+
1590 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1591 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1592 k*nDOF_trial_element+
1594 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1595 1*nQuadraturePoints_elementBoundary*nDOF_trial_element+
1596 k*nDOF_trial_element+
1597 j] += rightJacobian;
1931 double penalty_floor,
1932 int nExteriorElementBoundaries_global,
1933 int nElementBoundaries_element,
1934 int nQuadraturePoints_elementBoundary,
1936 int* exteriorElementBoundaries,
1937 int* elementBoundaryElements,
1938 int* elementBoundaryLocalElementBoundaries,
1942 double* bc_grad_phi,
1950 int ebNE,ebN,eN_global,ebN_element,k,J,I,nSpace2=nSpace*nSpace;
1951 double diffusiveVelocityComponent_I,penaltyFlux,max_a=0.0;
1952 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
1954 ebN = exteriorElementBoundaries[ebNE];
1955 eN_global = elementBoundaryElements[ebN*2+0];
1956 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
1957 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
1959 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
1962 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
1963 for(I=0;I<nSpace;I++)
1965 diffusiveVelocityComponent_I=0.0;
1966 for(J=0;J<nSpace;J++)
1968 diffusiveVelocityComponent_I
1970 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1971 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1976 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1977 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1979 max_a = fmax(max_a,a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
1980 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
1985 flux[ebN*nQuadraturePoints_elementBoundary+k]
1987 diffusiveVelocityComponent_I
1989 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
1990 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
1994 max_a = fmax(penalty_floor,max_a);
1995 penaltyFlux = penalty[ebN*nQuadraturePoints_elementBoundary+
1998 (
u[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
1999 ebN_element*nQuadraturePoints_elementBoundary+
2002 bc_u[ebNE*nQuadraturePoints_elementBoundary+
2004 if (scale_penalty) penaltyFlux *= max_a;
2005 flux[ebN*nQuadraturePoints_elementBoundary+k] += penaltyFlux;
2011 double penalty_floor,
2012 int nExteriorElementBoundaries_global,
2013 int nElementBoundaries_element,
2014 int nQuadraturePoints_elementBoundary,
2018 int* exteriorElementBoundaries,
2019 int* elementBoundaryElements,
2020 int* elementBoundaryLocalElementBoundaries,
2024 double* bc_grad_phi,
2032 int ebNE,ebN,eN_global,ebN_element,k,I,m,
nnz=rowptr[nSpace];
2033 double diffusiveVelocityComponent_I,penaltyFlux,max_a=0.0;
2034 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2036 ebN = exteriorElementBoundaries[ebNE];
2037 eN_global = elementBoundaryElements[ebN*2+0];
2038 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
2039 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2041 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
2042 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
2045 for(I=0;I<nSpace;I++)
2047 diffusiveVelocityComponent_I=0.0;
2048 for(m=rowptr[I];m<rowptr[I+1];m++)
2050 diffusiveVelocityComponent_I
2052 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
2053 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
2057 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2058 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2059 k*nSpace+colind[m]];
2060 max_a = fmax(max_a,a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
2061 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
2065 flux[ebN*nQuadraturePoints_elementBoundary+k]
2067 diffusiveVelocityComponent_I
2069 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2070 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2074 max_a = fmax(penalty_floor,max_a);
2075 penaltyFlux = penalty[ebN*nQuadraturePoints_elementBoundary+
2078 (
u[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
2079 ebN_element*nQuadraturePoints_elementBoundary+
2082 bc_u[ebNE*nQuadraturePoints_elementBoundary+
2084 if (scale_penalty) penaltyFlux *= max_a;
2085 flux[ebN*nQuadraturePoints_elementBoundary+k] += penaltyFlux;
2094 double penalty_floor,
2095 int nExteriorElementBoundaries_global,
2096 int nQuadraturePoints_elementBoundary,
2098 int* exteriorElementBoundaries,
2099 int* elementBoundaryElements,
2100 int* elementBoundaryLocalElementBoundaries,
2104 double* bc_grad_phi,
2112 int ebNE,k,J,I,nSpace2=nSpace*nSpace;
2113 double diffusiveVelocityComponent_I,penaltyFlux,max_a=0.0;
2114 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2116 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2118 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
2120 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
2122 for(I=0;I<nSpace;I++)
2124 diffusiveVelocityComponent_I=0.0;
2125 for(J=0;J<nSpace;J++)
2127 diffusiveVelocityComponent_I
2129 a[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
2134 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2136 max_a = fmax(max_a,a[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
2141 flux[ebNE*nQuadraturePoints_elementBoundary+k]
2143 diffusiveVelocityComponent_I
2145 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2149 max_a = fmax(max_a,penalty_floor);
2150 penaltyFlux = penalty[ebNE*nQuadraturePoints_elementBoundary+
2153 (
u[ebNE*nQuadraturePoints_elementBoundary+
2156 bc_u[ebNE*nQuadraturePoints_elementBoundary+
2158 if (scale_penalty) penaltyFlux *= max_a;
2159 flux[ebNE*nQuadraturePoints_elementBoundary+k] += penaltyFlux;
2166 double penalty_floor,
2167 int nExteriorElementBoundaries_global,
2168 int nQuadraturePoints_elementBoundary,
2172 int* exteriorElementBoundaries,
2173 int* elementBoundaryElements,
2174 int* elementBoundaryLocalElementBoundaries,
2178 double* bc_grad_phi,
2186 int ebNE,k,I,m,
nnz=rowptr[nSpace];
2187 double diffusiveVelocityComponent_I,penaltyFlux,max_a;
2188 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2190 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2192 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
2194 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
2196 for(I=0;I<nSpace;I++)
2198 diffusiveVelocityComponent_I=0.0;
2199 for(m=rowptr[I];m<rowptr[I+1];m++)
2201 diffusiveVelocityComponent_I
2203 a[ebNE*nQuadraturePoints_elementBoundary*
nnz+
2207 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2208 k*nSpace+colind[m]];
2209 max_a = fmax(max_a,a[ebNE*nQuadraturePoints_elementBoundary*
nnz+
2213 flux[ebNE*nQuadraturePoints_elementBoundary+k]
2215 diffusiveVelocityComponent_I
2217 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2221 max_a = fmax(penalty_floor,max_a);
2222 penaltyFlux = penalty[ebNE*nQuadraturePoints_elementBoundary+
2225 (
u[ebNE*nQuadraturePoints_elementBoundary+
2228 bc_u[ebNE*nQuadraturePoints_elementBoundary+
2231 if (scale_penalty) penaltyFlux *= max_a;
2232 flux[ebNE*nQuadraturePoints_elementBoundary+k] += penaltyFlux;
2524 int nElementBoundaries_element,
2525 int nQuadraturePoints_elementBoundary,
2527 int* exteriorElementBoundaries,
2528 int* elementBoundaryElements,
2529 int* elementBoundaryLocalElementBoundaries,
2533 double* bc_grad_phi,
2541 int ebNE,ebN,eN_global,ebN_element,k,J,I,nSpace2=nSpace*nSpace;
2542 double diffusiveVelocityComponent_I;
2543 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2545 ebN = exteriorElementBoundaries[ebNE];
2546 eN_global = elementBoundaryElements[ebN*2+0];
2547 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
2548 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2550 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
2551 for(I=0;I<nSpace;I++)
2553 diffusiveVelocityComponent_I=0.0;
2554 for(J=0;J<nSpace;J++)
2556 diffusiveVelocityComponent_I
2558 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
2559 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
2564 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2565 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2568 flux[ebN*nQuadraturePoints_elementBoundary+k]
2570 diffusiveVelocityComponent_I
2572 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2573 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2581 int nElementBoundaries_element,
2582 int nQuadraturePoints_elementBoundary,
2586 int* exteriorElementBoundaries,
2587 int* elementBoundaryElements,
2588 int* elementBoundaryLocalElementBoundaries,
2592 double* bc_grad_phi,
2600 int ebNE,ebN,eN_global,ebN_element,k,I,m,
nnz=rowptr[nSpace];
2601 double diffusiveVelocityComponent_I;
2602 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2604 ebN = exteriorElementBoundaries[ebNE];
2605 eN_global = elementBoundaryElements[ebN*2+0];
2606 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
2607 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2609 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
2610 for(I=0;I<nSpace;I++)
2612 diffusiveVelocityComponent_I=0.0;
2613 for(m=rowptr[I];m<rowptr[I+1];m++)
2615 diffusiveVelocityComponent_I
2617 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
2618 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
2622 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2623 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2624 k*nSpace+colind[m]];
2626 flux[ebN*nQuadraturePoints_elementBoundary+k]
2628 diffusiveVelocityComponent_I
2630 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2631 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2639 int nQuadraturePoints_elementBoundary,
2641 int* exteriorElementBoundaries,
2642 int* elementBoundaryElements,
2643 int* elementBoundaryLocalElementBoundaries,
2647 double* bc_grad_phi,
2655 int ebNE,k,J,I,nSpace2=nSpace*nSpace;
2656 double diffusiveVelocityComponent_I;
2657 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2659 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2661 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
2662 for(I=0;I<nSpace;I++)
2664 diffusiveVelocityComponent_I=0.0;
2665 for(J=0;J<nSpace;J++)
2667 diffusiveVelocityComponent_I
2669 a[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
2674 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2677 flux[ebNE*nQuadraturePoints_elementBoundary+k]
2679 diffusiveVelocityComponent_I
2681 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2689 int nQuadraturePoints_elementBoundary,
2693 int* exteriorElementBoundaries,
2694 int* elementBoundaryElements,
2695 int* elementBoundaryLocalElementBoundaries,
2699 double* bc_grad_phi,
2707 int ebNE,k,I,m,
nnz=rowptr[nSpace];
2708 double diffusiveVelocityComponent_I;
2709 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2711 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2713 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
2714 for(I=0;I<nSpace;I++)
2716 diffusiveVelocityComponent_I=0.0;
2717 for(m=rowptr[I];m<rowptr[I+1];m++)
2719 diffusiveVelocityComponent_I
2721 a[ebNE*nQuadraturePoints_elementBoundary*
nnz+
2725 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2726 k*nSpace+colind[m]];
2728 flux[ebNE*nQuadraturePoints_elementBoundary+k]
2730 diffusiveVelocityComponent_I
2732 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
2744 double penalty_floor,
2745 int nExteriorElementBoundaries_global,
2746 int nElementBoundaries_element,
2747 int nQuadraturePoints_elementBoundary,
2748 int nDOF_trial_element,
2751 int* exteriorElementBoundaries,
2752 int* elementBoundaryElements,
2753 int* elementBoundaryLocalElementBoundaries,
2763 double* fluxJacobian)
2765 int ebNE,ebN,eN_global,ebN_element,k,j,j_global,I,J,nSpace2=nSpace*nSpace;
2766 double Jacobian,diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2,max_a=0.0;
2767 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2769 ebN = exteriorElementBoundaries[ebNE];
2770 eN_global = elementBoundaryElements[ebN*2+0];
2771 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
2772 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2774 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
2777 for(j=0;j<nDOF_trial_element;j++)
2780 j_global = l2g[eN_global*nDOF_trial_element+j];
2781 for(I=0;I<nSpace;I++)
2783 diffusiveVelocityComponent_I_Jacobian=0.0;
2784 diffusiveVelocityComponent_I_Jacobian2=0.0;
2785 for(J=0;J<nSpace;J++)
2787 diffusiveVelocityComponent_I_Jacobian
2789 da[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
2790 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
2795 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2796 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2799 diffusiveVelocityComponent_I_Jacobian2
2801 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
2802 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
2807 grad_v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
2808 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
2809 k*nDOF_trial_element*nSpace+
2812 max_a = fmax(max_a,a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
2813 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
2820 (diffusiveVelocityComponent_I_Jacobian
2822 v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2823 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2824 k*nDOF_trial_element+
2827 diffusiveVelocityComponent_I_Jacobian2*
2830 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2831 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2835 max_a = fmax(penalty_floor,max_a);
2836 double penaltyJacobian = penalty[ebN*nQuadraturePoints_elementBoundary+
2839 v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2840 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2841 k*nDOF_trial_element+
2843 if (scale_penalty) penaltyJacobian *= max_a;
2844 Jacobian += penaltyJacobian;
2845 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2846 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2847 k*nDOF_trial_element+
2857 double penalty_floor,
2858 int nExteriorElementBoundaries_global,
2859 int nElementBoundaries_element,
2860 int nQuadraturePoints_elementBoundary,
2861 int nDOF_trial_element,
2866 int* exteriorElementBoundaries,
2867 int* elementBoundaryElements,
2868 int* elementBoundaryLocalElementBoundaries,
2878 double* fluxJacobian)
2880 int ebNE,ebN,eN_global,ebN_element,k,j,j_global,I,m,
nnz=rowptr[nSpace];
2881 double Jacobian,diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2,max_a=0.0;
2882 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2884 ebN = exteriorElementBoundaries[ebNE];
2885 eN_global = elementBoundaryElements[ebN*2+0];
2886 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
2887 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
2889 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
2892 for(j=0;j<nDOF_trial_element;j++)
2895 j_global = l2g[eN_global*nDOF_trial_element+j];
2896 for(I=0;I<nSpace;I++)
2898 diffusiveVelocityComponent_I_Jacobian=0.0;
2899 diffusiveVelocityComponent_I_Jacobian2=0.0;
2900 for(m=rowptr[I];m<rowptr[I+1];m++)
2902 diffusiveVelocityComponent_I_Jacobian
2904 da[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
2905 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
2909 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2910 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2913 diffusiveVelocityComponent_I_Jacobian2
2915 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
2916 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
2920 grad_v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
2921 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
2922 k*nDOF_trial_element*nSpace+
2925 max_a = fmax(max_a,a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
2926 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
2932 (diffusiveVelocityComponent_I_Jacobian
2934 v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2935 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2936 k*nDOF_trial_element+
2939 diffusiveVelocityComponent_I_Jacobian2*
2942 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
2943 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
2947 max_a = fmax(penalty_floor,max_a);
2948 double penaltyJacobian = penalty[ebN*nQuadraturePoints_elementBoundary+
2951 v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2952 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2953 k*nDOF_trial_element+
2956 if (scale_penalty) penaltyJacobian *= max_a;
2958 Jacobian += penaltyJacobian;
2960 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2961 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
2962 k*nDOF_trial_element+
2974 double penalty_floor,
2975 int nExteriorElementBoundaries_global,
2976 int nQuadraturePoints_elementBoundary,
2977 int nDOF_trial_element,
2980 int* exteriorElementBoundaries,
2981 int* elementBoundaryElements,
2982 int* elementBoundaryLocalElementBoundaries,
2992 double* fluxJacobian)
2994 int ebNE,ebN,eN_global,k,j,j_global,I,J,nSpace2=nSpace*nSpace;
2995 double Jacobian,diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2,max_a;
2996 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
2998 ebN = exteriorElementBoundaries[ebNE];
2999 eN_global = elementBoundaryElements[ebN*2+0];
3000 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3002 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] >= 1)
3004 for(j=0;j<nDOF_trial_element;j++)
3007 j_global = l2g[eN_global*nDOF_trial_element+j];
3009 for(I=0;I<nSpace;I++)
3011 diffusiveVelocityComponent_I_Jacobian=0.0;
3012 diffusiveVelocityComponent_I_Jacobian2=0.0;
3013 for(J=0;J<nSpace;J++)
3015 diffusiveVelocityComponent_I_Jacobian
3017 da[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
3022 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3025 diffusiveVelocityComponent_I_Jacobian2
3027 a[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
3032 grad_v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3033 k*nDOF_trial_element*nSpace+
3036 max_a = fmax(max_a,a[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
3044 (diffusiveVelocityComponent_I_Jacobian
3046 v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3047 k*nDOF_trial_element+
3050 diffusiveVelocityComponent_I_Jacobian2*
3053 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3057 max_a = fmax(penalty_floor,max_a);
3058 double penaltyJacobian = penalty[ebNE*nQuadraturePoints_elementBoundary+
3061 v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3062 k*nDOF_trial_element+
3064 if (scale_penalty) penaltyJacobian *= max_a;
3066 Jacobian += penaltyJacobian;
3068 fluxJacobian[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3069 k*nDOF_trial_element+
3078 double penalty_floor,
3079 int nExteriorElementBoundaries_global,
3080 int nQuadraturePoints_elementBoundary,
3081 int nDOF_trial_element,
3086 int* exteriorElementBoundaries,
3087 int* elementBoundaryElements,
3088 int* elementBoundaryLocalElementBoundaries,
3098 double* fluxJacobian)
3100 int ebNE,ebN,eN_global,k,j,j_global,I,m,
nnz=rowptr[nSpace];
3101 double Jacobian,diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2,max_a;
3102 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
3104 ebN = exteriorElementBoundaries[ebNE];
3105 eN_global = elementBoundaryElements[ebN*2+0];
3106 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3108 if(isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] >= 1)
3110 for(j=0;j<nDOF_trial_element;j++)
3113 j_global = l2g[eN_global*nDOF_trial_element+j];
3115 for(I=0;I<nSpace;I++)
3117 diffusiveVelocityComponent_I_Jacobian=0.0;
3118 diffusiveVelocityComponent_I_Jacobian2=0.0;
3119 for(m=rowptr[I];m<rowptr[I+1];m++)
3121 diffusiveVelocityComponent_I_Jacobian
3123 da[ebNE*nQuadraturePoints_elementBoundary*
nnz+
3127 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3130 diffusiveVelocityComponent_I_Jacobian2
3132 a[ebNE*nQuadraturePoints_elementBoundary*
nnz+
3136 grad_v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3137 k*nDOF_trial_element*nSpace+
3140 max_a = fmax(max_a,a[ebNE*nQuadraturePoints_elementBoundary*
nnz+
3147 (diffusiveVelocityComponent_I_Jacobian
3149 v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3150 k*nDOF_trial_element+
3153 diffusiveVelocityComponent_I_Jacobian2*
3156 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3160 max_a = fmax(penalty_floor,max_a);
3162 double penaltyJacobian = penalty[ebNE*nQuadraturePoints_elementBoundary+
3165 v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3166 k*nDOF_trial_element+
3168 if (scale_penalty) penaltyJacobian *= max_a;
3170 Jacobian += penaltyJacobian;
3172 fluxJacobian[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3173 k*nDOF_trial_element+
3594 int nElementBoundaries_element,
3595 int nQuadraturePoints_elementBoundary,
3596 int nDOF_trial_element,
3599 int* exteriorElementBoundaries,
3600 int* elementBoundaryElements,
3601 int* elementBoundaryLocalElementBoundaries,
3611 double* fluxJacobian)
3613 int ebNE,ebN,eN_global,ebN_element,k,j,j_global,I,J,nSpace2=nSpace*nSpace;
3614 double diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2;
3615 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
3617 ebN = exteriorElementBoundaries[ebNE];
3618 eN_global = elementBoundaryElements[ebN*2+0];
3619 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3620 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3622 for(j=0;j<nDOF_trial_element;j++)
3624 j_global = l2g[eN_global*nDOF_trial_element+j];
3625 for(I=0;I<nSpace;I++)
3627 diffusiveVelocityComponent_I_Jacobian=0.0;
3628 diffusiveVelocityComponent_I_Jacobian2=0.0;
3629 for(J=0;J<nSpace;J++)
3631 diffusiveVelocityComponent_I_Jacobian
3633 da[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
3634 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
3639 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3640 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3643 diffusiveVelocityComponent_I_Jacobian2
3645 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace2+
3646 ebN_element*nQuadraturePoints_elementBoundary*nSpace2+
3651 grad_v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3652 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3653 k*nDOF_trial_element*nSpace+
3657 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3658 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3659 k*nDOF_trial_element+
3661 (diffusiveVelocityComponent_I_Jacobian
3663 v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3664 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3665 k*nDOF_trial_element+
3668 diffusiveVelocityComponent_I_Jacobian2*
3671 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3672 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3681 int nElementBoundaries_element,
3682 int nQuadraturePoints_elementBoundary,
3683 int nDOF_trial_element,
3688 int* exteriorElementBoundaries,
3689 int* elementBoundaryElements,
3690 int* elementBoundaryLocalElementBoundaries,
3700 double* fluxJacobian)
3702 int ebNE,ebN,eN_global,ebN_element,k,j,j_global,I,m,
nnz=rowptr[nSpace];
3703 double diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2;
3704 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
3706 ebN = exteriorElementBoundaries[ebNE];
3707 eN_global = elementBoundaryElements[ebN*2+0];
3708 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3709 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3711 for(j=0;j<nDOF_trial_element;j++)
3713 j_global = l2g[eN_global*nDOF_trial_element+j];
3714 for(I=0;I<nSpace;I++)
3716 diffusiveVelocityComponent_I_Jacobian=0.0;
3717 diffusiveVelocityComponent_I_Jacobian2=0.0;
3718 for(m=rowptr[I];m<rowptr[I+1];m++)
3720 diffusiveVelocityComponent_I_Jacobian
3722 da[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
3723 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
3727 grad_phi[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3728 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3731 diffusiveVelocityComponent_I_Jacobian2
3733 a[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*
nnz+
3734 ebN_element*nQuadraturePoints_elementBoundary*
nnz+
3738 grad_v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3739 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3740 k*nDOF_trial_element*nSpace+
3744 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3745 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3746 k*nDOF_trial_element+
3748 (diffusiveVelocityComponent_I_Jacobian
3750 v[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3751 ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3752 k*nDOF_trial_element+
3755 diffusiveVelocityComponent_I_Jacobian2*
3758 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3759 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3769 int nQuadraturePoints_elementBoundary,
3770 int nDOF_trial_element,
3773 int* exteriorElementBoundaries,
3774 int* elementBoundaryElements,
3775 int* elementBoundaryLocalElementBoundaries,
3785 double* fluxJacobian)
3787 int ebNE,ebN,eN_global,k,j,j_global,I,J,nSpace2=nSpace*nSpace;
3788 double diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2;
3789 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
3791 ebN = exteriorElementBoundaries[ebNE];
3792 eN_global = elementBoundaryElements[ebN*2+0];
3793 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3795 for(j=0;j<nDOF_trial_element;j++)
3797 j_global = l2g[eN_global*nDOF_trial_element+j];
3798 for(I=0;I<nSpace;I++)
3800 diffusiveVelocityComponent_I_Jacobian=0.0;
3801 diffusiveVelocityComponent_I_Jacobian2=0.0;
3802 for(J=0;J<nSpace;J++)
3804 diffusiveVelocityComponent_I_Jacobian
3806 da[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
3811 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3814 diffusiveVelocityComponent_I_Jacobian2
3816 a[ebNE*nQuadraturePoints_elementBoundary*nSpace2+
3821 grad_v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3822 k*nDOF_trial_element*nSpace+
3826 fluxJacobian[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3827 k*nDOF_trial_element+
3829 (diffusiveVelocityComponent_I_Jacobian
3831 v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3832 k*nDOF_trial_element+
3835 diffusiveVelocityComponent_I_Jacobian2*
3838 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3847 int nQuadraturePoints_elementBoundary,
3848 int nDOF_trial_element,
3853 int* exteriorElementBoundaries,
3854 int* elementBoundaryElements,
3855 int* elementBoundaryLocalElementBoundaries,
3865 double* fluxJacobian)
3867 int ebNE,ebN,eN_global,k,j,j_global,I,m,
nnz=rowptr[nSpace];
3868 double diffusiveVelocityComponent_I_Jacobian,diffusiveVelocityComponent_I_Jacobian2;
3869 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
3871 ebN = exteriorElementBoundaries[ebNE];
3872 eN_global = elementBoundaryElements[ebN*2+0];
3873 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3875 for(j=0;j<nDOF_trial_element;j++)
3877 j_global = l2g[eN_global*nDOF_trial_element+j];
3878 for(I=0;I<nSpace;I++)
3880 diffusiveVelocityComponent_I_Jacobian=0.0;
3881 diffusiveVelocityComponent_I_Jacobian2=0.0;
3882 for(m=rowptr[I];m<rowptr[I+1];m++)
3884 diffusiveVelocityComponent_I_Jacobian
3886 da[ebNE*nQuadraturePoints_elementBoundary*
nnz+
3890 grad_phi[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3893 diffusiveVelocityComponent_I_Jacobian2
3895 a[ebNE*nQuadraturePoints_elementBoundary*
nnz+
3899 grad_v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element*nSpace+
3900 k*nDOF_trial_element*nSpace+
3904 fluxJacobian[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3905 k*nDOF_trial_element+
3907 (diffusiveVelocityComponent_I_Jacobian
3909 v[ebNE*nQuadraturePoints_elementBoundary*nDOF_trial_element+
3910 k*nDOF_trial_element+
3913 diffusiveVelocityComponent_I_Jacobian2*
3916 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
3928 int nElementBoundaries_element,
3929 int nQuadraturePoints_elementBoundary,
3931 int* interiorElementBoundaries,
3932 int* elementBoundaryElements,
3933 int* elementBoundaryLocalElementBoundaries,
3940 double* dflux_right)
3942 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,J;
3943 double left_speed,right_speed,left_flux,right_flux,shock_speed,flux_jump,u_jump;
3944 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
3946 ebN = interiorElementBoundaries[ebNI];
3947 left_eN_global = elementBoundaryElements[ebN*2+0];
3948 right_eN_global = elementBoundaryElements[ebN*2+1];
3949 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
3950 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
3951 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
3959 for(J=0;J<nSpace;J++)
3963 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3964 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3968 df[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3969 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3974 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3975 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3979 df[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3980 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3985 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3986 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3990 f[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3991 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
3996 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
3997 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4001 f[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4002 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4006 flux_jump = (right_flux - left_flux);
4007 u_jump = (
u[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4008 right_ebN_element*nQuadraturePoints_elementBoundary+
4011 u[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4012 left_ebN_element*nQuadraturePoints_elementBoundary+
4014 if (fabs(u_jump) > fabs(flux_jump*1.0e-16))
4015 shock_speed = flux_jump/u_jump;
4016 if (left_speed >= 0.0 && right_speed >= 0.0)
4018 flux[ebN*nQuadraturePoints_elementBoundary+
4020 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4022 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4025 else if (left_speed <= 0.0 && right_speed <= 0.0)
4027 flux[ebN*nQuadraturePoints_elementBoundary+
4029 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4031 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4034 else if (left_speed >= 0.0 && right_speed <= 0.0)
4036 if (shock_speed >= 0.0)
4038 flux[ebN*nQuadraturePoints_elementBoundary+
4040 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4042 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4047 flux[ebN*nQuadraturePoints_elementBoundary+
4049 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4051 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4058 printf(
"Transonic rarefaction detected in interior. This numerical flux treats transonic rarefactions incorrectly. left_speed= %12.5e right_speed= %12.5e \n",left_speed,right_speed);
4059 for (J=0; J < nSpace; J++)
4061 printf(
"n_l[%d]=%g df_l[%d]=%g df_r[%d]=%g \n",J,
4062 n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4063 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4067 df[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4068 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4072 df[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4073 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4078 flux[ebN*nQuadraturePoints_elementBoundary+
4079 k] = 0.5*(left_flux + right_flux);
4080 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4081 k] = 0.5*left_speed;
4082 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4083 k] = 0.5*right_speed;
4093 int nElementBoundaries_element,
4094 int nQuadraturePoints_elementBoundary,
4095 int nDOF_trial_element,
4096 int* interiorElementBoundaries,
4097 int* elementBoundaryElements,
4098 int* elementBoundaryLocalElementBoundaries,
4100 double* dflux_right,
4102 double* fluxJacobian)
4104 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,j;
4105 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
4107 ebN = interiorElementBoundaries[ebNI];
4108 left_eN_global = elementBoundaryElements[ebN*2+0];
4109 right_eN_global = elementBoundaryElements[ebN*2+1];
4110 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4111 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
4112 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4113 for(j=0;j<nDOF_trial_element;j++)
4115 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4116 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4117 k*nDOF_trial_element+
4120 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4123 v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4124 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4125 k*nDOF_trial_element+
4127 fluxJacobian[ebN*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4128 1*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4129 k*nDOF_trial_element+
4132 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4135 v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4136 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4137 k*nDOF_trial_element+
4146 int nElementBoundaries_element,
4147 int nQuadraturePoints_elementBoundary,
4148 int nDOF_trial_element,
4149 int* interiorElementBoundaries,
4150 int* elementBoundaryElements,
4151 int* elementBoundaryLocalElementBoundaries,
4153 double* dflux_right,
4155 double* fluxJacobian_2sided)
4157 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,j;
4158 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
4160 ebN = interiorElementBoundaries[ebNI];
4161 left_eN_global = elementBoundaryElements[ebN*2+0];
4162 right_eN_global = elementBoundaryElements[ebN*2+1];
4163 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4164 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
4165 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4166 for(j=0;j<nDOF_trial_element;j++)
4169 fluxJacobian_2sided[ebN*2*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4170 0*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4171 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4172 k*nDOF_trial_element+
4175 dflux_left[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4176 left_ebN_element*nQuadraturePoints_elementBoundary+
4179 v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4180 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4181 k*nDOF_trial_element+
4183 fluxJacobian_2sided[ebN*2*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4184 0*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4185 1*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4186 k*nDOF_trial_element+
4189 dflux_right[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4190 left_ebN_element*nQuadraturePoints_elementBoundary+
4193 v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4194 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4195 k*nDOF_trial_element+
4198 fluxJacobian_2sided[ebN*2*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4199 1*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4200 0*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4201 k*nDOF_trial_element+
4204 dflux_left[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4205 right_ebN_element*nQuadraturePoints_elementBoundary+
4208 v[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4209 left_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4210 k*nDOF_trial_element+
4212 fluxJacobian_2sided[ebN*2*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4213 1*2*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4214 1*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4215 k*nDOF_trial_element+
4218 dflux_right[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4219 right_ebN_element*nQuadraturePoints_elementBoundary+
4222 v[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4223 right_ebN_element*nQuadraturePoints_elementBoundary*nDOF_trial_element+
4224 k*nDOF_trial_element+
4234 int nElementBoundaries_element,
4235 int nQuadraturePoints_elementBoundary,
4237 int* interiorElementBoundaries,
4238 int* elementBoundaryElements,
4239 int* elementBoundaryLocalElementBoundaries,
4246 double* dflux_right)
4248 int ebNI,ebN,left_eN_global,right_eN_global,left_ebN_element,right_ebN_element,k,J;
4249 for(ebNI=0;ebNI<nInteriorElementBoundaries_global;ebNI++)
4251 ebN = interiorElementBoundaries[ebNI];
4252 left_eN_global = elementBoundaryElements[ebN*2+0];
4253 right_eN_global = elementBoundaryElements[ebN*2+1];
4254 left_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4255 right_ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+1];
4256 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4258 flux[ebN*nQuadraturePoints_elementBoundary+
4260 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4262 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4264 for(J=0;J<nSpace;J++)
4266 flux[ebN*nQuadraturePoints_elementBoundary+
4267 k] +=
n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4268 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4272 (
f[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4273 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4277 f[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4278 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4281 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4282 k] +=
n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4283 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4287 df[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4288 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4291 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4292 k] +=
n[left_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4293 left_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4297 df[right_eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4298 right_ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4302 flux[ebN*nQuadraturePoints_elementBoundary+
4304 dflux_left[ebN*nQuadraturePoints_elementBoundary+
4306 dflux_right[ebN*nQuadraturePoints_elementBoundary+
4316 int nElementBoundaries_element,
4317 int nQuadraturePoints_elementBoundary,
4319 int* exteriorElementBoundaries,
4320 int* elementBoundaryElements,
4321 int* elementBoundaryLocalElementBoundaries,
4329 int ebNE,ebN,eN_global,ebN_element,k,J;
4330 memset(inflowFlag,0,
sizeof(
int)*nExteriorElementBoundaries_global*nQuadraturePoints_elementBoundary);
4331 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4333 ebN = exteriorElementBoundaries[ebNE];
4334 eN_global = elementBoundaryElements[ebN*2+0];
4335 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4336 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4338 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
4339 dflux_left[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
4340 inflowFlag[ebNE*nQuadraturePoints_elementBoundary+k]=0;
4341 for(J=0;J<nSpace;J++)
4343 flux[ebN*nQuadraturePoints_elementBoundary+k]
4345 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4346 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4350 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4351 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4354 dflux_left[ebN*nQuadraturePoints_elementBoundary+k]
4356 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4357 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4361 df[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4362 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4366 if(dflux_left[ebN*nQuadraturePoints_elementBoundary+k] < 0.0)
4369 inflowFlag[ebNE*nQuadraturePoints_elementBoundary+k] = 1;
4370 flux[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
4371 dflux_left[ebN*nQuadraturePoints_elementBoundary+k] = 0.0;
4380 int nQuadraturePoints_elementBoundary,
4382 int* exteriorElementBoundaries,
4383 int* elementBoundaryElements,
4384 int* elementBoundaryLocalElementBoundaries,
4392 int ebNE,ebN,eN_global,ebN_element,k,J;
4393 memset(inflowFlag,0,
sizeof(
int)*nExteriorElementBoundaries_global*nQuadraturePoints_elementBoundary);
4394 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4396 ebN = exteriorElementBoundaries[ebNE];
4397 eN_global = elementBoundaryElements[ebN*2+0];
4398 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4399 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4401 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
4402 dflux_left[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
4403 inflowFlag[ebNE*nQuadraturePoints_elementBoundary+k]=0;
4404 for(J=0;J<nSpace;J++)
4406 flux[ebNE*nQuadraturePoints_elementBoundary+k]
4408 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4412 f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4415 dflux_left[ebNE*nQuadraturePoints_elementBoundary+k]
4417 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4421 df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4425 if(dflux_left[ebNE*nQuadraturePoints_elementBoundary+k] < 0.0)
4428 inflowFlag[ebNE*nQuadraturePoints_elementBoundary+k] = 1;
4429 flux[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
4430 dflux_left[ebNE*nQuadraturePoints_elementBoundary+k] = 0.0;
4440 int nElementBoundaries_element,
4441 int nQuadraturePoints_elementBoundary,
4443 int* exteriorElementBoundaries,
4444 int* elementBoundaryElements,
4445 int* elementBoundaryLocalElementBoundaries,
4458 int ebNE,ebN,eN_global,ebN_element,k,J;
4459 double left_speed,right_speed,left_flux,right_flux,shock_speed,flux_jump,u_jump;
4460 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4462 ebN = exteriorElementBoundaries[ebNE];
4463 eN_global = elementBoundaryElements[ebN*2+0];
4464 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4465 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4473 for(J=0;J<nSpace;J++)
4477 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4478 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4482 df[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4483 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4488 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4489 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4493 bc_df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4498 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4499 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4503 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4504 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4509 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4510 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4514 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4518 flux_jump = (right_flux - left_flux);
4519 u_jump = (bc_u[ebNE*nQuadraturePoints_elementBoundary+
4522 u[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
4523 ebN_element*nQuadraturePoints_elementBoundary+
4525 if (fabs(u_jump) > fabs(flux_jump*1.0e-16))
4526 shock_speed = flux_jump/u_jump;
4527 if (left_speed >= 0.0 && right_speed >= 0.0)
4529 flux[ebN*nQuadraturePoints_elementBoundary+
4531 dflux[ebN*nQuadraturePoints_elementBoundary+
4534 else if (left_speed <= 0.0 && right_speed <= 0.0)
4536 flux[ebN*nQuadraturePoints_elementBoundary+
4538 if (isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
4539 dflux[ebN*nQuadraturePoints_elementBoundary+
4542 dflux[ebN*nQuadraturePoints_elementBoundary+
4545 else if (left_speed >= 0.0 && right_speed <= 0.0)
4547 if (shock_speed >= 0.0)
4549 flux[ebN*nQuadraturePoints_elementBoundary+
4551 dflux[ebN*nQuadraturePoints_elementBoundary+
4556 flux[ebN*nQuadraturePoints_elementBoundary+
4558 if (isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
4559 dflux[ebN*nQuadraturePoints_elementBoundary+
4562 dflux[ebN*nQuadraturePoints_elementBoundary+
4569 printf(
"Transonic rarefaction detected on exterior. This numerical flux treats transonic rarefactions incorrectly. left_speed= %12.5e right_speed= %12.5e \n",left_speed,right_speed);
4570 for (J=0; J < nSpace; J++)
4572 printf(
"n_l[%d]=%g df_l[%d]=%g df_r[%d]=%g \n",J,
4573 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4574 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4578 df[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4579 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4583 bc_df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4588 flux[ebN*nQuadraturePoints_elementBoundary+
4589 k] = 0.5*(left_flux + right_flux);
4590 if (isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
4591 dflux[ebN*nQuadraturePoints_elementBoundary+
4592 k] = 0.5*left_speed;
4594 dflux[ebN*nQuadraturePoints_elementBoundary+
4595 k] = 0.5*(left_speed+right_speed);
4633 int nQuadraturePoints_elementBoundary,
4635 int* exteriorElementBoundaries,
4636 int* elementBoundaryElements,
4637 int* elementBoundaryLocalElementBoundaries,
4650 int ebNE,ebN,eN_global,k,J;
4651 double left_speed,right_speed,left_flux,right_flux,shock_speed,flux_jump,u_jump;
4652 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4654 ebN = exteriorElementBoundaries[ebNE];
4655 eN_global = elementBoundaryElements[ebN*2+0];
4656 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4663 for(J=0;J<nSpace;J++)
4667 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4671 df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4676 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4680 bc_df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4685 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4689 f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4694 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4698 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4702 flux_jump = (right_flux - left_flux);
4703 u_jump = (bc_u[ebNE*nQuadraturePoints_elementBoundary+
4706 u[ebNE*nQuadraturePoints_elementBoundary+
4708 if (fabs(u_jump) > fabs(flux_jump*1.0e-16))
4709 shock_speed = flux_jump/u_jump;
4710 if (left_speed >= 0.0 && right_speed >= 0.0)
4712 flux[ebNE*nQuadraturePoints_elementBoundary+
4714 dflux[ebNE*nQuadraturePoints_elementBoundary+
4717 else if (left_speed <= 0.0 && right_speed <= 0.0)
4719 flux[ebNE*nQuadraturePoints_elementBoundary+
4721 if (isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
4722 dflux[ebNE*nQuadraturePoints_elementBoundary+
4725 dflux[ebNE*nQuadraturePoints_elementBoundary+
4728 else if (left_speed >= 0.0 && right_speed <= 0.0)
4730 if (shock_speed >= 0.0)
4732 flux[ebNE*nQuadraturePoints_elementBoundary+
4734 dflux[ebNE*nQuadraturePoints_elementBoundary+
4739 flux[ebNE*nQuadraturePoints_elementBoundary+
4741 if (isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
4742 dflux[ebNE*nQuadraturePoints_elementBoundary+
4745 dflux[ebNE*nQuadraturePoints_elementBoundary+
4752 printf(
"Transonic rarefaction detected on exterior. This numerical flux treats transonic rarefactions incorrectly. left_speed= %12.5e right_speed= %12.5e \n",left_speed,right_speed);
4753 for (J=0; J < nSpace; J++)
4755 printf(
"n_l[%d]=%g df_l[%d]=%g df_r[%d]=%g \n",J,
4756 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4760 df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4764 bc_df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4769 flux[ebNE*nQuadraturePoints_elementBoundary+
4770 k] = 0.5*(left_flux + right_flux);
4771 if (isDOFBoundary[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
4772 dflux[ebNE*nQuadraturePoints_elementBoundary+
4773 k] = 0.5*left_speed;
4775 dflux[ebNE*nQuadraturePoints_elementBoundary+
4776 k] = 0.5*(left_speed+right_speed);
4790 int nElementBoundaries_element,
4791 int nQuadraturePoints_elementBoundary,
4793 int* exteriorElementBoundaries,
4794 int* elementBoundaryElements,
4795 int* elementBoundaryLocalElementBoundaries,
4808 int ebNE,ebN,eN_global,ebN_element,k,J;
4809 double left_speed,right_speed,left_flux,right_flux,shock_speed;
4810 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4812 ebN = exteriorElementBoundaries[ebNE];
4813 eN_global = elementBoundaryElements[ebN*2+0];
4814 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4815 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4823 for(J=0;J<nSpace;J++)
4827 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4828 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4832 df[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4833 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4838 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4839 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4843 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4844 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4848 flux[ebN*nQuadraturePoints_elementBoundary+
4850 dflux[ebN*nQuadraturePoints_elementBoundary+
4857 int nQuadraturePoints_elementBoundary,
4859 int* exteriorElementBoundaries,
4860 int* elementBoundaryElements,
4861 int* elementBoundaryLocalElementBoundaries,
4874 int ebNE,ebN,eN_global,k,J;
4875 double left_speed,right_speed,left_flux,right_flux,shock_speed;
4876 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4878 ebN = exteriorElementBoundaries[ebNE];
4879 eN_global = elementBoundaryElements[ebN*2+0];
4880 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4887 for(J=0;J<nSpace;J++)
4891 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4895 df[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4900 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4904 f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
4908 flux[ebNE*nQuadraturePoints_elementBoundary+
4910 dflux[ebNE*nQuadraturePoints_elementBoundary+
4919 int nElementBoundaries_element,
4920 int nQuadraturePoints_elementBoundary,
4922 int* exteriorElementBoundaries,
4923 int* elementBoundaryElements,
4924 int* elementBoundaryLocalElementBoundaries,
4925 int *isDOFBoundary_p,
4926 int *isDOFBoundary_u,
4927 int *isDOFBoundary_v,
4947 int ebNE,ebN,eN_global,ebN_element,k;
4948 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
4950 ebN = exteriorElementBoundaries[ebNE];
4951 eN_global = elementBoundaryElements[ebN*2+0];
4952 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
4953 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
4955 fluxpu[ebN*nQuadraturePoints_elementBoundary+
4957 fluxpv[ebN*nQuadraturePoints_elementBoundary+
4959 flux[ebN*nQuadraturePoints_elementBoundary+
4963 if (isDOFBoundary_p[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
4965 dfluxpu_dp[ebN*nQuadraturePoints_elementBoundary+
4968 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4969 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4973 dfpu_dp[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4974 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4977 fluxpu[ebN*nQuadraturePoints_elementBoundary+
4979 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4980 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4984 fpu[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4985 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4988 dfluxpv_dp[ebN*nQuadraturePoints_elementBoundary+
4991 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4992 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
4996 dfpv_dp[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
4997 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5000 fluxpv[ebN*nQuadraturePoints_elementBoundary+
5002 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5003 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5007 fpv[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5008 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5014 fluxpu[ebN*nQuadraturePoints_elementBoundary+
5016 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5017 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5021 bc_fpu[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5024 fluxpv[ebN*nQuadraturePoints_elementBoundary+
5026 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5027 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5031 bc_fpv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5036 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5038 dflux_du[ebN*nQuadraturePoints_elementBoundary+
5041 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5042 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5046 df_du[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5047 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5050 flux[ebN*nQuadraturePoints_elementBoundary+
5052 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5053 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5057 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5058 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5064 flux[ebN*nQuadraturePoints_elementBoundary+
5066 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5067 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5071 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5075 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5077 dflux_dv[ebN*nQuadraturePoints_elementBoundary+
5080 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5081 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5085 df_dv[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5086 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5089 flux[ebN*nQuadraturePoints_elementBoundary+
5091 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5092 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5096 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5097 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5103 flux[ebN*nQuadraturePoints_elementBoundary+
5105 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5106 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5110 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5118 int nElementBoundaries_element,
5119 int nQuadraturePoints_elementBoundary,
5121 int* exteriorElementBoundaries,
5122 int* elementBoundaryElements,
5123 int* elementBoundaryLocalElementBoundaries,
5124 int *isDOFBoundary_p,
5125 int *isDOFBoundary_u,
5126 int *isDOFBoundary_v,
5145 double* dflux_mass_du,
5146 double* dflux_mass_dv,
5147 double* dflux_umom_dp,
5148 double* dflux_umom_du,
5149 double* dflux_umom_dv,
5150 double* dflux_vmom_dp,
5151 double* dflux_vmom_du,
5152 double* dflux_vmom_dv)
5154 int ebNE,ebN,eN_global,ebN_element,k;
5155 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5157 ebN = exteriorElementBoundaries[ebNE];
5158 eN_global = elementBoundaryElements[ebN*2+0];
5159 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
5160 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5162 flux_umom[ebN*nQuadraturePoints_elementBoundary+
5164 flux_vmom[ebN*nQuadraturePoints_elementBoundary+
5166 flux_mass[ebN*nQuadraturePoints_elementBoundary+
5168 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5170 dflux_mass_du[ebN*nQuadraturePoints_elementBoundary+
5173 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5174 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5178 df_mass_du[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5179 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5182 flux_mass[ebN*nQuadraturePoints_elementBoundary+
5184 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5185 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5189 f_mass[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5190 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5193 dflux_umom_du[ebN*nQuadraturePoints_elementBoundary+
5196 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5197 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5201 df_umom_du[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5202 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5205 flux_umom[ebN*nQuadraturePoints_elementBoundary+
5208 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5209 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5213 f_umom[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5214 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5217 dflux_vmom_du[ebN*nQuadraturePoints_elementBoundary+
5220 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5221 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5225 df_vmom_du[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5226 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5229 flux_vmom[ebN*nQuadraturePoints_elementBoundary+
5232 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5233 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5237 f_vmom[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5238 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5244 flux_mass[ebN*nQuadraturePoints_elementBoundary+
5247 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5248 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5252 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5255 flux_umom[ebN*nQuadraturePoints_elementBoundary+
5258 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5259 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5263 bc_f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5266 flux_vmom[ebN*nQuadraturePoints_elementBoundary+
5269 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5270 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5274 bc_f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5278 if (isDOFBoundary_p[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
5280 dflux_umom_dp[ebN*nQuadraturePoints_elementBoundary+
5283 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5284 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5287 flux_umom[ebN*nQuadraturePoints_elementBoundary+
5290 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5291 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5295 (p[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5296 ebN_element*nQuadraturePoints_elementBoundary+
5299 bc_p[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5300 ebN_element*nQuadraturePoints_elementBoundary+
5302 dflux_vmom_dp[ebN*nQuadraturePoints_elementBoundary+
5305 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5306 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5309 flux_vmom[ebN*nQuadraturePoints_elementBoundary+
5311 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
5312 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
5316 (p[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5317 ebN_element*nQuadraturePoints_elementBoundary+
5320 bc_p[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary+
5321 ebN_element*nQuadraturePoints_elementBoundary+
5329 int nQuadraturePoints_elementBoundary,
5331 int* exteriorElementBoundaries,
5332 int* elementBoundaryElements,
5333 int* elementBoundaryLocalElementBoundaries,
5334 int *isDOFBoundary_p,
5335 int *isDOFBoundary_u,
5336 int *isDOFBoundary_v,
5358 double* dflux_mass_dp,
5359 double* dflux_mass_du,
5360 double* dflux_mass_dv,
5361 double* dflux_umom_dp,
5362 double* dflux_umom_du,
5363 double* dflux_umom_dv,
5364 double* dflux_vmom_dp,
5365 double* dflux_vmom_du,
5366 double* dflux_vmom_dv,
5370 double flowDirection;
5371 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5373 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5375 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5377 dflux_mass_du[ebNE*nQuadraturePoints_elementBoundary+
5379 dflux_mass_dv[ebNE*nQuadraturePoints_elementBoundary+
5382 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5384 dflux_umom_dp[ebNE*nQuadraturePoints_elementBoundary+
5386 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
5388 dflux_umom_dv[ebNE*nQuadraturePoints_elementBoundary+
5391 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5393 dflux_vmom_dp[ebNE*nQuadraturePoints_elementBoundary+
5395 dflux_vmom_du[ebNE*nQuadraturePoints_elementBoundary+
5397 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
5400 flowDirection=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5404 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5408 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5412 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5415 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5417 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5420 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5424 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5427 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5431 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5434 dflux_mass_du[ebNE*nQuadraturePoints_elementBoundary+
5437 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5441 df_mass_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5444 if (flowDirection >= 0.0)
5446 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5449 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5453 f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5456 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
5459 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5463 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5466 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5469 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5473 f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5476 dflux_vmom_du[ebNE*nQuadraturePoints_elementBoundary+
5479 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5483 df_vmom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5486 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
5489 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5493 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5500 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5503 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5507 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5510 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5514 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5517 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5520 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5524 bc_f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5527 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5530 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5534 bc_f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5537 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5538 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
5541 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5545 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5549 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5551 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5554 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5558 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5561 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5565 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5568 dflux_mass_dv[ebNE*nQuadraturePoints_elementBoundary+
5571 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5575 df_mass_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5578 if (flowDirection >= 0.0)
5580 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5583 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5587 f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5590 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
5593 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5597 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5600 dflux_umom_dv[ebNE*nQuadraturePoints_elementBoundary+
5603 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5607 df_umom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5610 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5613 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5617 f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5620 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
5622 +=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5626 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5633 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5636 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5640 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5643 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5647 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5651 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5654 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5658 bc_f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5661 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5662 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
5665 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5669 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5672 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5675 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5679 bc_f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5683 if (isDOFBoundary_p[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
5698 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5701 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5705 (bc_p[ebNE*nQuadraturePoints_elementBoundary+
5708 p[ebNE*nQuadraturePoints_elementBoundary+
5710 dflux_umom_dp[ebNE*nQuadraturePoints_elementBoundary+
5712 = -
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5715 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5717 +=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5721 (bc_p[ebNE*nQuadraturePoints_elementBoundary+
5724 p[ebNE*nQuadraturePoints_elementBoundary+
5726 dflux_vmom_dp[ebNE*nQuadraturePoints_elementBoundary+
5728 = -
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5736 int nQuadraturePoints_elementBoundary,
5738 int* exteriorElementBoundaries,
5739 int* elementBoundaryElements,
5740 int* elementBoundaryLocalElementBoundaries,
5741 int *isDOFBoundary_p,
5742 int *isDOFBoundary_u,
5743 int *isDOFBoundary_v,
5744 int *isDOFBoundary_w,
5775 double* dflux_mass_du,
5776 double* dflux_mass_dv,
5777 double* dflux_mass_dw,
5778 double* dflux_umom_dp,
5779 double* dflux_umom_du,
5780 double* dflux_umom_dv,
5781 double* dflux_umom_dw,
5782 double* dflux_vmom_dp,
5783 double* dflux_vmom_du,
5784 double* dflux_vmom_dv,
5785 double* dflux_vmom_dw,
5786 double* dflux_wmom_dp,
5787 double* dflux_wmom_du,
5788 double* dflux_wmom_dv,
5789 double* dflux_wmom_dw,
5793 double flowDirection;
5794 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
5796 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
5798 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5800 dflux_mass_du[ebNE*nQuadraturePoints_elementBoundary+
5802 dflux_mass_dv[ebNE*nQuadraturePoints_elementBoundary+
5804 dflux_mass_dw[ebNE*nQuadraturePoints_elementBoundary+
5807 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5809 dflux_umom_dp[ebNE*nQuadraturePoints_elementBoundary+
5811 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
5813 dflux_umom_dv[ebNE*nQuadraturePoints_elementBoundary+
5815 dflux_umom_dw[ebNE*nQuadraturePoints_elementBoundary+
5818 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5820 dflux_vmom_dp[ebNE*nQuadraturePoints_elementBoundary+
5822 dflux_vmom_du[ebNE*nQuadraturePoints_elementBoundary+
5824 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
5826 dflux_vmom_dw[ebNE*nQuadraturePoints_elementBoundary+
5829 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
5831 dflux_wmom_dp[ebNE*nQuadraturePoints_elementBoundary+
5833 dflux_wmom_du[ebNE*nQuadraturePoints_elementBoundary+
5835 dflux_wmom_dv[ebNE*nQuadraturePoints_elementBoundary+
5837 dflux_wmom_dw[ebNE*nQuadraturePoints_elementBoundary+
5840 flowDirection=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5844 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5848 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5852 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5855 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5859 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5862 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
5864 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5867 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5871 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5874 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5878 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5881 dflux_mass_du[ebNE*nQuadraturePoints_elementBoundary+
5884 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5888 df_mass_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5891 if (flowDirection >= 0.0)
5893 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5896 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5900 f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5903 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
5906 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5910 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5913 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
5916 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5920 f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5923 dflux_vmom_du[ebNE*nQuadraturePoints_elementBoundary+
5926 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5930 df_vmom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5933 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
5936 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5940 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5943 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
5946 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5950 f_wmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5953 dflux_wmom_du[ebNE*nQuadraturePoints_elementBoundary+
5956 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5960 df_wmom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5963 dflux_wmom_dw[ebNE*nQuadraturePoints_elementBoundary+
5966 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5970 df_wmom_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5977 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
5980 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5984 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5987 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5991 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
5994 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
5997 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6001 bc_f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6004 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
6007 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6011 bc_f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6014 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6015 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
6018 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6022 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6025 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
6028 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6032 bc_f_wmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6035 if (isDOFBoundary_w[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6036 dflux_wmom_dw[ebNE*nQuadraturePoints_elementBoundary+
6039 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6043 df_wmom_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6047 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6049 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
6052 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6056 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6059 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6063 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6066 dflux_mass_dv[ebNE*nQuadraturePoints_elementBoundary+
6069 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6073 df_mass_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6076 if (flowDirection >= 0.0)
6078 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
6081 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6085 f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6088 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
6091 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6095 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6098 dflux_umom_dv[ebNE*nQuadraturePoints_elementBoundary+
6101 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6105 df_umom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6108 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
6111 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6115 f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6118 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
6120 +=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6124 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6127 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
6130 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6134 f_wmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6137 dflux_wmom_dw[ebNE*nQuadraturePoints_elementBoundary+
6140 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6144 df_wmom_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6147 dflux_wmom_dv[ebNE*nQuadraturePoints_elementBoundary+
6150 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6154 df_wmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6161 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
6164 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6168 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6171 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6175 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6178 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
6181 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6185 bc_f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6188 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6189 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
6192 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6196 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6199 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
6202 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6206 bc_f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6209 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
6212 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6216 bc_f_wmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6219 if (isDOFBoundary_w[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6220 dflux_wmom_dw[ebNE*nQuadraturePoints_elementBoundary+
6223 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6227 df_wmom_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6231 if (isDOFBoundary_w[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6233 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
6236 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6240 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6243 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6247 f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6250 dflux_mass_dw[ebNE*nQuadraturePoints_elementBoundary+
6253 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6257 df_mass_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6260 if (flowDirection >= 0.0)
6262 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
6265 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6269 f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6272 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
6275 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6279 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6282 dflux_umom_dw[ebNE*nQuadraturePoints_elementBoundary+
6285 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6289 df_umom_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6292 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
6295 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6299 f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6302 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
6304 +=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6308 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6311 dflux_vmom_dw[ebNE*nQuadraturePoints_elementBoundary+
6313 +=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6317 df_vmom_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6320 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
6323 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6327 f_wmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6330 dflux_wmom_dw[ebNE*nQuadraturePoints_elementBoundary+
6333 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6337 df_wmom_dw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6344 flux_mass[ebNE*nQuadraturePoints_elementBoundary+
6347 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6351 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6354 velocity[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6358 bc_f_mass[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6361 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
6364 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6368 bc_f_umom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6371 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6372 dflux_umom_du[ebNE*nQuadraturePoints_elementBoundary+
6375 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6379 df_umom_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6382 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
6385 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6389 bc_f_vmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6392 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6393 dflux_vmom_dv[ebNE*nQuadraturePoints_elementBoundary+
6396 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6400 df_vmom_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6403 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
6406 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6410 bc_f_wmom[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6414 if (isDOFBoundary_p[ebNE*nQuadraturePoints_elementBoundary+k] == 1)
6416 flux_umom[ebNE*nQuadraturePoints_elementBoundary+
6419 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6423 (bc_p[ebNE*nQuadraturePoints_elementBoundary+
6426 p[ebNE*nQuadraturePoints_elementBoundary+
6428 dflux_umom_dp[ebNE*nQuadraturePoints_elementBoundary+
6430 = -
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6433 flux_vmom[ebNE*nQuadraturePoints_elementBoundary+
6435 +=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6439 (bc_p[ebNE*nQuadraturePoints_elementBoundary+
6442 p[ebNE*nQuadraturePoints_elementBoundary+
6444 dflux_vmom_dp[ebNE*nQuadraturePoints_elementBoundary+
6446 = -
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6449 flux_wmom[ebNE*nQuadraturePoints_elementBoundary+
6451 +=
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6455 (bc_p[ebNE*nQuadraturePoints_elementBoundary+
6458 p[ebNE*nQuadraturePoints_elementBoundary+
6460 dflux_wmom_dp[ebNE*nQuadraturePoints_elementBoundary+
6462 = -
n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6473 int nQuadraturePoints_elementBoundary,
6475 int* exteriorElementBoundaries,
6476 int* elementBoundaryElements,
6477 int* elementBoundaryLocalElementBoundaries,
6478 int *isDOFBoundary_p,
6479 int *isDOFBoundary_u,
6480 int *isDOFBoundary_v,
6501 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
6503 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
6505 fluxpu[ebNE*nQuadraturePoints_elementBoundary+
6507 fluxpv[ebNE*nQuadraturePoints_elementBoundary+
6509 flux[ebNE*nQuadraturePoints_elementBoundary+
6513 if (isDOFBoundary_p[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6515 dfluxpu_dp[ebNE*nQuadraturePoints_elementBoundary+
6518 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6522 dfpu_dp[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6525 fluxpu[ebNE*nQuadraturePoints_elementBoundary+
6527 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6531 fpu[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6534 dfluxpv_dp[ebNE*nQuadraturePoints_elementBoundary+
6537 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6541 dfpv_dp[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6544 fluxpv[ebNE*nQuadraturePoints_elementBoundary+
6546 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6550 fpv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6556 fluxpu[ebNE*nQuadraturePoints_elementBoundary+
6558 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6562 bc_fpu[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6565 fluxpv[ebNE*nQuadraturePoints_elementBoundary+
6567 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6571 bc_fpv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6576 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6578 dflux_du[ebNE*nQuadraturePoints_elementBoundary+
6581 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6585 df_du[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6588 flux[ebNE*nQuadraturePoints_elementBoundary+
6590 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6594 f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6600 flux[ebNE*nQuadraturePoints_elementBoundary+
6602 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6606 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6610 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6612 dflux_dv[ebNE*nQuadraturePoints_elementBoundary+
6615 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6619 df_dv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6622 flux[ebNE*nQuadraturePoints_elementBoundary+
6624 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6628 f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6634 flux[ebNE*nQuadraturePoints_elementBoundary+
6636 n[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6640 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6651 int nElementBoundaries_element,
6652 int nQuadraturePoints_elementBoundary,
6654 int* exteriorElementBoundaries,
6655 int* elementBoundaryElements,
6656 int* elementBoundaryLocalElementBoundaries,
6657 int *isDOFBoundary_p,
6658 int *isDOFBoundary_u,
6659 int *isDOFBoundary_v,
6660 int *isDOFBoundary_w,
6687 int ebNE,ebN,eN_global,ebN_element,k;
6688 for(ebNE=0;ebNE<nExteriorElementBoundaries_global;ebNE++)
6690 ebN = exteriorElementBoundaries[ebNE];
6691 eN_global = elementBoundaryElements[ebN*2+0];
6692 ebN_element = elementBoundaryLocalElementBoundaries[ebN*2+0];
6693 for(k=0;k<nQuadraturePoints_elementBoundary;k++)
6695 fluxpu[ebN*nQuadraturePoints_elementBoundary+
6697 fluxpv[ebN*nQuadraturePoints_elementBoundary+
6699 fluxpw[ebN*nQuadraturePoints_elementBoundary+
6701 flux[ebN*nQuadraturePoints_elementBoundary+
6705 if (isDOFBoundary_p[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6707 dfluxpu_dp[ebN*nQuadraturePoints_elementBoundary+
6710 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6711 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6715 dfpu_dp[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6716 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6719 fluxpu[ebN*nQuadraturePoints_elementBoundary+
6721 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6722 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6726 fpu[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6727 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6730 dfluxpv_dp[ebN*nQuadraturePoints_elementBoundary+
6733 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6734 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6738 dfpv_dp[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6739 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6742 fluxpv[ebN*nQuadraturePoints_elementBoundary+
6744 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6745 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6749 fpv[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6750 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6753 dfluxpw_dp[ebN*nQuadraturePoints_elementBoundary+
6756 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6757 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6761 dfpw_dp[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6762 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6765 fluxpw[ebN*nQuadraturePoints_elementBoundary+
6767 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6768 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6772 fpw[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6773 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6779 fluxpu[ebN*nQuadraturePoints_elementBoundary+
6781 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6782 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6786 bc_fpu[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6789 fluxpv[ebN*nQuadraturePoints_elementBoundary+
6791 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6792 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6796 bc_fpv[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6799 fluxpw[ebN*nQuadraturePoints_elementBoundary+
6801 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6802 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6806 bc_fpw[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6811 if (isDOFBoundary_u[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6813 dflux_du[ebN*nQuadraturePoints_elementBoundary+
6816 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6817 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6821 df_du[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6822 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6825 flux[ebN*nQuadraturePoints_elementBoundary+
6827 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6828 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6832 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6833 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6839 flux[ebN*nQuadraturePoints_elementBoundary+
6841 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6842 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6846 bc_f[ebNE*nQuadraturePoints_elementBoundary*nSpace+
6850 if (isDOFBoundary_v[ebNE*nQuadraturePoints_elementBoundary+k] != 1)
6852 dflux_dv[ebN*nQuadraturePoints_elementBoundary+
6855 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6856 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6860 df_dv[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6861 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6864 flux[ebN*nQuadraturePoints_elementBoundary+
6866 n[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6867 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6871 f[eN_global*nElementBoundaries_element*nQuadraturePoints_elementBoundary*nSpace+
6872 ebN_element*nQuadraturePoints_elementBoundary*nSpace+
6878 flux[ebN*nQuadraturePoints_elementBoun