36 template<
unsigned DIM>
44 template<
unsigned DIM>
48 template<
unsigned DIM>
53 template<
unsigned DIM>
58 template<
unsigned DIM>
68 template<
unsigned DIM>
73 const unsigned& which_one)
76 unsigned n_dof = ndof();
79 if ((which_one == 0) || (which_one == 1))
82 press_mass_diag.assign(n_dof, 0.0);
86 if ((which_one == 0) || (which_one == 2))
89 veloc_mass_diag.assign(n_dof, 0.0);
93 unsigned n_node = nnode();
102 for (
unsigned i = 0;
i < DIM;
i++)
105 u_nodal_index[
i] = this->u_index_nst(
i);
112 unsigned n_pres = npres_nst();
118 unsigned n_intpt = integral_pt()->nweight();
124 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
127 double w = integral_pt()->weight(ipt);
130 double J = J_eulerian_at_knot(ipt);
133 for (
unsigned i = 0;
i < DIM + 1;
i++)
136 s[
i] = integral_pt()->knot(ipt,
i);
143 if ((which_one == 0) || (which_one == 2))
146 shape_at_knot(ipt, psi);
149 for (
unsigned l = 0; l < n_node; l++)
152 for (
unsigned i = 0;
i < DIM;
i++)
155 local_eqn = nodal_local_eqn(l, u_nodal_index[
i]);
161 veloc_mass_diag[local_eqn] += pow(psi[l], 2) *
W;
168 if ((which_one == 0) || (which_one == 1))
171 pshape_nst(
s, psi_p);
174 for (
unsigned l = 0; l < n_pres; l++)
177 local_eqn = p_local_eqn(l);
183 press_mass_diag[local_eqn] += pow(psi_p[l], 2) *
W;
194 template<
unsigned DIM>
198 norm.resize(DIM + 1, 0.0);
204 unsigned n_intpt = integral_pt()->nweight();
207 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
210 for (
unsigned i = 0;
i < DIM + 1;
i++)
212 s[
i] = integral_pt()->knot(ipt,
i);
216 double w = integral_pt()->weight(ipt);
219 double J = J_eulerian(
s);
225 for (
unsigned i = 0;
i < DIM;
i++)
228 norm[
i] += pow(interpolated_u_nst(
s,
i), 2) *
W;
232 norm[DIM] += pow(interpolated_p_nst(
s), 2) *
W;
242 template<
unsigned DIM>
244 std::ostream& outfile,
257 double interpolated_t = 0.0;
266 unsigned n_intpt = integral_pt()->nweight();
272 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
275 for (
unsigned i = 0;
i < DIM + 1;
i++)
277 s[
i] = integral_pt()->knot(ipt,
i);
281 for (
unsigned i = 0;
i < DIM;
i++)
284 x[
i] = interpolated_x(
s,
i);
288 interpolated_t = interpolated_x(
s, DIM);
291 (*exact_soln_pt)(interpolated_t, x, exact_soln);
294 double w = integral_pt()->weight(ipt);
297 double J = J_eulerian(
s);
303 for (
unsigned i = 0;
i < DIM;
i++)
306 norm += exact_soln[
i] * exact_soln[
i] *
W;
309 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
314 norm += pow(exact_soln[DIM], 2) *
W;
317 error += pow(exact_soln[DIM] - interpolated_p_nst(
s), 2) *
W;
328 outfile << tecplot_zone_string(n_plot);
331 unsigned num_plot_points = nplot_points(n_plot);
334 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
337 get_s_plot(i_plot, n_plot,
s);
340 for (
unsigned i = 0;
i < DIM;
i++)
343 x[
i] = interpolated_x(
s,
i);
347 interpolated_t = interpolated_x(
s, DIM);
350 (*exact_soln_pt)(interpolated_t, x, exact_soln);
353 for (
unsigned i = 0;
i < DIM;
i++)
356 outfile << x[
i] <<
" ";
360 outfile << interpolated_t <<
" ";
363 for (
unsigned i = 0;
i < DIM;
i++)
366 outfile << exact_soln[
i] - interpolated_u_nst(
s,
i) <<
" ";
370 outfile << exact_soln[DIM] - interpolated_p_nst(
s) <<
" ";
373 outfile << std::endl;
377 write_tecplot_zone_footer(outfile, n_plot);
386 template<
unsigned DIM>
388 std::ostream& outfile,
395 error.resize(DIM + 1, 0.0);
398 norm.resize(DIM + 1, 0.0);
401 double interpolated_t = 0.0;
410 unsigned n_intpt = integral_pt()->nweight();
416 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
419 for (
unsigned i = 0;
i < DIM + 1;
i++)
421 s[
i] = integral_pt()->knot(ipt,
i);
425 for (
unsigned i = 0;
i < DIM;
i++)
428 x[
i] = interpolated_x(
s,
i);
432 interpolated_t = interpolated_x(
s, DIM);
435 (*exact_soln_pt)(interpolated_t, x, exact_soln);
438 double w = integral_pt()->weight(ipt);
441 double J = J_eulerian(
s);
447 for (
unsigned i = 0;
i < DIM;
i++)
450 norm[
i] += exact_soln[
i] * exact_soln[
i] *
W;
453 error[
i] += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
458 norm[DIM] += pow(exact_soln[DIM], 2) *
W;
461 error[DIM] += pow(exact_soln[DIM] - interpolated_p_nst(
s), 2) *
W;
472 outfile << tecplot_zone_string(n_plot);
475 unsigned num_plot_points = nplot_points(n_plot);
478 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
481 get_s_plot(i_plot, n_plot,
s);
484 for (
unsigned i = 0;
i < DIM;
i++)
487 x[
i] = interpolated_x(
s,
i);
491 interpolated_t = interpolated_x(
s, DIM);
494 (*exact_soln_pt)(interpolated_t, x, exact_soln);
497 for (
unsigned i = 0;
i < DIM;
i++)
500 outfile << x[
i] <<
" ";
504 outfile << interpolated_t <<
" ";
507 for (
unsigned i = 0;
i < DIM;
i++)
510 outfile << exact_soln[
i] - interpolated_u_nst(
s,
i) <<
" ";
514 outfile << exact_soln[DIM] - interpolated_p_nst(
s) <<
" ";
517 outfile << std::endl;
521 write_tecplot_zone_footer(outfile, n_plot);
530 template<
unsigned DIM>
544 double interpolated_t = 0.0;
553 unsigned n_intpt = integral_pt()->nweight();
559 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
562 for (
unsigned i = 0;
i < DIM + 1;
i++)
564 s[
i] = integral_pt()->knot(ipt,
i);
568 for (
unsigned i = 0;
i < DIM;
i++)
571 x[
i] = interpolated_x(
s,
i);
575 interpolated_t = interpolated_x(
s, DIM);
578 (*exact_soln_pt)(interpolated_t, x, exact_soln);
581 double w = integral_pt()->weight(ipt);
584 double J = J_eulerian(
s);
590 for (
unsigned i = 0;
i < DIM;
i++)
593 norm += exact_soln[
i] * exact_soln[
i] *
W;
596 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
606 template<
unsigned DIM>
625 unsigned n_intpt = integral_pt()->nweight();
631 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
634 for (
unsigned i = 0;
i < DIM + 1;
i++)
636 s[
i] = integral_pt()->knot(ipt,
i);
640 for (
unsigned i = 0;
i < DIM;
i++)
643 x[
i] = interpolated_x(
s,
i);
647 (*exact_soln_pt)(x, exact_soln);
650 double w = integral_pt()->weight(ipt);
653 double J = J_eulerian(
s);
659 for (
unsigned i = 0;
i < DIM;
i++)
662 norm += exact_soln[
i] * exact_soln[
i] *
W;
665 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
676 template<
unsigned DIM>
678 std::ostream& outfile,
690 double interpolated_t = 0.0;
699 unsigned n_intpt = integral_pt()->nweight();
702 outfile <<
"ZONE" << std::endl;
708 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
711 for (
unsigned i = 0;
i < DIM + 1;
i++)
713 s[
i] = integral_pt()->knot(ipt,
i);
717 for (
unsigned i = 0;
i < DIM;
i++)
720 x[
i] = interpolated_x(
s,
i);
724 interpolated_t = interpolated_x(
s, DIM);
727 (*exact_soln_pt)(x, exact_soln);
730 double w = integral_pt()->weight(ipt);
733 double J = J_eulerian(
s);
739 for (
unsigned i = 0;
i < DIM;
i++)
742 norm += exact_soln[
i] * exact_soln[
i] *
W;
745 error += pow(exact_soln[
i] - interpolated_u_nst(
s,
i), 2) *
W;
749 for (
unsigned i = 0;
i < DIM;
i++)
752 outfile << x[
i] <<
" ";
756 outfile << interpolated_t <<
" ";
759 for (
unsigned i = 0;
i < DIM;
i++)
762 outfile << exact_soln[
i] - interpolated_u_nst(
s,
i) <<
" ";
766 outfile << std::endl;
776 template<
unsigned DIM>
778 std::ostream& outfile,
779 const unsigned& n_plot,
783 double interpolated_t = 0.0;
792 outfile << tecplot_zone_string(n_plot);
798 unsigned num_plot_points = nplot_points(n_plot);
801 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
804 get_s_plot(i_plot, n_plot,
s);
807 for (
unsigned i = 0;
i < DIM;
i++)
810 x[
i] = interpolated_x(
s,
i);
814 interpolated_t = interpolated_x(
s, DIM);
817 (*exact_soln_pt)(x, exact_soln);
820 for (
unsigned i = 0;
i < DIM;
i++)
823 outfile << x[
i] <<
" ";
827 outfile << interpolated_t <<
" ";
830 for (
unsigned i = 0;
i < exact_soln.size();
i++)
833 outfile << exact_soln[
i] <<
" ";
837 outfile << std::endl;
841 write_tecplot_zone_footer(outfile, n_plot);
850 template<
unsigned DIM>
852 std::ostream& outfile,
853 const unsigned& n_plot,
858 double interpolated_t = 0.0;
867 outfile << tecplot_zone_string(n_plot);
873 unsigned num_plot_points = nplot_points(n_plot);
876 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
879 get_s_plot(i_plot, n_plot,
s);
882 for (
unsigned i = 0;
i < DIM;
i++)
885 x[
i] = interpolated_x(
s,
i);
889 interpolated_t = interpolated_x(
s, DIM);
892 (*exact_soln_pt)(interpolated_t, x, exact_soln);
895 for (
unsigned i = 0;
i < DIM;
i++)
898 outfile << x[
i] <<
" ";
902 outfile << interpolated_t <<
" ";
905 for (
unsigned i = 0;
i < exact_soln.size();
i++)
908 outfile << exact_soln[
i] <<
" ";
912 outfile << std::endl;
916 write_tecplot_zone_footer(outfile, n_plot);
927 template<
unsigned DIM>
929 const unsigned& n_plot,
933 unsigned n_node = nnode();
942 outfile << tecplot_zone_string(n_plot);
945 unsigned num_plot_points = nplot_points(n_plot);
948 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
951 get_s_plot(i_plot, n_plot,
s);
954 for (
unsigned i = 0;
i < DIM + 1;
i++)
957 outfile << interpolated_x(
s,
i) <<
" ";
961 for (
unsigned i = 0;
i < DIM;
i++)
964 outfile << interpolated_u_nst(
s,
i) <<
" ";
968 outfile << std::endl;
972 write_tecplot_zone_footer(outfile, n_plot);
981 template<
unsigned DIM>
983 const unsigned& n_plot)
986 unsigned n_node = nnode();
995 outfile << tecplot_zone_string(n_plot);
998 unsigned num_plot_points = nplot_points(n_plot);
1001 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1004 get_s_plot(i_plot, n_plot,
s);
1007 for (
unsigned i = 0;
i < DIM + 1;
i++)
1010 outfile << interpolated_x(
s,
i) <<
" ";
1014 for (
unsigned i = 0;
i < DIM;
i++)
1017 outfile << interpolated_u_nst(
s,
i) <<
" ";
1021 outfile << interpolated_p_nst(
s) <<
" ";
1024 outfile << std::endl;
1028 outfile << std::endl;
1031 write_tecplot_zone_footer(outfile, n_plot);
1040 template<
unsigned DIM>
1042 const unsigned& n_plot)
1048 fprintf(file_pt,
"%s", tecplot_zone_string(n_plot).c_str());
1051 unsigned num_plot_points = nplot_points(n_plot);
1054 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1057 get_s_plot(i_plot, n_plot,
s);
1060 for (
unsigned i = 0;
i < DIM + 1;
i++)
1063 fprintf(file_pt,
"%g ", interpolated_x(
s,
i));
1067 for (
unsigned i = 0;
i < DIM;
i++)
1070 fprintf(file_pt,
"%g ", interpolated_u_nst(
s,
i));
1074 fprintf(file_pt,
"%g \n", interpolated_p_nst(
s));
1078 fprintf(file_pt,
"\n");
1081 write_tecplot_zone_footer(file_pt, n_plot);
1090 template<
unsigned DIM>
1092 const unsigned& n_plot)
1098 outfile << tecplot_zone_string(n_plot);
1101 unsigned n_node = nnode();
1107 DShape dpsifdx(n_node, DIM);
1110 unsigned num_plot_points = nplot_points(n_plot);
1113 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1116 get_s_plot(i_plot, n_plot,
s);
1119 dshape_eulerian(
s, psif, dpsifdx);
1137 for (
unsigned i = 0;
i < DIM;
i++)
1140 unsigned u_nodal_index = u_index_nst(
i);
1143 for (
unsigned l = 0; l < n_node; l++)
1146 double u_value = nodal_value(l, u_nodal_index);
1149 du_dt[
i] += u_value * dpsifdx(l, DIM);
1152 mesh_velocity[
i] += nodal_position(l,
i) * dpsifdx(l, DIM);
1155 for (
unsigned j = 0; j < DIM; j++)
1158 interpolated_dudx(
i, j) += u_value * dpsifdx(l, j);
1167 for (
unsigned i = 0;
i < DIM;
i++)
1170 du_dt_ALE[
i] = du_dt[
i];
1173 for (
unsigned k = 0; k < DIM; k++)
1176 du_dt_ALE[
i] -= mesh_velocity[k] * interpolated_dudx(
i, k);
1181 for (
unsigned i = 0;
i < DIM + 1;
i++)
1184 outfile << interpolated_x(
s,
i) <<
" ";
1188 for (
unsigned i = 0;
i < DIM;
i++)
1191 outfile << interpolated_u_nst(
s,
i) <<
" ";
1195 outfile << interpolated_p_nst(
s) <<
" ";
1198 for (
unsigned i = 0;
i < DIM;
i++)
1201 outfile << du_dt_ALE[
i] <<
" ";
1205 outfile << dissipation(
s) <<
" ";
1208 outfile << std::endl;
1212 write_tecplot_zone_footer(outfile, n_plot);
1222 template<
unsigned DIM>
1224 std::ostream& outfile,
const unsigned& n_plot)
1236 vorticity.resize(1);
1242 vorticity.resize(3);
1248 "Can't output vorticity in 1D - in fact, what ";
1249 error_message +=
"do you mean\nby the 1D Navier-Stokes equations?\n";
1251 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1255 outfile << tecplot_zone_string(n_plot);
1258 unsigned num_plot_points = nplot_points(n_plot);
1261 for (
unsigned i_plot = 0; i_plot < num_plot_points; i_plot++)
1264 get_s_plot(i_plot, n_plot,
s);
1267 for (
unsigned i = 0;
i < DIM + 1;
i++)
1270 outfile << interpolated_x(
s,
i) <<
" ";
1274 get_vorticity(
s, vorticity);
1280 outfile << vorticity[0];
1286 outfile << vorticity[0] <<
" " << vorticity[1] <<
" " << vorticity[2]
1291 outfile << std::endl;
1295 outfile << std::endl;
1298 write_tecplot_zone_footer(outfile, n_plot);
1304 template<
unsigned DIM>
1311 unsigned n_intpt = integral_pt()->nweight();
1317 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1320 for (
unsigned i = 0;
i < DIM + 1;
i++)
1323 s[
i] = integral_pt()->knot(ipt,
i);
1327 double w = integral_pt()->weight(ipt);
1330 double J = J_eulerian(
s);
1336 strain_rate(
s, strainrate);
1339 double local_diss = 0.0;
1342 for (
unsigned i = 0;
i < DIM;
i++)
1345 for (
unsigned j = 0; j < DIM; j++)
1348 local_diss += 2.0 * strainrate(
i, j) * strainrate(
i, j);
1353 diss += local_diss * w * J;
1365 template<
unsigned DIM>
1373 strain_rate(
s, strainrate);
1376 double press = interpolated_p_nst(
s);
1379 for (
unsigned i = 0;
i < DIM;
i++)
1382 traction[
i] = -press *
N[
i];
1385 for (
unsigned j = 0; j < DIM; j++)
1388 traction[
i] += 2.0 * strainrate(
i, j) *
N[j];
1399 template<
unsigned DIM>
1414 strain_rate(
s, strainrate);
1417 double press = interpolated_p_nst(
s);
1420 for (
unsigned i = 0;
i < DIM;
i++)
1423 traction_p[
i] = -press *
N[
i];
1426 for (
unsigned j = 0; j < DIM; j++)
1429 traction_visc[
i] += 2.0 * strainrate(
i, j) *
N[j];
1433 traction_visc_n[
i] = traction_visc[
i] *
N[
i];
1436 traction_visc_t[
i] = traction_visc[
i] - traction_visc_n[
i];
1443 template<
unsigned DIM>
1449 strain_rate(
s, strainrate);
1452 double local_diss = 0.0;
1453 for (
unsigned i = 0;
i < DIM;
i++)
1455 for (
unsigned j = 0; j < DIM; j++)
1457 local_diss += 2.0 * strainrate(
i, j) * strainrate(
i, j);
1467 template<
unsigned DIM>
1472 if ((strainrate.
ncol() != DIM) || (strainrate.
nrow() != DIM))
1474 std::ostringstream error_message;
1475 error_message <<
"The strain rate has incorrect dimensions "
1476 << strainrate.
ncol() <<
" x " << strainrate.
nrow()
1477 <<
" not " << DIM << std::endl;
1479 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1487 unsigned n_node = nnode();
1493 DShape dpsidx(n_node, DIM + 1);
1496 dshape_eulerian(
s, psi, dpsidx);
1502 for (
unsigned i = 0;
i < DIM;
i++)
1505 unsigned u_nodal_index = u_index_nst(
i);
1508 for (
unsigned j = 0; j < DIM; j++)
1511 for (
unsigned l = 0; l < n_node; l++)
1514 dudx(
i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1520 for (
unsigned i = 0;
i < DIM;
i++)
1523 for (
unsigned j = 0; j < DIM; j++)
1526 strainrate(
i, j) = 0.5 * (dudx(
i, j) + dudx(j,
i));
1540 if (vorticity.size() != 1)
1542 std::ostringstream error_message;
1543 error_message <<
"Computation of vorticity in 2D requires a 1D vector\n"
1544 <<
"which contains the only non-zero component of the\n"
1545 <<
"vorticity vector. You've passed a vector of size "
1546 << vorticity.size() << std::endl;
1549 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1560 unsigned n_node = nnode();
1566 DShape dpsidx(n_node, DIM + 1);
1569 dshape_eulerian(
s, psi, dpsidx);
1575 for (
unsigned i = 0;
i < DIM;
i++)
1578 unsigned u_nodal_index = u_index_nst(
i);
1581 for (
unsigned j = 0; j < DIM; j++)
1584 for (
unsigned l = 0; l < n_node; l++)
1587 dudx(
i, j) += nodal_value(l, u_nodal_index) * dpsidx(l, j);
1593 vorticity[0] = dudx(1, 0) - dudx(0, 1);
1602 double& vorticity)
const
1608 get_vorticity(
s, vort);
1611 vorticity = vort[0];
1622 template<
unsigned DIM>
1626 double kin_en = 0.0;
1629 unsigned n_intpt = integral_pt()->nweight();
1635 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1638 for (
unsigned i = 0;
i < DIM + 1;
i++)
1640 s[
i] = integral_pt()->knot(ipt,
i);
1644 double w = integral_pt()->weight(ipt);
1647 double J = J_eulerian(
s);
1650 double veloc_squared = 0.0;
1651 for (
unsigned i = 0;
i < DIM;
i++)
1653 veloc_squared += interpolated_u_nst(
s,
i) * interpolated_u_nst(
s,
i);
1656 kin_en += 0.5 * veloc_squared * w * J;
1667 template<
unsigned DIM>
1671 double d_kin_en_dt = 0.0;
1674 const unsigned n_intpt = integral_pt()->nweight();
1677 const unsigned n_node = this->nnode();
1681 DShape dpsidx(n_node, DIM + 1);
1684 unsigned u_index[DIM];
1685 for (
unsigned i = 0;
i < DIM;
i++)
1687 u_index[
i] = this->u_index_nst(
i);
1691 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1694 double J = dshape_eulerian_at_knot(ipt, psi, dpsidx);
1697 double w = integral_pt()->weight(ipt);
1704 for (
unsigned l = 0; l < n_node; l++)
1707 for (
unsigned i = 0;
i < DIM;
i++)
1709 interpolated_u[
i] += nodal_value(l, u_index[
i]) * psi(l);
1710 interpolated_dudt[
i] += du_dt_nst(l, u_index[
i]) * psi(l);
1715 if (!ALE_is_disabled)
1721 for (
unsigned l = 0; l < n_node; l++)
1723 for (
unsigned i = 0;
i < DIM;
i++)
1725 mesh_velocity[
i] += this->dnodal_position_dt(l,
i) * psi(l);
1727 for (
unsigned j = 0; j < DIM; j++)
1729 interpolated_dudx(
i, j) +=
1730 this->nodal_value(l, u_index[
i]) * dpsidx(l, j);
1736 for (
unsigned i = 0;
i < DIM;
i++)
1738 for (
unsigned k = 0; k < DIM; k++)
1740 interpolated_dudt[
i] -= mesh_velocity[k] * interpolated_dudx(
i, k);
1748 for (
unsigned i = 0;
i < DIM;
i++)
1750 sum += interpolated_u[
i] * interpolated_dudt[
i];
1753 d_kin_en_dt += sum * w * J;
1764 template<
unsigned DIM>
1768 double press_int = 0;
1771 unsigned n_intpt = integral_pt()->nweight();
1777 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1780 for (
unsigned i = 0;
i < DIM + 1;
i++)
1782 s[
i] = integral_pt()->knot(ipt,
i);
1786 double w = integral_pt()->weight(ipt);
1789 double J = J_eulerian(
s);
1795 double press = interpolated_p_nst(
s);
1798 press_int += press *
W;
1810 template<
unsigned DIM>
1815 const unsigned& flag)
1818 OOMPH_CURRENT_FUNCTION,
1819 OOMPH_EXCEPTION_LOCATION);
1822 if (ndof() == 0)
return;
1825 unsigned n_node = nnode();
1828 unsigned n_pres = npres_nst();
1831 unsigned u_nodal_index[DIM];
1832 for (
unsigned i = 0;
i < DIM;
i++)
1834 u_nodal_index[
i] = u_index_nst(
i);
1839 DShape dpsidx(n_node, DIM + 1);
1843 Shape testp(n_pres);
1844 DShape dpsip(n_pres, DIM);
1845 DShape dtestp(n_pres, DIM);
1848 unsigned n_intpt = integral_pt()->nweight();
1855 double scaled_re = re() * density_ratio();
1859 int local_unknown = 0;
1862 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1865 for (
unsigned i = 0;
i < DIM + 1;
i++)
1867 s[
i] = integral_pt()->knot(ipt,
i);
1871 double w = integral_pt()->weight(ipt);
1875 double J = this->dshape_eulerian_at_knot(ipt, psif, dpsidx);
1878 this->dpshape_and_dptest_eulerian_nst(
s, psip, dpsip, testp, dtestp);
1891 for (
unsigned l = 0; l < n_pres; l++)
1893 for (
unsigned i = 0;
i < DIM;
i++)
1895 interpolated_dpdx[
i] += p_nst(l) * dpsip(l,
i);
1900 for (
unsigned i = 0;
i < DIM;
i++)
1902 interpolated_u[
i] = interpolated_u_nst(
s,
i);
1906 for (
unsigned i = 0;
i < DIM + 1;
i++)
1908 x[
i] = interpolated_x(
s,
i);
1912 double source = 0.0;
1913 if (Press_adv_diff_source_fct_pt != 0)
1915 source = Press_adv_diff_source_fct_pt(x);
1919 for (
unsigned l = 0; l < n_pres; l++)
1921 local_eqn = p_local_eqn(l);
1926 residuals[local_eqn] -= source * testp[l] *
W;
1928 for (
unsigned k = 0; k < DIM; k++)
1930 residuals[local_eqn] +=
1931 interpolated_dpdx[k] *
1932 (scaled_re * interpolated_u[k] * testp[l] + dtestp(l, k)) *
W;
1939 for (
unsigned l2 = 0; l2 < n_pres; l2++)
1941 local_unknown = p_local_eqn(l2);
1944 if (local_unknown >= 0)
1946 if ((
int(eqn_number(local_eqn)) != Pinned_fp_pressure_eqn) &&
1947 (
int(eqn_number(local_unknown)) != Pinned_fp_pressure_eqn))
1949 for (
unsigned k = 0; k < DIM; k++)
1951 jacobian(local_eqn, local_unknown) +=
1953 (scaled_re * interpolated_u[k] * testp[l] +
1960 if ((
int(eqn_number(local_eqn)) == Pinned_fp_pressure_eqn) &&
1961 (int(eqn_number(local_unknown)) ==
1962 Pinned_fp_pressure_eqn))
1965 jacobian(local_eqn, local_unknown) = 1.0;
1977 unsigned nrobin = Pressure_advection_diffusion_robin_element_pt.size();
1978 for (
unsigned e = 0;
e < nrobin;
e++)
1980 Pressure_advection_diffusion_robin_element_pt[
e]
1981 ->fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
1982 residuals, jacobian, flag);
1992 template<
unsigned DIM>
1997 const unsigned& flag)
2000 if (ndof() == 0)
return;
2003 unsigned n_node = nnode();
2006 unsigned n_pres = npres_nst();
2009 unsigned u_nodal_index[DIM];
2012 for (
unsigned i = 0;
i < DIM;
i++)
2015 u_nodal_index[
i] = u_index_nst(
i);
2022 Shape testf(n_node);
2025 DShape dpsifdx(n_node, DIM + 1);
2028 DShape dtestfdx(n_node, DIM + 1);
2034 Shape testp(n_pres);
2037 unsigned n_intpt = integral_pt()->nweight();
2046 double scaled_re = re() * density_ratio();
2049 double scaled_re_st = re() * st() * density_ratio();
2052 double scaled_dre_st_dst = re() * density_ratio();
2055 double scaled_re_inv_fr = re_invfr() * density_ratio();
2058 double visc_ratio = viscosity_ratio();
2067 int local_unknown = 0;
2070 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
2073 for (
unsigned i = 0;
i < DIM + 1;
i++)
2076 s[
i] = integral_pt()->knot(ipt,
i);
2080 double w = integral_pt()->weight(ipt);
2083 double J = dshape_and_dtest_eulerian_at_knot_nst(
2084 ipt, psif, dpsifdx, testf, dtestfdx);
2087 pshape_nst(
s, psip, testp);
2093 double interpolated_t = 0.0;
2096 double interpolated_p = 0.0;
2114 for (
unsigned l = 0; l < n_pres; l++)
2117 interpolated_p += p_nst(l) * psip[l];
2124 for (
unsigned l = 0; l < n_node; l++)
2127 interpolated_t += raw_nodal_position(l, DIM) * psif(l);
2130 for (
unsigned i = 0;
i < DIM;
i++)
2133 double u_value = raw_nodal_value(l, u_nodal_index[
i]);
2136 interpolated_u[
i] += u_value * psif[l];
2139 interpolated_x[
i] += raw_nodal_position(l,
i) * psif[l];
2142 interpolated_dudt[
i] += u_value * dpsifdx(l, DIM);
2145 for (
unsigned j = 0; j < DIM; j++)
2148 interpolated_dudx(
i, j) += u_value * dpsifdx(l, j);
2154 if (!ALE_is_disabled)
2157 for (
unsigned l = 0; l < n_node; l++)
2160 for (
unsigned i = 0;
i < DIM;
i++)
2165 mesh_velocity[
i] += raw_nodal_position(l,
i) * dpsifdx(l, DIM);
2174 get_body_force_nst(interpolated_t, ipt,
s, interpolated_x, body_force);
2177 double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2186 for (
unsigned l = 0; l < n_node; l++)
2189 for (
unsigned i = 0;
i < DIM;
i++)
2192 local_eqn = nodal_local_eqn(l, u_nodal_index[
i]);
2198 residuals[local_eqn] += body_force[
i] * testf[l] *
W;
2201 residuals[local_eqn] += scaled_re_inv_fr * testf[l] * G[
i] *
W;
2204 residuals[local_eqn] += interpolated_p * dtestfdx(l,
i) *
W;
2207 residuals[local_eqn] -=
2208 scaled_re_st * interpolated_dudt[
i] * testf[l] *
W;
2211 if (!ALE_is_disabled)
2214 for (
unsigned k = 0; k < DIM; k++)
2217 residuals[local_eqn] +=
2218 (scaled_re_st * mesh_velocity[k] * interpolated_dudx(
i, k) *
2224 for (
unsigned k = 0; k < DIM; k++)
2227 residuals[local_eqn] -= (scaled_re * interpolated_u[k] *
2228 interpolated_dudx(
i, k) * testf[l] *
W);
2232 for (
unsigned k = 0; k < DIM; k++)
2238 residuals[local_eqn] -= ((interpolated_dudx(
i, k) +
2239 Gamma[
i] * interpolated_dudx(k,
i)) *
2240 visc_ratio * dtestfdx(l, k) *
W);
2250 for (
unsigned l2 = 0; l2 < n_node; l2++)
2253 for (
unsigned i2 = 0; i2 < DIM; i2++)
2256 local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2259 if (local_unknown >= 0)
2262 jacobian(local_eqn, local_unknown) -=
2263 (visc_ratio * Gamma[
i] * dpsifdx(l2,
i) *
2264 dtestfdx(l, i2) *
W);
2267 jacobian(local_eqn, local_unknown) -=
2268 (scaled_re * psif[l2] * interpolated_dudx(
i, i2) *
2281 mass_matrix(local_eqn, local_unknown) +=
2282 (scaled_re_st * psif[l2] * testf[l] *
W);
2286 jacobian(local_eqn, local_unknown) -=
2287 (scaled_re_st * dpsifdx(l2, DIM) * testf[l] *
W);
2290 if (!ALE_is_disabled)
2293 for (
unsigned k = 0; k < DIM; k++)
2296 jacobian(local_eqn, local_unknown) +=
2297 (scaled_re_st * mesh_velocity[k] * dpsifdx(l2, k) *
2303 for (
unsigned k = 0; k < DIM; k++)
2306 jacobian(local_eqn, local_unknown) -=
2307 (scaled_re * interpolated_u[k] * dpsifdx(l2, k) *
2312 for (
unsigned k = 0; k < DIM; k++)
2315 jacobian(local_eqn, local_unknown) -=
2316 (visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) *
W);
2324 for (
unsigned l2 = 0; l2 < n_pres; l2++)
2328 local_unknown = p_local_eqn(l2);
2331 if (local_unknown >= 0)
2334 jacobian(local_eqn, local_unknown) +=
2335 psip[l2] * dtestfdx(l,
i) *
W;
2344 if (Strouhal_is_stored_as_external_data)
2347 unsigned data_index = 0;
2350 unsigned value_index = 0;
2355 this->external_local_eqn(data_index, value_index);
2358 if (local_unknown >= 0)
2361 jacobian(local_eqn, local_unknown) -=
2362 (scaled_dre_st_dst * interpolated_dudt[
i] * testf[l] *
W);
2365 if (!ALE_is_disabled)
2368 for (
unsigned k = 0; k < DIM; k++)
2371 jacobian(local_eqn, local_unknown) +=
2372 (scaled_dre_st_dst * mesh_velocity[k] *
2373 interpolated_dudx(
i, k) * testf[l] *
W);
2387 for (
unsigned l = 0; l < n_pres; l++)
2390 local_eqn = p_local_eqn(l);
2396 residuals[local_eqn] -= source * testp[l] *
W;
2399 for (
unsigned k = 0; k < DIM; k++)
2402 residuals[local_eqn] += interpolated_dudx(k, k) * testp[l] *
W;
2412 for (
unsigned l2 = 0; l2 < n_node; l2++)
2415 for (
unsigned i2 = 0; i2 < DIM; i2++)
2418 local_unknown = nodal_local_eqn(l2, u_nodal_index[i2]);
2421 if (local_unknown >= 0)
2424 jacobian(local_eqn, local_unknown) +=
2425 dpsifdx(l2, i2) * testp[l] *
W;
2441 template<
unsigned DIM>
2444 double*
const& parameter_pt,
2448 const unsigned& flag)
2452 OOMPH_CURRENT_FUNCTION,
2453 OOMPH_EXCEPTION_LOCATION);
2460 template<
unsigned DIM>
2469 OOMPH_CURRENT_FUNCTION,
2470 OOMPH_EXCEPTION_LOCATION);
2480 template<
unsigned DIM>
2485 throw OomphLibError(
"Space-time update needs to be checked!",
2486 OOMPH_CURRENT_FUNCTION,
2487 OOMPH_EXCEPTION_LOCATION);
2496 const unsigned n_node = nnode();
2499 const unsigned n_pres = npres_nst();
2502 unsigned u_nodal_index[DIM];
2503 for (
unsigned i = 0;
i < DIM;
i++)
2505 u_nodal_index[
i] = u_index_nst(
i);
2510 Shape testf(n_node);
2511 DShape dpsifdx(n_node, DIM + 1);
2512 DShape dtestfdx(n_node, DIM + 1);
2515 Shape psip(n_pres), testp(n_pres);
2533 const unsigned n_intpt = integral_pt()->nweight();
2540 double scaled_re = re() * density_ratio();
2541 double scaled_re_st = re_st() * density_ratio();
2542 double scaled_re_inv_fr = re_invfr() * density_ratio();
2543 double visc_ratio = viscosity_ratio();
2552 bool element_has_node_with_aux_node_update_fct =
false;
2553 for (
unsigned q = 0; q < n_node; q++)
2555 Node* nod_pt = node_pt(q);
2560 element_has_node_with_aux_node_update_fct =
true;
2564 for (
unsigned i = 0;
i < DIM;
i++)
2566 u_ref[
i] = *(nod_pt->
value_pt(u_nodal_index[
i]));
2570 for (
unsigned p = 0; p < DIM; p++)
2573 double backup = nod_pt->
x(p);
2577 nod_pt->
x(p) += eps_fd;
2584 (*(nod_pt->
value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
2587 nod_pt->
x(p) = backup;
2599 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
2602 for (
unsigned i = 0;
i < DIM + 1;
i++)
2604 s[
i] = integral_pt()->knot(ipt,
i);
2608 const double w = integral_pt()->weight(ipt);
2611 const double J = dshape_and_dtest_eulerian_at_knot_nst(ipt,
2621 pshape_nst(
s, psip, testp);
2624 double interpolated_t = 0.0;
2627 double interpolated_p = 0.0;
2645 for (
unsigned l = 0; l < n_pres; l++)
2648 interpolated_p += p_nst(l) * psip[l];
2655 for (
unsigned l = 0; l < n_node; l++)
2658 interpolated_t += raw_nodal_position(l, DIM) * psif(l);
2661 for (
unsigned i = 0;
i < DIM;
i++)
2664 double u_value = raw_nodal_value(l, u_nodal_index[
i]);
2667 interpolated_u[
i] += u_value * psif[l];
2670 interpolated_x[
i] += raw_nodal_position(l,
i) * psif[l];
2673 interpolated_dudt[
i] += u_value * dpsifdx(l, DIM);
2676 for (
unsigned j = 0; j < DIM; j++)
2679 interpolated_dudx(
i, j) += u_value * dpsifdx(l, j);
2685 if (!ALE_is_disabled)
2688 for (
unsigned l = 0; l < n_node; l++)
2691 for (
unsigned i = 0;
i < DIM;
i++)
2694 mesh_velocity[
i] += raw_nodal_position(l,
i) * dpsifdx(l, DIM);
2701 for (
unsigned q = 0; q < n_node; q++)
2704 for (
unsigned p = 0; p < DIM; p++)
2706 for (
unsigned i = 0;
i < DIM;
i++)
2708 for (
unsigned k = 0; k < DIM; k++)
2712 for (
unsigned j = 0; j < n_node; j++)
2714 aux += raw_nodal_value(j, u_nodal_index[
i]) *
2715 d_dpsifdx_dX(p, q, j, k);
2719 d_dudx_dX(p, q,
i, k) = aux;
2729 get_body_force_nst(interpolated_t, ipt,
s, interpolated_x, body_force);
2732 const double source = get_source_nst(interpolated_t, ipt, interpolated_x);
2741 get_body_force_gradient_nst(
2742 interpolated_t, ipt,
s, interpolated_x, d_body_force_dx);
2748 get_source_gradient_nst(
2749 interpolated_t, ipt, interpolated_x, source_gradient);
2753 const double pos_time_weight =
2754 (node_pt(0)->position_time_stepper_pt()->weight(1, 0));
2758 const double val_time_weight =
2759 node_pt(0)->time_stepper_pt()->weight(1, 0);
2768 for (
unsigned l = 0; l < n_node; l++)
2771 for (
unsigned i = 0;
i < DIM;
i++)
2775 local_eqn = nodal_local_eqn(l, u_nodal_index[
i]);
2781 for (
unsigned p = 0; p < DIM; p++)
2784 for (
unsigned q = 0; q < n_node; q++)
2790 double sum = body_force[
i] * testf[l];
2793 sum += scaled_re_inv_fr * testf[l] * G[
i];
2796 sum += interpolated_p * dtestfdx(l,
i);
2799 for (
unsigned k = 0; k < DIM; k++)
2805 sum -= (visc_ratio *
2806 (interpolated_dudx(
i, k) +
2807 Gamma[
i] * interpolated_dudx(k,
i)) *
2812 sum -= scaled_re_st * interpolated_dudt[
i] * testf[l];
2815 for (
unsigned k = 0; k < DIM; k++)
2818 sum -= scaled_re * interpolated_u[k] *
2819 interpolated_dudx(
i, k) * testf[l];
2823 if (!ALE_is_disabled)
2826 for (
unsigned k = 0; k < DIM; k++)
2829 sum += scaled_re_st * mesh_velocity[k] *
2830 interpolated_dudx(
i, k) * testf[l];
2836 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2837 sum * dJ_dX(p, q) * w;
2843 sum = d_body_force_dx(
i, p) * psif(q) * testf(l);
2846 sum += interpolated_p * d_dtestfdx_dX(p, q, l,
i);
2849 for (
unsigned k = 0; k < DIM; k++)
2851 sum -= (visc_ratio * ((interpolated_dudx(
i, k) +
2852 Gamma[
i] * interpolated_dudx(k,
i)) *
2853 d_dtestfdx_dX(p, q, l, k) +
2854 (d_dudx_dX(p, q,
i, k) +
2855 Gamma[
i] * d_dudx_dX(p, q, k,
i)) *
2860 for (
unsigned k = 0; k < DIM; k++)
2862 double tmp = scaled_re * interpolated_u[k];
2863 if (!ALE_is_disabled)
2865 tmp -= scaled_re_st * mesh_velocity[k];
2867 sum -= tmp * d_dudx_dX(p, q,
i, k) * testf(l);
2870 if (!ALE_is_disabled)
2872 sum += scaled_re_st * pos_time_weight * psif(q) *
2873 interpolated_dudx(
i, p) * testf(l);
2877 dresidual_dnodal_coordinates(local_eqn, p, q) += sum * J * w;
2881 if (element_has_node_with_aux_node_update_fct)
2884 -visc_ratio * Gamma[
i] * dpsifdx(q,
i) * dtestfdx(l, p) -
2885 scaled_re * psif(q) * interpolated_dudx(
i, p) * testf(l);
2888 sum -= scaled_re_st * val_time_weight * psif(q) * testf(l);
2889 for (
unsigned k = 0; k < DIM; k++)
2891 sum -= visc_ratio * dpsifdx(q, k) * dtestfdx(l, k);
2892 double tmp = scaled_re * interpolated_u[k];
2893 if (!ALE_is_disabled)
2894 tmp -= scaled_re_st * mesh_velocity[k];
2895 sum -= tmp * dpsifdx(q, k) * testf(l);
2898 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2899 sum * d_U_dX(p, q) * J * w;
2912 for (
unsigned l = 0; l < n_pres; l++)
2914 local_eqn = p_local_eqn(l);
2920 for (
unsigned p = 0; p < DIM; p++)
2923 for (
unsigned q = 0; q < n_node; q++)
2929 double aux = -source;
2932 for (
unsigned k = 0; k < DIM; k++)
2934 aux += interpolated_dudx(k, k);
2938 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2939 aux * dJ_dX(p, q) * testp[l] * w;
2945 aux = -source_gradient[p] * psif(q);
2946 for (
unsigned k = 0; k < DIM; k++)
2948 aux += d_dudx_dX(p, q, k, k);
2951 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2952 aux * testp[l] * J * w;
2956 if (element_has_node_with_aux_node_update_fct)
2958 aux = dpsifdx(q, p) * testp(l);
2959 dresidual_dnodal_coordinates(local_eqn, p, q) +=
2960 aux * d_U_dX(p, q) * J * w;
2980 3, 2, 3, 2, 2, 2, 3, 2, 3, 3, 2, 3, 2, 2,
2981 2, 3, 2, 3, 3, 2, 3, 2, 2, 2, 3, 2, 3};
2986 0, 2, 6, 8, 9, 11, 15, 17, 18, 20, 24, 26};
2997 template<
unsigned DIM>
2999 std::set<std::pair<Data*, unsigned>>& paired_load_data)
3002 unsigned u_index[DIM];
3005 for (
unsigned i = 0;
i < DIM;
i++)
3008 u_index[
i] = this->u_index_nst(
i);
3012 unsigned n_node = this->nnode();
3015 for (
unsigned n = 0; n < n_node; n++)
3019 for (
unsigned i = 0;
i < DIM;
i++)
3022 paired_load_data.insert(std::make_pair(this->node_pt(n), u_index[
i]));
3027 this->identify_pressure_data(paired_load_data);
3039 template<
unsigned DIM>
3041 std::set<std::pair<Data*, unsigned>>& paired_pressure_data)
3044 unsigned p_index =
static_cast<unsigned>(this->p_nodal_index_nst());
3047 unsigned n_pres = npres_nst();
3050 for (
unsigned l = 0; l < n_pres; l++)
3054 paired_pressure_data.insert(
3055 std::make_pair(this->node_pt(Pconv[l]), p_index));
3067 template<
unsigned DIM>
3069 std::list<std::pair<unsigned long, unsigned>>& dof_lookup_list)
const
3072 unsigned n_node = this->nnode();
3075 std::pair<unsigned, unsigned> dof_lookup;
3078 for (
unsigned n = 0; n < n_node; n++)
3081 unsigned n_value = this->required_nvalue(n);
3084 for (
unsigned i_value = 0; i_value < n_value; i_value++)
3087 int local_eqn_number = this->nodal_local_eqn(n, i_value);
3092 if (local_eqn_number >= 0)
3096 dof_lookup.first = this->eqn_number(local_eqn_number);
3106 dof_lookup.second = 0;
3109 dof_lookup_list.push_front(dof_lookup);
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
double * value_pt(const unsigned &i) const
Return the pointer to the i-the stored value. Typically this is required when direct access to the st...
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 .
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
double & x(const unsigned &i)
Return the i-th nodal coordinate.
void perform_auxiliary_node_update_fct()
Execute auxiliary update function (if any) – this can be used to update any nodal values following th...
bool has_auxiliary_node_update_fct_pt()
Boolean to indicate if node has a pointer to and auxiliary update function.
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
/////////////////////////////////////////////////////////////////////////// /////////////////////////...
void identify_load_data(std::set< std::pair< Data *, unsigned >> &paired_load_data)
Add to the set paired_load_data pairs containing.
void identify_pressure_data(std::set< std::pair< Data *, unsigned >> &paired_pressure_data)
Add to the set paired_pressure_data pairs containing.
void get_dof_numbers_for_unknowns(std::list< std::pair< unsigned long, unsigned >> &dof_lookup_list) const
Create a list of pairs for all unknowns in this element, so that the first entry in each pair contain...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
////////////////////////////////////////////////////////////////// //////////////////////////////////...
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
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,...
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 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....
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....
void compute_norm(Vector< double > &norm)
Compute the vector norm of the FEM solution.
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....
static double Default_Physical_Constant_Value
Static default value for the physical constants (all initialised to zero)
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...
double kin_energy() const
Get integral of kinetic energy over element.
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...
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 ...
static int Pressure_not_stored_at_node
Static "magic" number that indicates that the pressure is not stored at a node.
void output(std::ostream &outfile)
Output function: x,y,t,u,v,p in tecplot format. The default number of plot points is five.
double d_kin_energy_dt() const
Get integral of time derivative of kinetic energy over element.
double dissipation() const
Return integral of dissipation over element.
static Vector< double > Default_Gravity_vector
Static default value for the gravity vector.
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 pressure_integral() const
Integral of pressure 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 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....
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates....
void get_vorticity(const Vector< double > &s, Vector< double > &vorticity) const
Compute the vorticity vector at local coordinate s.
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 double Default_Physical_Ratio_Value
Static default value for the physical ratios (all are initialised to one)
void output_vorticity(std::ostream &outfile, const unsigned &nplot)
Output function: x,y,t,omega in tecplot format. nplot points in each coordinate direction.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...