36 template<
unsigned DIM>
44 template<
unsigned DIM>
50 template<
unsigned DIM>
51 double SpaceTimeNavierStokesMixedOrderEquations<
55 template<
unsigned DIM>
56 double SpaceTimeNavierStokesMixedOrderEquations<
57 DIM>::Default_Physical_Ratio_Value = 1.0;
60 template<
unsigned DIM>
61 Vector<double> SpaceTimeNavierStokesMixedOrderEquations<
62 DIM>::Default_Gravity_vector(DIM, 0.0);
70 template<
unsigned DIM>
75 const unsigned& which_one)
78 unsigned n_dof = ndof();
81 if ((which_one == 0) || (which_one == 1))
84 press_mass_diag.assign(n_dof, 0.0);
88 if ((which_one == 0) || (which_one == 2))
91 veloc_mass_diag.assign(n_dof, 0.0);
95 unsigned n_node = nnode();
104 for (
unsigned i = 0;
i < DIM;
i++)
107 u_nodal_index[
i] = this->u_index_nst(
i);
114 unsigned n_pres = npres_nst();
120 unsigned n_intpt = integral_pt()->nweight();
126 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
129 double w = integral_pt()->weight(ipt);
132 double J = J_eulerian_at_knot(ipt);
135 for (
unsigned i = 0;
i < DIM + 1;
i++)
138 s[
i] = integral_pt()->knot(ipt,
i);
145 if ((which_one == 0) || (which_one == 2))
148 shape_at_knot(ipt, psi);
151 for (
unsigned l = 0; l < n_node; l++)
154 for (
unsigned i = 0;
i < DIM;
i++)
157 local_eqn = nodal_local_eqn(l, u_nodal_index[
i]);
163 veloc_mass_diag[local_eqn] += pow(psi[l], 2) *
W;
170 if ((which_one == 0) || (which_one == 1))
173 pshape_nst(
s, psi_p);
176 for (
unsigned l = 0; l < n_pres; l++)
179 local_eqn = p_local_eqn(l);
185 press_mass_diag[local_eqn] += pow(psi_p[l], 2) *
W;
196 template<
unsigned DIM>
201 norm.resize(DIM + 1, 0.0);
207 unsigned n_intpt = integral_pt()->nweight();
210 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
213 for (
unsigned i = 0;
i < DIM + 1;
i++)
215 s[
i] = integral_pt()->knot(ipt,
i);
219 double w = integral_pt()->weight(ipt);
222 double J = J_eulerian(
s);
228 for (
unsigned i = 0;
i < DIM;
i++)
231 norm[
i] += pow(interpolated_u_nst(
s,
i), 2) *
W;
235 norm[DIM] += pow(interpolated_p_nst(
s), 2) *
W;
245 template<
unsigned DIM>
247 std::ostream& outfile,
260 double interpolated_t = 0.0;
269 unsigned n_intpt = integral_pt()->nweight();
275 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
278 for (
unsigned i = 0;
i < DIM + 1;
i++)
280 s[
i] = integral_pt()->knot(ipt,
i);
284 for (
unsigned i = 0;
i < DIM;
i++)
287 x[
i] = interpolated_x(
s,
i);
291 interpolated_t = interpolated_x(
s, DIM);
294 (*exact_soln_pt)(interpolated_t, x, exact_soln);
297 double w = integral_pt()->weight(ipt);
300 double J = J_eulerian(
s);
306 for (
unsigned i = 0;
i < DIM;
i++)
309 norm += exact_soln[
i] * exact_soln[
i] *
W;
312 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
317 norm += pow(exact_soln[DIM], 2) *
W;
320 error += pow(exact_soln[DIM] - interpolated_p_nst(
s), 2) *
W;
331 outfile << tecplot_zone_string(n_plot);
334 unsigned num_plot_points = nplot_points(n_plot);
337 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
340 get_s_plot(i_plot, n_plot,
s);
343 for (
unsigned i = 0;
i < DIM;
i++)
346 x[
i] = interpolated_x(
s,
i);
350 interpolated_t = interpolated_x(
s, DIM);
353 (*exact_soln_pt)(interpolated_t, x, exact_soln);
356 for (
unsigned i = 0;
i < DIM;
i++)
359 outfile << x[
i] <<
" ";
363 outfile << interpolated_t <<
" ";
366 for (
unsigned i = 0;
i < DIM;
i++)
369 outfile << exact_soln[
i] - interpolated_u_nst(
s,
i) <<
" ";
373 outfile << exact_soln[DIM] - interpolated_p_nst(
s) <<
" ";
376 outfile << std::endl;
380 write_tecplot_zone_footer(outfile, n_plot);
389 template<
unsigned DIM>
391 std::ostream& outfile,
398 error.resize(DIM + 1, 0.0);
401 norm.resize(DIM + 1, 0.0);
404 double interpolated_t = 0.0;
413 unsigned n_intpt = integral_pt()->nweight();
419 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
422 for (
unsigned i = 0;
i < DIM + 1;
i++)
424 s[
i] = integral_pt()->knot(ipt,
i);
428 for (
unsigned i = 0;
i < DIM;
i++)
431 x[
i] = interpolated_x(
s,
i);
435 interpolated_t = interpolated_x(
s, DIM);
438 (*exact_soln_pt)(interpolated_t, x, exact_soln);
441 double w = integral_pt()->weight(ipt);
444 double J = J_eulerian(
s);
450 for (
unsigned i = 0;
i < DIM;
i++)
453 norm[
i] += exact_soln[
i] * exact_soln[
i] *
W;
456 error[
i] += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
461 norm[DIM] += pow(exact_soln[DIM], 2) *
W;
464 error[DIM] += pow(exact_soln[DIM] - interpolated_p_nst(
s), 2) *
W;
475 outfile << tecplot_zone_string(n_plot);
478 unsigned num_plot_points = nplot_points(n_plot);
481 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
484 get_s_plot(i_plot, n_plot,
s);
487 for (
unsigned i = 0;
i < DIM;
i++)
490 x[
i] = interpolated_x(
s,
i);
494 interpolated_t = interpolated_x(
s, DIM);
497 (*exact_soln_pt)(interpolated_t, x, exact_soln);
500 for (
unsigned i = 0;
i < DIM;
i++)
503 outfile << x[
i] <<
" ";
507 outfile << interpolated_t <<
" ";
510 for (
unsigned i = 0;
i < DIM;
i++)
513 outfile << exact_soln[
i] - interpolated_u_nst(
s,
i) <<
" ";
517 outfile << exact_soln[DIM] - interpolated_p_nst(
s) <<
" ";
520 outfile << std::endl;
524 write_tecplot_zone_footer(outfile, n_plot);
533 template<
unsigned DIM>
547 double interpolated_t = 0.0;
556 unsigned n_intpt = integral_pt()->nweight();
562 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
565 for (
unsigned i = 0;
i < DIM + 1;
i++)
567 s[
i] = integral_pt()->knot(ipt,
i);
571 for (
unsigned i = 0;
i < DIM;
i++)
574 x[
i] = interpolated_x(
s,
i);
578 interpolated_t = interpolated_x(
s, DIM);
581 (*exact_soln_pt)(interpolated_t, x, exact_soln);
584 double w = integral_pt()->weight(ipt);
587 double J = J_eulerian(
s);
593 for (
unsigned i = 0;
i < DIM;
i++)
596 norm += exact_soln[
i] * exact_soln[
i] *
W;
599 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
609 template<
unsigned DIM>
628 unsigned n_intpt = integral_pt()->nweight();
634 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
637 for (
unsigned i = 0;
i < DIM + 1;
i++)
639 s[
i] = integral_pt()->knot(ipt,
i);
643 for (
unsigned i = 0;
i < DIM;
i++)
646 x[
i] = interpolated_x(
s,
i);
650 (*exact_soln_pt)(x, exact_soln);
653 double w = integral_pt()->weight(ipt);
656 double J = J_eulerian(
s);
662 for (
unsigned i = 0;
i < DIM;
i++)
665 norm += exact_soln[
i] * exact_soln[
i] *
W;
668 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
679 template<
unsigned DIM>
681 std::ostream& outfile,
693 double interpolated_t = 0.0;
702 unsigned n_intpt = integral_pt()->nweight();
705 outfile <<
"ZONE" << std::endl;
711 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
714 for (
unsigned i = 0;
i < DIM + 1;
i++)
716 s[
i] = integral_pt()->knot(ipt,
i);
720 for (
unsigned i = 0;
i < DIM;
i++)
723 x[
i] = interpolated_x(
s,
i);
727 interpolated_t = interpolated_x(
s, DIM);
730 (*exact_soln_pt)(x, exact_soln);
733 double w = integral_pt()->weight(ipt);
736 double J = J_eulerian(
s);
742 for (
unsigned i = 0;
i < DIM;
i++)
745 norm += exact_soln[
i] * exact_soln[
i] *
W;
748 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
752 for (
unsigned i = 0;
i < DIM;
i++)
755 outfile << x[
i] <<
" ";
759 outfile << interpolated_t <<
" ";
762 for (
unsigned i = 0;
i < DIM;
i++)
765 outfile << exact_soln[
i] - interpolated_u_nst(
s,
i) <<
" ";
769 outfile << std::endl;
779 template<
unsigned DIM>
781 std::ostream& outfile,
782 const unsigned& n_plot,
786 double interpolated_t = 0.0;
795 outfile << tecplot_zone_string(n_plot);
801 unsigned num_plot_points = nplot_points(n_plot);
804 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
807 get_s_plot(i_plot, n_plot,
s);
810 for (
unsigned i = 0;
i < DIM;
i++)
813 x[
i] = interpolated_x(
s,
i);
817 interpolated_t = interpolated_x(
s, DIM);
820 (*exact_soln_pt)(x, exact_soln);
823 for (
unsigned i = 0;
i < DIM;
i++)
826 outfile << x[
i] <<
" ";
830 outfile << interpolated_t <<
" ";
833 for (
unsigned i = 0;
i < exact_soln.size();
i++)
836 outfile << exact_soln[
i] <<
" ";
840 outfile << std::endl;
844 write_tecplot_zone_footer(outfile, n_plot);
853 template<
unsigned DIM>
855 std::ostream& outfile,
856 const unsigned& n_plot,
861 double interpolated_t = 0.0;
870 outfile << tecplot_zone_string(n_plot);
876 unsigned num_plot_points = nplot_points(n_plot);
879 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
882 get_s_plot(i_plot, n_plot,
s);
885 for (
unsigned i = 0;
i < DIM;
i++)
888 x[
i] = interpolated_x(
s,
i);
892 interpolated_t = interpolated_x(
s, DIM);
895 (*exact_soln_pt)(interpolated_t, x, exact_soln);
898 for (
unsigned i = 0;
i < DIM;
i++)
901 outfile << x[
i] <<
" ";
905 outfile << interpolated_t <<
" ";
908 for (
unsigned i = 0;
i < exact_soln.size();
i++)
911 outfile << exact_soln[
i] <<
" ";
915 outfile << std::endl;
919 write_tecplot_zone_footer(outfile, n_plot);
930 template<
unsigned DIM>
932 std::ostream& outfile,
const unsigned& n_plot,
const unsigned&
t)
935 unsigned n_node = nnode();
944 outfile << tecplot_zone_string(n_plot);
947 unsigned num_plot_points = nplot_points(n_plot);
950 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
953 get_s_plot(i_plot, n_plot,
s);
956 for (
unsigned i = 0;
i < DIM + 1;
i++)
959 outfile << interpolated_x(
s,
i) <<
" ";
963 for (
unsigned i = 0;
i < DIM;
i++)
966 outfile << interpolated_u_nst(
s,
i) <<
" ";
970 outfile << std::endl;
974 write_tecplot_zone_footer(outfile, n_plot);
983 template<
unsigned DIM>
985 std::ostream& outfile,
const unsigned& n_plot)
988 unsigned n_node = nnode();
997 outfile << tecplot_zone_string(n_plot);
1000 unsigned num_plot_points = nplot_points(n_plot);
1003 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1006 get_s_plot(i_plot, n_plot,
s);
1009 for (
unsigned i = 0;
i < DIM + 1;
i++)
1012 outfile << interpolated_x(
s,
i) <<
" ";
1016 for (
unsigned i = 0;
i < DIM;
i++)
1019 outfile << interpolated_u_nst(
s,
i) <<
" ";
1023 outfile << interpolated_p_nst(
s) <<
" ";
1026 outfile << std::endl;
1030 outfile << std::endl;
1033 write_tecplot_zone_footer(outfile, n_plot);
1042 template<
unsigned DIM>
1044 FILE* file_pt,
const unsigned& n_plot)
1050 fprintf(file_pt,
"%s", tecplot_zone_string(n_plot).c_str());
1053 unsigned num_plot_points = nplot_points(n_plot);
1056 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1059 get_s_plot(i_plot, n_plot,
s);
1062 for (
unsigned i = 0;
i < DIM + 1;
i++)
1065 fprintf(file_pt,
"%g ", interpolated_x(
s,
i));
1069 for (
unsigned i = 0;
i < DIM;
i++)
1072 fprintf(file_pt,
"%g ", interpolated_u_nst(
s,
i));
1076 fprintf(file_pt,
"%g \n", interpolated_p_nst(
s));
1080 fprintf(file_pt,
"\n");
1083 write_tecplot_zone_footer(file_pt, n_plot);
1092 template<
unsigned DIM>
1094 std::ostream& outfile,
const unsigned& n_plot)
1100 outfile << tecplot_zone_string(n_plot);
1103 unsigned n_node = nnode();
1109 DShape dpsifdx(n_node, DIM);
1112 unsigned num_plot_points = nplot_points(n_plot);
1115 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1118 get_s_plot(i_plot, n_plot,
s);
1122 OOMPH_CURRENT_FUNCTION,
1123 OOMPH_EXCEPTION_LOCATION);
1126 dshape_eulerian(
s, psif, dpsifdx);
1144 for (
unsigned i = 0;
i < DIM;
i++)
1147 unsigned u_nodal_index = u_index_nst(
i);
1150 for (
unsigned l = 0; l < n_node; l++)
1153 double u_value = nodal_value(l, u_nodal_index);
1156 du_dt[
i] += u_value * dpsifdx(l, DIM);
1159 mesh_velocity[
i] += nodal_position(l,
i) * dpsifdx(l, DIM);
1162 for (
unsigned j = 0; j < DIM; j++)
1165 interpolated_dudx(
i, j) += u_value * dpsifdx(l, j);
1174 for (
unsigned i = 0;
i < DIM;
i++)
1177 du_dt_ALE[
i] = du_dt[
i];
1180 for (
unsigned k = 0; k < DIM; k++)
1183 du_dt_ALE[
i] -= mesh_velocity[k] * interpolated_dudx(
i, k);
1188 for (
unsigned i = 0;
i < DIM + 1;
i++)
1191 outfile << interpolated_x(
s,
i) <<
" ";
1195 for (
unsigned i = 0;
i < DIM;
i++)
1198 outfile << interpolated_u_nst(
s,
i) <<
" ";
1202 outfile << interpolated_p_nst(
s) <<
" ";
1205 for (
unsigned i = 0;
i < DIM;
i++)
1208 outfile << du_dt_ALE[
i] <<
" ";
1212 outfile << dissipation(
s) <<
" ";
1215 outfile << std::endl;
1219 write_tecplot_zone_footer(outfile, n_plot);
1229 template<
unsigned DIM>
1231 std::ostream& outfile,
const unsigned& n_plot)
1243 vorticity.resize(1);
1249 vorticity.resize(3);
1255 "Can't output vorticity in 1D - in fact, what ";
1256 error_message +=
"do you mean\nby the 1D Navier-Stokes equations?\n";
1258 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1262 outfile << tecplot_zone_string(n_plot);
1265 unsigned num_plot_points = nplot_points(n_plot);
1268 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1271 get_s_plot(i_plot, n_plot,
s);
1274 for (
unsigned i = 0;
i < DIM + 1;
i++)
1277 outfile << interpolated_x(
s,
i) <<
" ";
1281 get_vorticity(
s, vorticity);
1287 outfile << vorticity[0];
1293 outfile << vorticity[0] <<
" " << vorticity[1] <<
" " << vorticity[2]
1298 outfile << std::endl;
1302 outfile << std::endl;
1305 write_tecplot_zone_footer(outfile, n_plot);
1311 template<
unsigned DIM>
1318 unsigned n_intpt = integral_pt()->nweight();
1324 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1327 for (
unsigned i = 0;
i < DIM + 1;
i++)
1330 s[
i] = integral_pt()->knot(ipt,
i);
1334 double w = integral_pt()->weight(ipt);
1337 double J = J_eulerian(
s);
1343 strain_rate(
s, strainrate);
1346 double local_diss = 0.0;
1349 for (
unsigned i = 0;
i < DIM;
i++)
1352 for (
unsigned j = 0; j < DIM; j++)
1355 local_diss += 2.0 * strainrate(
i, j) * strainrate(
i, j);
1360 diss += local_diss * w * J;
1372 template<
unsigned DIM>
1380 strain_rate(
s, strainrate);
1383 double press = interpolated_p_nst(
s);
1386 for (
unsigned i = 0;
i < DIM;
i++)
1389 traction[
i] = -press *
N[
i];
1392 for (
unsigned j = 0; j < DIM; j++)
1395 traction[
i] += 2.0 * strainrate(
i, j) *
N[j];
1406 template<
unsigned DIM>
1421 strain_rate(
s, strainrate);
1424 double press = interpolated_p_nst(
s);
1427 for (
unsigned i = 0;
i < DIM;
i++)
1430 traction_p[
i] = -press *
N[
i];
1433 for (
unsigned j = 0; j < DIM; j++)
1436 traction_visc[
i] += 2.0 * strainrate(
i, j) *
N[j];
1440 traction_visc_n[
i] = traction_visc[
i] *
N[
i];
1443 traction_visc_t[
i] = traction_visc[
i] - traction_visc_n[
i];
1450 template<
unsigned DIM>
1456 strain_rate(
s, strainrate);
1459 double local_diss = 0.0;
1460 for (
unsigned i = 0;
i < DIM;
i++)
1462 for (
unsigned j = 0; j < DIM; j++)
1464 local_diss += 2.0 * strainrate(
i, j) * strainrate(
i, j);
1474 template<
unsigned DIM>
1479 if ((strainrate.
ncol() != DIM) || (strainrate.
nrow() != DIM))
1481 std::ostringstream error_message;
1482 error_message <<
"The strain rate has incorrect dimensions "
1483 << strainrate.
ncol() <<
" x " << strainrate.
nrow()
1484 <<
" not " << DIM << std::endl;
1486 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1494 unsigned n_node = nnode();
1500 DShape dpsidx(n_node, DIM + 1);
1503 dshape_eulerian(
s, psi, dpsidx);
1509 for (
unsigned i = 0;
i < DIM;
i++)
1512 unsigned u_nodal_index = u_index_nst(
i);
1515 for (
unsigned j = 0; j < DIM; j++)
1518 for (
unsigned l = 0; l < n_node; l++)
1521 dudx(
i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1527 for (
unsigned i = 0;
i < DIM;
i++)
1530 for (
unsigned j = 0; j < DIM; j++)
1533 strainrate(
i, j) = 0.5 * (dudx(
i, j) + dudx(j,
i));
1547 if (vorticity.size() != 1)
1550 std::ostringstream error_message;
1553 error_message <<
"Computation of vorticity in 2D requires a 1D vector\n"
1554 <<
"which contains the only non-zero component of the\n"
1555 <<
"vorticity vector. You've passed a vector of size "
1556 << vorticity.size() << std::endl;
1560 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1571 unsigned n_node = nnode();
1577 DShape dpsidx(n_node, dim + 1);
1580 dshape_eulerian(
s, psi, dpsidx);
1586 for (
unsigned i = 0;
i < dim;
i++)
1589 unsigned u_nodal_index = u_index_nst(
i);
1592 for (
unsigned j = 0; j < dim; j++)
1595 for (
unsigned l = 0; l < n_node; l++)
1598 dudx(
i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1604 vorticity[0] = dudx(1, 0) - dudx(0, 1);
1619 get_vorticity(
s, vort);
1622 vorticity = vort[0];
1633 template<
unsigned DIM>
1637 double kin_en = 0.0;
1640 unsigned n_intpt = integral_pt()->nweight();
1646 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1649 for (
unsigned i = 0;
i < DIM + 1;
i++)
1651 s[
i] = integral_pt()->knot(ipt,
i);
1655 double w = integral_pt()->weight(ipt);
1658 double J = J_eulerian(
s);
1661 double veloc_squared = 0.0;
1662 for (
unsigned i = 0;
i < DIM;
i++)
1664 veloc_squared += interpolated_u_nst(
s,
i) * interpolated_u_nst(
s,
i);
1667 kin_en += 0.5 * veloc_squared * w * J;
1678 template<
unsigned DIM>
1682 double d_kin_en_dt = 0.0;
1685 const unsigned n_intpt = integral_pt()->nweight();
1688 const unsigned n_node = this->nnode();
1692 DShape dpsidx(n_node, DIM + 1);
1695 unsigned u_index[DIM];
1696 for (
unsigned i = 0;
i < DIM;
i++)
1698 u_index[
i] = this->u_index_nst(
i);
1702 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1705 double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
1708 double w = integral_pt()->weight(ipt);
1715 for (
unsigned l = 0; l < n_node; l++)
1718 for (
unsigned i = 0;
i < DIM;
i++)
1720 interpolated_u[
i] += nodal_value(l, u_index[
i]) * psi(l);
1721 interpolated_dudt[
i] += du_dt_nst(l, u_index[
i]) * psi(l);
1726 if (!ALE_is_disabled)
1732 for (
unsigned l = 0; l < n_node; l++)
1734 for (
unsigned i = 0;
i < DIM;
i++)
1736 mesh_velocity[
i] += this->dnodal_position_dt(l,
i) * psi(l);
1738 for (
unsigned j = 0; j < DIM; j++)
1740 interpolated_dudx(
i, j) +=
1741 this->nodal_value(l, u_index[
i]) * dpsidx(l, j);
1747 for (
unsigned i = 0;
i < DIM;
i++)
1749 for (
unsigned k = 0; k < DIM; k++)
1751 interpolated_dudt[
i] -= mesh_velocity[k] * interpolated_dudx(
i, k);
1759 for (
unsigned i = 0;
i < DIM;
i++)
1761 sum += interpolated_u[
i] * interpolated_dudt[
i];
1764 d_kin_en_dt += sum * w * J;
1775 template<
unsigned DIM>
1780 double press_int = 0;
1783 unsigned n_intpt = integral_pt()->nweight();
1789 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1792 for (
unsigned i = 0;
i < DIM + 1;
i++)
1794 s[
i] = integral_pt()->knot(ipt,
i);
1798 double w = integral_pt()->weight(ipt);
1801 double J = J_eulerian(
s);
1807 double press = interpolated_p_nst(
s);
1810 press_int += press *
W;
1822 template<
unsigned DIM>
1827 const unsigned& flag)
1831 OOMPH_CURRENT_FUNCTION,
1832 OOMPH_EXCEPTION_LOCATION);
1841 template<
unsigned DIM>
1846 const unsigned& flag)
1849 if (ndof() == 0)
return;
1852 unsigned n_node = nnode();
1855 unsigned n_pres = npres_nst();
1858 unsigned u_nodal_index[DIM];
1861 for (
unsigned i = 0;
i < DIM;
i++)
1864 u_nodal_index[
i] = u_index_nst(
i);
1871 Shape testf(n_node);
1874 DShape dpsifdx(n_node, DIM + 1);
1877 DShape dtestfdx(n_node, DIM + 1);
1883 Shape testp(n_pres);
1886 unsigned n_intpt = integral_pt()->nweight();
1895 double scaled_re = re() * density_ratio();
1898 double scaled_re_st = re_st() * density_ratio();
1902 double scaled_dre_st_dre_st = density_ratio();
1905 double scaled_re_inv_fr = re_invfr() * density_ratio();
1908 double visc_ratio = viscosity_ratio();
1917 int local_unknown = 0;
1920 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1923 for (
unsigned i = 0;
i < DIM + 1;
i++)
1926 s[
i] = integral_pt()->knot(ipt,
i);
1930 double w = integral_pt()->weight(ipt);
1933 double J = dshape_and_dtest_eulerian_at_knot_nst(
1934 ipt, psif, dpsifdx, testf, dtestfdx);
1937 pshape_nst(
s, psip, testp);
1943 double interpolated_t = 0.0;
1946 double interpolated_p = 0.0;
1964 for (
unsigned l = 0; l < n_pres; l++)
1967 interpolated_p += p_nst(l) * psip[l];
1974 for (
unsigned l = 0; l < n_node; l++)
1977 interpolated_t += raw_nodal_position(l, DIM) * psif(l);
1980 for (
unsigned i = 0;
i < DIM;
i++)
1983 double u_value = raw_nodal_value(l, u_nodal_index[
i]);
1986 interpolated_u[
i] += u_value * psif[l];
1989 interpolated_x[
i] += raw_nodal_position(l,
i) * psif[l];
1992 interpolated_dudt[
i] += u_value * dpsifdx(l, DIM);
1995 for (
unsigned j = 0; j < DIM; j++)
1998 interpolated_dudx(
i, j) += u_value * dpsifdx(l, j);
2004 if (!ALE_is_disabled)
2007 for (
unsigned l = 0; l < n_node; l++)
2010 for (
unsigned i = 0;
i < DIM;
i++)
2013 mesh_velocity[
i] += raw_nodal_position(l,
i) * dpsifdx(l, DIM);
2022 get_body_force_nst(interpolated_t, ipt,
s, interpolated_x, body_force);
2025 double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2034 for (
unsigned l = 0; l < n_node; l++)
2037 for (
unsigned i = 0;
i < DIM;
i++)
2040 local_eqn = nodal_local_eqn(l, u_nodal_index[
i]);
2046 residuals[local_eqn] += body_force[
i] * testf[l] *
W;
2049 residuals[local_eqn] += scaled_re_inv_fr * testf[l] * G[
i] *
W;
2052 residuals[local_eqn] += interpolated_p * dtestfdx(l,
i) *
W;
2055 residuals[local_eqn] -=
2056 scaled_re_st * interpolated_dudt[
i] * testf[l] *
W;
2059 if (!ALE_is_disabled)
2062 for (
unsigned k = 0; k < DIM; k++)
2065 residuals[local_eqn] +=
2066 (scaled_re_st * mesh_velocity[k] * interpolated_dudx(
i, k) *
2072 for (
unsigned k = 0; k < DIM; k++)
2075 residuals[local_eqn] -= (scaled_re * interpolated_u[k] *
2076 interpolated_dudx(
i, k) * testf[l] *
W);
2080 for (
unsigned k = 0; k < DIM; k++)
2086 residuals[local_eqn] -= ((interpolated_dudx(
i, k) +
2087 Gamma[
i] * interpolated_dudx(k,
i)) *
2088 visc_ratio * dtestfdx(l, k) *
W);
2098 for (
unsigned l2 = 0; l2 < n_node; l2++)
2101 for (
unsigned i2 = 0; i2 < DIM; i2++)
2104 local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2107 if (local_unknown >= 0)
2110 jacobian(local_eqn, local_unknown) -=
2111 (visc_ratio * Gamma[
i] * dpsifdx(l2,
i) *
2112 dtestfdx(l, i2) *
W);
2115 jacobian(local_eqn, local_unknown) -=
2116 (scaled_re * psif[l2] * interpolated_dudx(
i, i2) *
2129 mass_matrix(local_eqn, local_unknown) +=
2130 (scaled_re_st * psif[l2] * testf[l] *
W);
2134 jacobian(local_eqn, local_unknown) -=
2135 (scaled_re_st * dpsifdx(l2, DIM) * testf[l] *
W);
2138 if (!ALE_is_disabled)
2141 for (
unsigned k = 0; k < DIM; k++)
2144 jacobian(local_eqn, local_unknown) +=
2145 (scaled_re_st * mesh_velocity[k] * dpsifdx(l2, k) *
2151 for (
unsigned k = 0; k < DIM; k++)
2154 jacobian(local_eqn, local_unknown) -=
2155 (scaled_re * interpolated_u[k] * dpsifdx(l2, k) *
2160 for (
unsigned k = 0; k < DIM; k++)
2163 jacobian(local_eqn, local_unknown) -=
2164 (visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) *
W);
2172 for (
unsigned l2 = 0; l2 < n_pres; l2++)
2176 local_unknown = p_local_eqn(l2);
2179 if (local_unknown >= 0)
2182 jacobian(local_eqn, local_unknown) +=
2183 psip[l2] * dtestfdx(l,
i) *
W;
2192 if (ReynoldsStrouhal_is_stored_as_external_data)
2195 unsigned data_index = 0;
2198 unsigned value_index = 0;
2203 this->external_local_eqn(data_index, value_index);
2206 if (local_unknown >= 0)
2209 jacobian(local_eqn, local_unknown) -=
2210 (scaled_dre_st_dre_st * interpolated_dudt[
i] * testf[l] *
2214 if (!ALE_is_disabled)
2217 for (
unsigned k = 0; k < DIM; k++)
2220 jacobian(local_eqn, local_unknown) +=
2221 (scaled_dre_st_dre_st * mesh_velocity[k] *
2222 interpolated_dudx(
i, k) * testf[l] *
W);
2236 for (
unsigned l = 0; l < n_pres; l++)
2239 local_eqn = p_local_eqn(l);
2245 residuals[local_eqn] -= source * testp[l] *
W;
2248 for (
unsigned k = 0; k < DIM; k++)
2251 residuals[local_eqn] += interpolated_dudx(k, k) * testp[l] *
W;
2261 for (
unsigned l2 = 0; l2 < n_node; l2++)
2264 for (
unsigned i2 = 0; i2 < DIM; i2++)
2267 local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2270 if (local_unknown >= 0)
2273 jacobian(local_eqn, local_unknown) +=
2274 dpsifdx(l2, i2) * testp[l] *
W;
2290 template<
unsigned DIM>
2293 double*
const& parameter_pt,
2297 const unsigned& flag)
2301 OOMPH_CURRENT_FUNCTION,
2302 OOMPH_EXCEPTION_LOCATION);
2309 template<
unsigned DIM>
2318 OOMPH_CURRENT_FUNCTION,
2319 OOMPH_EXCEPTION_LOCATION);
2329 template<
unsigned DIM>
2332 dresidual_dnodal_coordinates)
2335 throw OomphLibError(
"Space-time update needs to be checked!",
2336 OOMPH_CURRENT_FUNCTION,
2337 OOMPH_EXCEPTION_LOCATION);
2351 {3, 2, 3, 2, 2, 2, 3, 2, 3, 2, 2, 2, 2, 2,
2352 2, 2, 2, 2, 3, 2, 3, 2, 2, 2, 3, 2, 3};
2357 0, 2, 6, 8, 18, 20, 24, 26};
2368 template<
unsigned DIM>
2370 std::set<std::pair<Data*, unsigned>>& paired_load_data)
2373 unsigned u_index[DIM];
2376 for (
unsigned i = 0;
i < DIM;
i++)
2379 u_index[
i] = this->u_index_nst(
i);
2383 unsigned n_node = this->nnode();
2386 for (
unsigned n = 0; n < n_node; n++)
2390 for (
unsigned i = 0;
i < DIM;
i++)
2393 paired_load_data.insert(std::make_pair(this->node_pt(n), u_index[
i]));
2398 this->identify_pressure_data(paired_load_data);
2410 template<
unsigned DIM>
2412 std::set<std::pair<Data*, unsigned>>& paired_pressure_data)
2415 unsigned p_index =
static_cast<unsigned>(this->p_nodal_index_nst());
2418 unsigned n_pres = npres_nst();
2421 for (
unsigned l = 0; l < n_pres; l++)
2425 paired_pressure_data.insert(
2426 std::make_pair(this->node_pt(Pconv[l]), p_index));
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
unsigned long nrow() const
Return the number of rows of the matrix.
unsigned long ncol() const
Return the number of columns of the matrix.
void initialise(const T &val)
Initialize all values in the matrix to val.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
void(* UnsteadyExactSolutionFctPt)(const double &, const Vector< double > &, Vector< double > &)
Function pointer for function that computes Vector-valued time-dependent function as .
An OomphLibError object which should be thrown when an run-time error is encountered....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void identify_pressure_data(std::set< std::pair< Data *, unsigned >> &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
void identify_load_data(std::set< std::pair< Data *, unsigned >> &paired_load_data)
Add to the set paired_load_data pairs containing.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void fill_in_contribution_to_hessian_vector_products(Vector< double > const &Y, DenseMatrix< double > const &C, DenseMatrix< double > &product)
Compute the hessian tensor vector products required to perform continuation of bifurcations analytica...
virtual void fill_in_generic_dresidual_contribution_nst(double *const ¶meter_pt, Vector< double > &dres_dparam, DenseMatrix< double > &djac_dparam, DenseMatrix< double > &dmass_matrix_dparam, const unsigned &flag)
Compute the derivatives of the residuals for the Navier-Stokes equations with respect to a parameter ...
void full_output(std::ostream &outfile)
Full output function: x,y,[z],u,v,[w],p,du/dt,dv/dt,[dw/dt],dissipation in tecplot format....
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
void output_veloc(std::ostream &outfile, const unsigned &nplot, const unsigned &t)
Output function: x,y,t,u,v in tecplot format. Use n_plot points in each coordinate direction at times...
void output_fct(std::ostream &outfile, const unsigned &nplot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact solution specified via function pointer at a given number of plot points....
virtual void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
void get_pressure_and_velocity_mass_matrix_diagonal(Vector< double > &press_mass_diag, Vector< double > &veloc_mass_diag, const unsigned &which_one=0)
Compute the diagonal of the velocity/pressure mass matrices. If which one=0, both are computed,...
virtual void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &flag)
Compute the residuals for the Navier-Stokes equations. Flag=1 (or 0): do (or don't) compute the Jacob...
double kin_energy() const
Get integral of kinetic energy over element.
void compute_norm(Vector< double > &norm)
Compute the vector norm of the FEM solution.
void strain_rate(const Vector< double > &s, DenseMatrix< double > &strain_rate) const
Strain-rate tensor: 1/2 (du_i/dx_j+du_j/dx_i)
static Vector< double > Gamma
Vector to decide whether the stress-divergence form is used or not N.B. This needs to be public so th...
void output_vorticity(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,t,omega in tecplot format. nplot points in each coordinate direction.
double dissipation() const
Return integral of dissipation over element.
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
void output(std::ostream &outfile)
Output function: x,y,t,u,v,p in tecplot format. The default number of plot points is five.
void get_traction(const Vector< double > &s, const Vector< double > &N, Vector< double > &traction)
Compute traction (on the viscous scale) exerted onto the fluid at local coordinate s....
double pressure_integral() const
Integral of pressure over element.
void compute_error(std::ostream &outfile, FiniteElement::UnsteadyExactSolutionFctPt exact_soln_pt, const double &time, double &error, double &norm)
Validate against exact solution at given time Solution is provided via function pointer....
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double Default_Physical_Constant_Value
Default value for physical constants.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...