41 template<
unsigned SPATIAL_DIM>
44 Vector<double>& residuals,
45 DenseMatrix<double>& jacobian,
49 unsigned n_node = nnode();
52 unsigned u_nodal_index = this->u_index_ust_heat();
61 DShape dpsidx(n_node, SPATIAL_DIM + 1);
64 DShape dtestdx(n_node, SPATIAL_DIM + 1);
67 unsigned n_intpt = integral_pt()->nweight();
70 Vector<double>
s(SPATIAL_DIM + 1, 0.0);
73 double alpha_local = this->alpha();
76 double beta_local = this->beta();
82 int local_unknown = 0;
85 HangInfo* hang_info_pt = 0;
88 HangInfo* hang_info2_pt = 0;
91 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
94 for (
unsigned i = 0;
i < SPATIAL_DIM + 1;
i++)
97 s[
i] = integral_pt()->knot(ipt,
i);
101 double w = integral_pt()->weight(ipt);
104 double J = this->dshape_and_dtest_eulerian_at_knot_ust_heat(
105 ipt, psi, dpsidx, test, dtestdx);
111 double interpolated_t = 0.0;
114 double interpolated_u = 0.0;
117 double interpolated_dudt = 0.0;
120 Vector<double> interpolated_x(SPATIAL_DIM, 0.0);
123 Vector<double> interpolated_dudx(SPATIAL_DIM + 1, 0.0);
126 Vector<double> mesh_velocity(SPATIAL_DIM, 0.0);
132 for (
unsigned l = 0; l < n_node; l++)
135 double u_value = nodal_value(l, u_nodal_index);
138 interpolated_u += u_value * psi(l);
141 for (
unsigned j = 0; j < SPATIAL_DIM; j++)
144 interpolated_x[j] += nodal_position(l, j) * psi(l);
147 interpolated_dudx[j] += u_value * dpsidx(l, j);
151 interpolated_t += nodal_position(l, SPATIAL_DIM) * psi(l);
154 interpolated_dudt += u_value * dpsidx(l, SPATIAL_DIM);
158 if (!(this->ALE_is_disabled))
161 for (
unsigned j = 0; j < SPATIAL_DIM; j++)
164 for (
unsigned l = 0; l < n_node; l++)
167 mesh_velocity[j] += nodal_position(l, j) * dpsidx(l, SPATIAL_DIM);
176 this->get_source_ust_heat(interpolated_t, ipt, interpolated_x, source);
182 for (
unsigned l = 0; l < n_node; l++)
185 unsigned n_master = 1;
188 double hang_weight = 1.0;
191 bool is_node_hanging = this->node_pt(l)->is_hanging();
198 hang_info_pt = this->node_pt(l)->hanging_pt();
201 n_master = hang_info_pt->nmaster();
211 for (
unsigned m = 0; m < n_master; m++)
217 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
221 hang_weight = hang_info_pt->master_weight(m);
227 local_eqn = this->nodal_local_eqn(l, u_nodal_index);
237 residuals[local_eqn] +=
238 ((source + alpha_local * interpolated_dudt) * test(l) *
W *
242 if (!(this->ALE_is_disabled))
245 for (
unsigned k = 0; k < SPATIAL_DIM; k++)
248 residuals[local_eqn] -=
249 (alpha_local * mesh_velocity[k] * interpolated_dudx[k] *
250 test(l) *
W * hang_weight);
255 for (
unsigned k = 0; k < SPATIAL_DIM; k++)
258 residuals[local_eqn] += (beta_local * interpolated_dudx[k] *
259 dtestdx(l, k) *
W * hang_weight);
269 unsigned n_master2 = 1;
272 double hang_weight2 = 1.0;
275 for (
unsigned l2 = 0; l2 < n_node; l2++)
278 bool is_node2_hanging = this->node_pt(l2)->is_hanging();
281 if (is_node2_hanging)
284 hang_info2_pt = this->node_pt(l2)->hanging_pt();
287 n_master2 = hang_info2_pt->nmaster();
298 for (
unsigned m2 = 0; m2 < n_master2; m2++)
301 if (is_node2_hanging)
304 local_unknown = this->local_hang_eqn(
305 hang_info2_pt->master_node_pt(m2), u_nodal_index);
308 hang_weight2 = hang_info2_pt->master_weight(m2);
314 local_unknown = this->nodal_local_eqn(l2, u_nodal_index);
321 if (local_unknown >= 0)
324 jacobian(local_eqn, local_unknown) +=
325 (alpha_local * test(l) * dpsidx(l2, SPATIAL_DIM) *
W *
326 hang_weight * hang_weight2);
329 for (
unsigned k = 0; k < SPATIAL_DIM; k++)
333 jacobian(local_eqn, local_unknown) +=
334 (beta_local * dpsidx(l2, k) * dtestdx(l, k) *
W *
335 hang_weight * hang_weight2);
339 if (!(this->ALE_is_disabled))
342 for (
unsigned k = 0; k < SPATIAL_DIM; k++)
345 jacobian(local_eqn, local_unknown) -=
346 (alpha_local * mesh_velocity[k] * dpsidx(l2, k) *
347 test(l) *
W * hang_weight * hang_weight2);
364 template class RefineableQUnsteadyHeatSpaceTimeElement<2, 2>;
365 template class RefineableQUnsteadyHeatSpaceTimeElement<2, 3>;
366 template class RefineableQUnsteadyHeatSpaceTimeElement<2, 4>;
virtual void fill_in_generic_residual_contribution_ust_heat(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=0: compute residu...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...