41 template<
unsigned DIM>
44 Vector<double>& press_mass_diag,
45 Vector<double>& veloc_mass_diag,
46 const unsigned& which_one)
49 unsigned n_dof = ndof();
52 if ((which_one == 0) || (which_one == 1))
55 press_mass_diag.assign(n_dof, 0.0);
59 if ((which_one == 0) || (which_one == 2))
62 veloc_mass_diag.assign(n_dof, 0.0);
66 HangInfo* hang_info_pt = 0;
69 unsigned n_master = 1;
72 double hang_weight = 1.0;
75 unsigned n_node = nnode();
78 Vector<double>
s(DIM + 1, 0.0);
81 Vector<unsigned> u_nodal_index(DIM, 0.0);
84 for (
unsigned i = 0;
i < DIM;
i++)
87 u_nodal_index[
i] = this->u_index_nst(
i);
94 unsigned n_pres = this->npres_nst();
101 int p_index = this->p_nodal_index_nst();
105 bool pressure_dof_is_hanging[n_pres];
111 for (
unsigned l = 0; l < n_pres; l++)
114 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
121 for (
unsigned l = 0; l < n_pres; l++)
124 pressure_dof_is_hanging[l] =
false;
129 unsigned n_intpt = integral_pt()->nweight();
135 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
138 double w = integral_pt()->weight(ipt);
141 double J = J_eulerian_at_knot(ipt);
144 for (
unsigned i = 0;
i < DIM + 1;
i++)
147 s[
i] = integral_pt()->knot(ipt,
i);
154 if ((which_one == 0) || (which_one == 2))
157 shape_at_knot(ipt, psi);
160 for (
unsigned l = 0; l < n_node; l++)
163 bool is_node_hanging = node_pt(l)->is_hanging();
169 hang_info_pt = node_pt(l)->hanging_pt();
172 n_master = hang_info_pt->nmaster();
182 for (
unsigned m = 0; m < n_master; m++)
185 for (
unsigned i = 0;
i < DIM;
i++)
191 local_eqn = this->local_hang_eqn(
192 hang_info_pt->master_node_pt(m), u_nodal_index[
i]);
195 hang_weight = hang_info_pt->master_weight(m);
201 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i]);
211 veloc_mass_diag[local_eqn] += pow(psi[l] * hang_weight, 2) *
W;
219 if ((which_one == 0) || (which_one == 1))
222 this->pshape_nst(
s, psi_p);
225 for (
unsigned l = 0; l < n_pres; l++)
228 if (pressure_dof_is_hanging[l])
231 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
234 n_master = hang_info_pt->nmaster();
244 for (
unsigned m = 0; m < n_master; m++)
247 if (pressure_dof_is_hanging[l])
251 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
254 hang_weight = hang_info_pt->master_weight(m);
260 local_eqn = this->p_local_eqn(l);
270 press_mass_diag[local_eqn] += pow(psi_p[l] * hang_weight, 2) *
W;
284 template<
unsigned DIM>
287 Vector<double>& residuals,
288 DenseMatrix<double>& jacobian,
289 const unsigned& flag)
291 OomphLibWarning(
"I'm not sure this is correct yet...",
292 OOMPH_CURRENT_FUNCTION,
293 OOMPH_EXCEPTION_LOCATION);
296 if (ndof() == 0)
return;
299 unsigned n_node = nnode();
302 unsigned n_pres = this->npres_nst();
305 unsigned u_nodal_index[DIM];
306 for (
unsigned i = 0;
i < DIM;
i++)
308 u_nodal_index[
i] = this->u_index_nst(
i);
314 int p_index = this->p_nodal_index_nst();
318 bool pressure_dof_is_hanging[n_pres];
323 for (
unsigned l = 0; l < n_pres; ++l)
325 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
333 "Pressure advection diffusion does not work in this case\n",
334 OOMPH_CURRENT_FUNCTION,
335 OOMPH_EXCEPTION_LOCATION);
337 for (
unsigned l = 0; l < n_pres; ++l)
339 pressure_dof_is_hanging[l] =
false;
345 DShape dpsidx(n_node, DIM);
348 Shape psip(n_pres), testp(n_pres);
349 DShape dpsip(n_pres, DIM);
350 DShape dtestp(n_pres, DIM);
353 unsigned n_intpt = integral_pt()->nweight();
356 Vector<double>
s(DIM);
360 double scaled_re = this->re() * this->density_ratio();
363 int local_eqn = 0, local_unknown = 0;
366 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
369 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
372 for (
unsigned i = 0;
i < DIM;
i++)
s[
i] = integral_pt()->knot(ipt,
i);
375 double w = integral_pt()->weight(ipt);
379 double J = this->dshape_eulerian_at_knot(ipt, psif, dpsidx);
382 this->dpshape_and_dptest_eulerian_nst(
s, psip, dpsip, testp, dtestp);
389 Vector<double> interpolated_u(DIM, 0.0);
390 Vector<double> interpolated_x(DIM, 0.0);
391 Vector<double> interpolated_dpdx(DIM, 0.0);
394 for (
unsigned l = 0; l < n_pres; l++)
396 for (
unsigned i = 0;
i < DIM;
i++)
398 interpolated_dpdx[
i] += this->p_nst(l) * dpsip(l,
i);
405 for (
unsigned l = 0; l < n_node; l++)
408 for (
unsigned i = 0;
i < DIM;
i++)
411 double u_value = nodal_value(l, u_nodal_index[
i]);
412 interpolated_u[
i] += u_value * psif[l];
413 interpolated_x[
i] += nodal_position(l,
i) * psif[l];
419 if (this->Press_adv_diff_source_fct_pt != 0)
421 source = this->Press_adv_diff_source_fct_pt(interpolated_x);
426 unsigned n_master = 1;
427 double hang_weight = 1.0;
431 for (
unsigned l = 0; l < n_pres; l++)
434 if (pressure_dof_is_hanging[l])
438 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
441 n_master = hang_info_pt->nmaster();
450 for (
unsigned m = 0; m < n_master; m++)
454 if (pressure_dof_is_hanging[l])
458 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
460 hang_weight = hang_info_pt->master_weight(m);
464 local_eqn = this->p_local_eqn(l);
471 residuals[local_eqn] -= source * testp[l] *
W * hang_weight;
472 for (
unsigned k = 0; k < DIM; k++)
474 residuals[local_eqn] +=
475 interpolated_dpdx[k] *
476 (scaled_re * interpolated_u[k] * testp[l] + dtestp(l, k)) *
W *
484 unsigned n_master2 = 1;
485 double hang_weight2 = 1.0;
488 for (
unsigned l2 = 0; l2 < n_pres; l2++)
491 if (pressure_dof_is_hanging[l2])
494 this->pressure_node_pt(l2)->hanging_pt(p_index);
497 n_master2 = hang_info2_pt->nmaster();
506 for (
unsigned m2 = 0; m2 < n_master2; m2++)
510 if (pressure_dof_is_hanging[l2])
513 local_unknown = this->local_hang_eqn(
514 hang_info2_pt->master_node_pt(m2), p_index);
516 hang_weight2 = hang_info2_pt->master_weight(m2);
520 local_unknown = this->p_local_eqn(l2);
525 if (local_unknown >= 0)
527 if ((
int(eqn_number(local_eqn)) !=
528 this->Pinned_fp_pressure_eqn) &&
529 (
int(eqn_number(local_unknown)) !=
530 this->Pinned_fp_pressure_eqn))
532 for (
unsigned k = 0; k < DIM; k++)
534 jacobian(local_eqn, local_unknown) +=
536 (scaled_re * interpolated_u[k] * testp[l] +
538 W * hang_weight * hang_weight2;
543 if ((
int(eqn_number(local_eqn)) ==
544 this->Pinned_fp_pressure_eqn) &&
545 (int(eqn_number(local_unknown)) ==
546 this->Pinned_fp_pressure_eqn))
548 jacobian(local_eqn, local_unknown) = 1.0;
562 this->Pressure_advection_diffusion_robin_element_pt.size();
563 for (
unsigned e = 0;
e < nrobin;
e++)
565 this->Pressure_advection_diffusion_robin_element_pt[
e]
566 ->fill_in_generic_residual_contribution_fp_press_adv_diff_robin_bc(
567 residuals, jacobian, flag);
578 template<
unsigned DIM>
581 DenseMatrix<double>& jacobian,
582 DenseMatrix<double>& mass_matrix,
583 const unsigned& flag)
586 if (ndof() == 0)
return;
589 unsigned n_node = nnode();
592 unsigned n_pres = this->npres_nst();
595 unsigned u_nodal_index[DIM];
598 for (
unsigned i = 0;
i < DIM;
i++)
601 u_nodal_index[
i] = this->u_index_nst(
i);
606 int p_index = this->p_nodal_index_nst();
610 bool pressure_dof_is_hanging[n_pres];
616 for (
unsigned l = 0; l < n_pres; ++l)
619 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
626 for (
unsigned l = 0; l < n_pres; ++l)
629 pressure_dof_is_hanging[l] =
false;
640 DShape dpsifdx(n_node, DIM + 1);
643 DShape dtestfdx(n_node, DIM + 1);
652 unsigned n_intpt = integral_pt()->nweight();
655 Vector<double>
s(DIM + 1, 0.0);
661 double scaled_re = this->re() * this->density_ratio();
664 double scaled_re_st = this->re() * this->st() * this->density_ratio();
667 double scaled_dre_st_dst = this->re() * this->density_ratio();
670 double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
673 double visc_ratio = this->viscosity_ratio();
676 Vector<double> G = this->g();
682 int local_unknown = 0;
685 HangInfo* hang_info_pt = 0;
688 HangInfo* hang_info2_pt = 0;
691 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
694 for (
unsigned i = 0;
i < DIM + 1;
i++)
697 s[
i] = integral_pt()->knot(ipt,
i);
701 double w = integral_pt()->weight(ipt);
704 double J = this->dshape_and_dtest_eulerian_at_knot_nst(
705 ipt, psif, dpsifdx, testf, dtestfdx);
708 this->pshape_nst(
s, psip, testp);
714 double interpolated_t = 0.0;
717 double interpolated_p = 0.0;
720 Vector<double> interpolated_x(DIM, 0.0);
723 Vector<double> interpolated_u(DIM, 0.0);
726 Vector<double> interpolated_dudt(DIM, 0.0);
729 Vector<double> mesh_velocity(DIM, 0.0);
732 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
735 for (
unsigned l = 0; l < n_pres; l++)
738 interpolated_p += this->p_nst(l) * psip[l];
745 for (
unsigned l = 0; l < n_node; l++)
748 interpolated_t += this->nodal_position(l, DIM) * psif(l);
751 for (
unsigned i = 0;
i < DIM;
i++)
754 double u_value = this->nodal_value(l, u_nodal_index[
i]);
757 interpolated_u[
i] += u_value * psif[l];
760 interpolated_x[
i] += this->nodal_position(l,
i) * psif[l];
763 interpolated_dudt[
i] += u_value * dpsifdx(l, DIM);
766 for (
unsigned j = 0; j < DIM; j++)
769 interpolated_dudx(
i, j) += u_value * dpsifdx(l, j);
775 if (!(this->ALE_is_disabled))
778 for (
unsigned l = 0; l < n_node; l++)
781 for (
unsigned i = 0;
i < DIM;
i++)
784 mesh_velocity[
i] += this->nodal_position(l,
i) * dpsifdx(l, DIM);
790 Vector<double> body_force(DIM, 0.0);
793 this->get_body_force_nst(
794 interpolated_t, ipt,
s, interpolated_x, body_force);
797 double source = this->get_source_nst(interpolated_t, ipt, interpolated_x);
806 unsigned n_master = 1;
809 double hang_weight = 1.0;
812 for (
unsigned l = 0; l < n_node; l++)
815 bool is_node_hanging = node_pt(l)->is_hanging();
821 hang_info_pt = node_pt(l)->hanging_pt();
824 n_master = hang_info_pt->nmaster();
834 for (
unsigned m = 0; m < n_master; m++)
837 for (
unsigned i = 0;
i < DIM;
i++)
843 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
847 hang_weight = hang_info_pt->master_weight(m);
853 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i]);
863 residuals[local_eqn] +=
864 body_force[
i] * testf[l] *
W * hang_weight;
867 residuals[local_eqn] +=
868 scaled_re_inv_fr * testf[l] * G[
i] *
W * hang_weight;
871 residuals[local_eqn] +=
872 interpolated_p * dtestfdx(l,
i) *
W * hang_weight;
875 residuals[local_eqn] -= (scaled_re_st * interpolated_dudt[
i] *
876 testf[l] *
W * hang_weight);
879 if (!(this->ALE_is_disabled))
882 for (
unsigned k = 0; k < DIM; k++)
885 residuals[local_eqn] +=
886 (scaled_re_st * mesh_velocity[k] * interpolated_dudx(
i, k) *
887 testf[l] *
W * hang_weight);
892 for (
unsigned k = 0; k < DIM; k++)
895 residuals[local_eqn] -=
896 (scaled_re * interpolated_u[k] * interpolated_dudx(
i, k) *
897 testf[l] *
W * hang_weight);
901 for (
unsigned k = 0; k < DIM; k++)
907 residuals[local_eqn] -=
908 ((interpolated_dudx(
i, k) +
909 this->Gamma[
i] * interpolated_dudx(k,
i)) *
910 visc_ratio * dtestfdx(l, k) *
W * hang_weight);
920 unsigned n_master2 = 1;
923 double hang_weight2 = 1.0;
926 for (
unsigned l2 = 0; l2 < n_node; l2++)
929 bool is_node2_hanging = node_pt(l2)->is_hanging();
932 if (is_node2_hanging)
935 hang_info2_pt = node_pt(l2)->hanging_pt();
938 n_master2 = hang_info2_pt->nmaster();
948 for (
unsigned m2 = 0; m2 < n_master2; m2++)
951 for (
unsigned i2 = 0; i2 < DIM; i2++)
954 if (is_node2_hanging)
957 local_unknown = this->local_hang_eqn(
958 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
961 hang_weight2 = hang_info2_pt->master_weight(m2);
968 this->nodal_local_eqn(l2, u_nodal_index[i2]);
975 if (local_unknown >= 0)
978 jacobian(local_eqn, local_unknown) -=
979 (visc_ratio * this->Gamma[
i] * dpsifdx(l2,
i) *
980 dtestfdx(l, i2) *
W * hang_weight * hang_weight2);
983 jacobian(local_eqn, local_unknown) -=
984 (scaled_re * psif[l2] * interpolated_dudx(
i, i2) *
985 testf[l] *
W * hang_weight * hang_weight2);
997 mass_matrix(local_eqn, local_unknown) +=
998 (scaled_re_st * psif[l2] * testf[l] *
W *
999 hang_weight * hang_weight2);
1003 jacobian(local_eqn, local_unknown) -=
1004 (scaled_re_st * dpsifdx(l2, DIM) * testf[l] *
W *
1005 hang_weight * hang_weight2);
1008 for (
unsigned k = 0; k < DIM; k++)
1011 jacobian(local_eqn, local_unknown) -=
1012 (scaled_re * interpolated_u[k] * dpsifdx(l2, k) *
1013 testf[l] *
W * hang_weight * hang_weight2);
1017 if (!(this->ALE_is_disabled))
1020 for (
unsigned k = 0; k < DIM; k++)
1023 jacobian(local_eqn, local_unknown) +=
1024 (scaled_re_st * mesh_velocity[k] *
1025 dpsifdx(l2, k) * testf[l] *
W * hang_weight *
1031 for (
unsigned k = 0; k < DIM; k++)
1034 jacobian(local_eqn, local_unknown) -=
1035 (visc_ratio * dpsifdx(l2, k) * dtestfdx(l, k) *
1036 W * hang_weight * hang_weight2);
1045 for (
unsigned l2 = 0; l2 < n_pres; l2++)
1048 if (pressure_dof_is_hanging[l2])
1053 this->pressure_node_pt(l2)->hanging_pt(p_index);
1056 n_master2 = hang_info2_pt->nmaster();
1066 for (
unsigned m2 = 0; m2 < n_master2; m2++)
1069 if (pressure_dof_is_hanging[l2])
1072 local_unknown = this->local_hang_eqn(
1073 hang_info2_pt->master_node_pt(m2), p_index);
1076 hang_weight2 = hang_info2_pt->master_weight(m2);
1082 local_unknown = this->p_local_eqn(l2);
1089 if (local_unknown >= 0)
1091 jacobian(local_eqn, local_unknown) +=
1092 (psip[l2] * dtestfdx(l,
i) *
W * hang_weight *
1103 if (this->Strouhal_is_stored_as_external_data)
1106 unsigned data_index = 0;
1109 unsigned value_index = 0;
1114 this->external_local_eqn(data_index, value_index);
1117 if (local_unknown >= 0)
1120 jacobian(local_eqn, local_unknown) -=
1121 (scaled_dre_st_dst * interpolated_dudt[
i] * testf[l] *
W *
1125 if (!this->ALE_is_disabled)
1128 for (
unsigned k = 0; k < DIM; k++)
1131 jacobian(local_eqn, local_unknown) +=
1132 (scaled_dre_st_dst * mesh_velocity[k] *
1133 interpolated_dudx(
i, k) * testf[l] *
W *
1149 for (
unsigned l = 0; l < n_pres; l++)
1152 if (pressure_dof_is_hanging[l])
1155 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1158 n_master = hang_info_pt->nmaster();
1168 for (
unsigned m = 0; m < n_master; m++)
1171 if (pressure_dof_is_hanging[l])
1175 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1178 hang_weight = hang_info_pt->master_weight(m);
1184 local_eqn = this->p_local_eqn(l);
1194 residuals[local_eqn] -= source * testp[l] *
W * hang_weight;
1197 for (
unsigned k = 0; k < DIM; k++)
1200 residuals[local_eqn] +=
1201 interpolated_dudx(k, k) * testp[l] *
W * hang_weight;
1211 unsigned n_master2 = 1;
1214 double hang_weight2 = 1.0;
1217 for (
unsigned l2 = 0; l2 < n_node; l2++)
1220 bool is_node2_hanging = node_pt(l2)->is_hanging();
1223 if (is_node2_hanging)
1226 hang_info2_pt = node_pt(l2)->hanging_pt();
1229 n_master2 = hang_info2_pt->nmaster();
1239 for (
unsigned m2 = 0; m2 < n_master2; m2++)
1242 for (
unsigned i2 = 0; i2 < DIM; i2++)
1245 if (is_node2_hanging)
1248 local_unknown = this->local_hang_eqn(
1249 hang_info2_pt->master_node_pt(m2), u_nodal_index[i2]);
1252 hang_weight2 = hang_info2_pt->master_weight(m2);
1259 this->nodal_local_eqn(l2, u_nodal_index[i2]);
1266 if (local_unknown >= 0)
1269 jacobian(local_eqn, local_unknown) +=
1270 (dpsifdx(l2, i2) * testp[l] *
W * hang_weight *
1290 template<
unsigned DIM>
1291 void RefineableSpaceTimeNavierStokesEquations<
1292 DIM>::get_dresidual_dnodal_coordinates(RankThreeTensor<double>&
1293 dresidual_dnodal_coordinates)
1296 throw OomphLibError(
"Space-time update needs to be checked!",
1297 OOMPH_CURRENT_FUNCTION,
1298 OOMPH_EXCEPTION_LOCATION);
1307 const unsigned n_node = nnode();
1310 double time = node_pt(0)->time_stepper_pt()->time_pt()->time();
1313 const unsigned n_pres = this->npres_nst();
1316 unsigned u_nodal_index[DIM];
1317 for (
unsigned i = 0;
i < DIM;
i++)
1319 u_nodal_index[
i] = this->u_index_nst(
i);
1324 const int p_index = this->p_nodal_index_nst();
1328 bool pressure_dof_is_hanging[n_pres];
1334 for (
unsigned l = 0; l < n_pres; ++l)
1336 pressure_dof_is_hanging[l] = pressure_node_pt(l)->is_hanging(p_index);
1342 for (
unsigned l = 0; l < n_pres; ++l)
1344 pressure_dof_is_hanging[l] =
false;
1349 Shape psif(n_node), testf(n_node);
1350 DShape dpsifdx(n_node, DIM), dtestfdx(n_node, DIM);
1353 Shape psip(n_pres), testp(n_pres);
1356 const unsigned n_shape_controlling_node = nshape_controlling_nodes();
1359 RankFourTensor<double> d_dpsifdx_dX(
1360 DIM, n_shape_controlling_node, n_node, DIM);
1361 RankFourTensor<double> d_dtestfdx_dX(
1362 DIM, n_shape_controlling_node, n_node, DIM);
1365 DenseMatrix<double> dJ_dX(DIM, n_shape_controlling_node);
1368 RankFourTensor<double> d_dudx_dX(DIM, n_shape_controlling_node, DIM, DIM);
1373 DenseMatrix<double> d_U_dX(DIM, n_shape_controlling_node, 0.0);
1376 const unsigned n_intpt = integral_pt()->nweight();
1379 Vector<double>
s(DIM);
1383 double scaled_re = this->re() * this->density_ratio();
1384 double scaled_re_st = this->re_st() * this->density_ratio();
1385 double scaled_re_inv_fr = this->re_invfr() * this->density_ratio();
1386 double visc_ratio = this->viscosity_ratio();
1387 Vector<double> G = this->g();
1395 bool element_has_node_with_aux_node_update_fct =
false;
1397 std::map<Node*, unsigned> local_shape_controlling_node_lookup =
1398 shape_controlling_node_lookup();
1401 for (std::map<Node*, unsigned>::iterator it =
1402 local_shape_controlling_node_lookup.begin();
1403 it != local_shape_controlling_node_lookup.end();
1407 Node* nod_pt = it->first;
1410 unsigned q = it->second;
1413 if (nod_pt->has_auxiliary_node_update_fct_pt())
1415 element_has_node_with_aux_node_update_fct =
true;
1418 Vector<double> u_ref(DIM);
1419 for (
unsigned i = 0;
i < DIM;
i++)
1421 u_ref[
i] = *(nod_pt->value_pt(u_nodal_index[
i]));
1425 for (
unsigned p = 0; p < DIM; p++)
1428 double backup = nod_pt->x(p);
1432 nod_pt->x(p) += eps_fd;
1435 nod_pt->perform_auxiliary_node_update_fct();
1439 (*(nod_pt->value_pt(u_nodal_index[p])) - u_ref[p]) / eps_fd;
1442 nod_pt->x(p) = backup;
1445 nod_pt->perform_auxiliary_node_update_fct();
1454 HangInfo* hang_info_pt = 0;
1457 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
1460 for (
unsigned i = 0;
i < DIM;
i++)
1462 s[
i] = integral_pt()->knot(ipt,
i);
1466 const double w = integral_pt()->weight(ipt);
1470 this->dshape_and_dtest_eulerian_at_knot_nst(ipt,
1480 this->pshape_nst(
s, psip, testp);
1484 double interpolated_p = 0.0;
1485 Vector<double> interpolated_u(DIM, 0.0);
1486 Vector<double> interpolated_x(DIM, 0.0);
1487 Vector<double> mesh_velocity(DIM, 0.0);
1488 Vector<double> dudt(DIM, 0.0);
1489 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
1492 for (
unsigned l = 0; l < n_pres; l++)
1494 interpolated_p += this->p_nst(l) * psip[l];
1500 for (
unsigned l = 0; l < n_node; l++)
1503 for (
unsigned i = 0;
i < DIM;
i++)
1506 const double u_value = nodal_value(l, u_nodal_index[
i]);
1507 interpolated_u[
i] += u_value * psif[l];
1508 interpolated_x[
i] += nodal_position(l,
i) * psif[l];
1509 dudt[
i] += this->du_dt_nst(l,
i) * psif[l];
1512 for (
unsigned j = 0; j < DIM; j++)
1514 interpolated_dudx(
i, j) += u_value * dpsifdx(l, j);
1519 if (!this->ALE_is_disabled)
1522 for (
unsigned l = 0; l < n_node; l++)
1525 for (
unsigned i = 0;
i < DIM;
i++)
1527 mesh_velocity[
i] += this->dnodal_position_dt(l,
i) * psif[l];
1535 for (
unsigned q = 0; q < n_shape_controlling_node; q++)
1538 for (
unsigned p = 0; p < DIM; p++)
1540 for (
unsigned i = 0;
i < DIM;
i++)
1542 for (
unsigned k = 0; k < DIM; k++)
1545 for (
unsigned j = 0; j < n_node; j++)
1548 nodal_value(j, u_nodal_index[
i]) * d_dpsifdx_dX(p, q, j, k);
1550 d_dudx_dX(p, q,
i, k) = aux;
1558 const double pos_time_weight =
1559 node_pt(0)->position_time_stepper_pt()->weight(1, 0);
1560 const double val_time_weight =
1561 node_pt(0)->time_stepper_pt()->weight(1, 0);
1564 Vector<double> body_force(DIM);
1565 this->get_body_force_nst(time, ipt,
s, interpolated_x, body_force);
1568 const double source = this->get_source_nst(time, ipt, interpolated_x);
1571 DenseMatrix<double> d_body_force_dx(DIM, DIM, 0.0);
1572 this->get_body_force_gradient_nst(
1573 time, ipt,
s, interpolated_x, d_body_force_dx);
1576 Vector<double> source_gradient(DIM, 0.0);
1577 this->get_source_gradient_nst(time, ipt, interpolated_x, source_gradient);
1587 unsigned n_master = 1;
1588 double hang_weight = 1.0;
1591 for (
unsigned l = 0; l < n_node; l++)
1594 bool is_node_hanging = node_pt(l)->is_hanging();
1597 if (is_node_hanging)
1599 hang_info_pt = node_pt(l)->hanging_pt();
1602 n_master = hang_info_pt->nmaster();
1611 for (
unsigned m = 0; m < n_master; m++)
1614 for (
unsigned i = 0;
i < DIM;
i++)
1618 if (is_node_hanging)
1621 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
1624 hang_weight = hang_info_pt->master_weight(m);
1630 local_eqn = this->nodal_local_eqn(l, u_nodal_index[
i]);
1640 for (
unsigned p = 0; p < DIM; p++)
1643 for (
unsigned q = 0; q < n_shape_controlling_node; q++)
1649 double sum = body_force[
i] * testf[l];
1652 sum += scaled_re_inv_fr * testf[l] * G[
i];
1655 sum += interpolated_p * dtestfdx(l,
i);
1661 for (
unsigned k = 0; k < DIM; k++)
1664 (interpolated_dudx(
i, k) +
1665 this->Gamma[
i] * interpolated_dudx(k,
i)) *
1672 sum -= scaled_re_st * dudt[
i] * testf[l];
1675 for (
unsigned k = 0; k < DIM; k++)
1677 double tmp = scaled_re * interpolated_u[k];
1678 if (!this->ALE_is_disabled)
1680 tmp -= scaled_re_st * mesh_velocity[k];
1682 sum -= tmp * interpolated_dudx(
i, k) * testf[l];
1687 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1688 sum * dJ_dX(p, q) * w * hang_weight;
1694 sum = d_body_force_dx(
i, p) * psif(q) * testf(l);
1697 sum += interpolated_p * d_dtestfdx_dX(p, q, l,
i);
1700 for (
unsigned k = 0; k < DIM; k++)
1703 visc_ratio * ((interpolated_dudx(
i, k) +
1704 this->Gamma[
i] * interpolated_dudx(k,
i)) *
1705 d_dtestfdx_dX(p, q, l, k) +
1706 (d_dudx_dX(p, q,
i, k) +
1707 this->Gamma[
i] * d_dudx_dX(p, q, k,
i)) *
1712 for (
unsigned k = 0; k < DIM; k++)
1714 double tmp = scaled_re * interpolated_u[k];
1715 if (!this->ALE_is_disabled)
1717 tmp -= scaled_re_st * mesh_velocity[k];
1719 sum -= tmp * d_dudx_dX(p, q,
i, k) * testf(l);
1721 if (!this->ALE_is_disabled)
1723 sum += scaled_re_st * pos_time_weight * psif(q) *
1724 interpolated_dudx(
i, p) * testf(l);
1728 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1729 sum * J * w * hang_weight;
1737 if (element_has_node_with_aux_node_update_fct)
1740 for (
unsigned q_local = 0; q_local < n_node; q_local++)
1744 unsigned n_master2 = 1;
1745 double hang_weight2 = 1.0;
1746 HangInfo* hang_info2_pt = 0;
1749 bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1751 Node* actual_shape_controlling_node_pt = node_pt(q_local);
1754 if (is_node_hanging2)
1756 hang_info2_pt = node_pt(q_local)->hanging_pt();
1759 n_master2 = hang_info2_pt->nmaster();
1768 for (
unsigned mm = 0; mm < n_master2; mm++)
1770 if (is_node_hanging2)
1772 actual_shape_controlling_node_pt =
1773 hang_info2_pt->master_node_pt(mm);
1774 hang_weight2 = hang_info2_pt->master_weight(mm);
1778 unsigned q = local_shape_controlling_node_lookup
1779 [actual_shape_controlling_node_pt];
1782 for (
unsigned p = 0; p < DIM; p++)
1784 double sum = -visc_ratio * this->Gamma[
i] *
1785 dpsifdx(q_local,
i) * dtestfdx(l, p) -
1786 scaled_re * psif(q_local) *
1787 interpolated_dudx(
i, p) * testf(l);
1790 sum -= scaled_re_st * val_time_weight * psif(q_local) *
1792 for (
unsigned k = 0; k < DIM; k++)
1795 visc_ratio * dpsifdx(q_local, k) * dtestfdx(l, k);
1796 double tmp = scaled_re * interpolated_u[k];
1797 if (!this->ALE_is_disabled)
1799 tmp -= scaled_re_st * mesh_velocity[k];
1801 sum -= tmp * dpsifdx(q_local, k) * testf(l);
1805 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1806 sum * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1823 for (
unsigned l = 0; l < n_pres; l++)
1826 if (pressure_dof_is_hanging[l])
1830 hang_info_pt = this->pressure_node_pt(l)->hanging_pt(p_index);
1833 n_master = hang_info_pt->nmaster();
1842 for (
unsigned m = 0; m < n_master; m++)
1846 if (pressure_dof_is_hanging[l])
1850 this->local_hang_eqn(hang_info_pt->master_node_pt(m), p_index);
1852 hang_weight = hang_info_pt->master_weight(m);
1856 local_eqn = this->p_local_eqn(l);
1864 for (
unsigned p = 0; p < DIM; p++)
1867 for (
unsigned q = 0; q < n_shape_controlling_node; q++)
1873 double aux = -source;
1876 for (
unsigned k = 0; k < DIM; k++)
1878 aux += interpolated_dudx(k, k);
1882 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1883 aux * dJ_dX(p, q) * testp[l] * w * hang_weight;
1890 aux = -source_gradient[p] * psif(q);
1891 for (
unsigned k = 0; k < DIM; k++)
1893 aux += d_dudx_dX(p, q, k, k);
1896 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1897 aux * testp[l] * J * w * hang_weight;
1904 if (element_has_node_with_aux_node_update_fct)
1907 for (
unsigned q_local = 0; q_local < n_node; q_local++)
1911 unsigned n_master2 = 1;
1912 double hang_weight2 = 1.0;
1913 HangInfo* hang_info2_pt = 0;
1916 bool is_node_hanging2 = node_pt(q_local)->is_hanging();
1918 Node* actual_shape_controlling_node_pt = node_pt(q_local);
1921 if (is_node_hanging2)
1923 hang_info2_pt = node_pt(q_local)->hanging_pt();
1926 n_master2 = hang_info2_pt->nmaster();
1935 for (
unsigned mm = 0; mm < n_master2; mm++)
1937 if (is_node_hanging2)
1939 actual_shape_controlling_node_pt =
1940 hang_info2_pt->master_node_pt(mm);
1941 hang_weight2 = hang_info2_pt->master_weight(mm);
1945 unsigned q = local_shape_controlling_node_lookup
1946 [actual_shape_controlling_node_pt];
1949 for (
unsigned p = 0; p < DIM; p++)
1951 double aux = dpsifdx(q_local, p) * testp(l);
1952 dresidual_dnodal_coordinates(local_eqn, p, q) +=
1953 aux * d_U_dX(p, q) * J * w * hang_weight * hang_weight2;
1969 template class RefineableSpaceTimeNavierStokesEquations<2>;
1970 template class RefineableQTaylorHoodSpaceTimeElement<2>;
static double Default_fd_jacobian_step
Double used for the default finite difference step in elemental jacobian calculations.
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,...
void fill_in_generic_pressure_advection_diffusion_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &compute_jacobian_flag)
Compute the residuals for the associated pressure advection diffusion problem. Used by the Fp precond...
void fill_in_generic_residual_contribution_nst(Vector< double > &residuals, DenseMatrix< double > &jacobian, DenseMatrix< double > &mass_matrix, const unsigned &compute_jacobian_flag)
Add the elements contribution to elemental residual vector and/or Jacobian matrix....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...