42 template<
unsigned DIM,
unsigned NNODE_1D>
54 template <
unsigned DIM>
61 const unsigned n_node = nnode();
64 Shape psi(n_node), test(n_node);
65 DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
68 const unsigned u_nodal_index = u_index_poisson();
71 const unsigned n_intpt = integral_pt()->nweight();
74 int local_eqn=0, local_unknown=0;
77 for(
unsigned ipt=0;ipt<n_intpt;ipt++)
80 double w = integral_pt()->weight(ipt);
83 double J = dshape_and_dtest_eulerian_at_knot_poisson(ipt,psi,dpsidx,
91 double interpolated_u=0.0;
98 for(
unsigned l=0;l<n_node;l++)
101 double u_value = raw_nodal_value(l,u_nodal_index);
102 interpolated_u += u_value*psi(l);
104 for(
unsigned j=0;j<DIM;j++)
106 interpolated_x[j] += raw_nodal_position(l,j)*psi(l);
107 interpolated_dudx[j] += u_value*dpsidx(l,j);
115 get_source_poisson(ipt,interpolated_x,source);
121 for(
unsigned l=0;l<n_node;l++)
124 local_eqn = nodal_local_eqn(l,u_nodal_index);
129 residuals[local_eqn] += source*test(l)*
W;
132 for(
unsigned k=0;k<DIM;k++)
134 residuals[local_eqn] += interpolated_dudx[k]*dtestdx(l,k)*
W;
142 for(
unsigned l2=0;l2<n_node;l2++)
144 local_unknown = nodal_local_eqn(l2,u_nodal_index);
146 if(local_unknown >= 0)
149 for(
unsigned i=0;
i<DIM;
i++)
151 jacobian(local_eqn,local_unknown)
152 += dpsidx(l2,
i)*dtestdx(l,
i)*
W;
171 template <
unsigned DIM>
174 dresidual_dnodal_coordinates)
177 const unsigned n_node = nnode();
180 Shape psi(n_node), test(n_node);
181 DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
198 const unsigned u_nodal_index = u_index_poisson();
201 const unsigned n_intpt = integral_pt()->nweight();
207 for(
unsigned ipt=0;ipt<n_intpt;ipt++)
210 double w = integral_pt()->weight(ipt);
215 const double J = dshape_and_dtest_eulerian_at_knot_poisson(
216 ipt,psi,dpsidx,d_dpsidx_dX,test,dtestdx,d_dtestdx_dX,dJ_dX);
226 for(
unsigned l=0;l<n_node;l++)
229 double u_value = raw_nodal_value(l,u_nodal_index);
232 for(
unsigned i=0;
i<DIM;
i++)
234 interpolated_x[
i] += raw_nodal_position(l,
i)*psi(l);
235 interpolated_dudx[
i] += u_value*dpsidx(l,
i);
240 for(
unsigned q=0;q<n_node;q++)
243 for(
unsigned p=0;p<DIM;p++)
245 for(
unsigned i=0;
i<DIM;
i++)
248 for(
unsigned j=0;j<n_node;j++)
250 aux += raw_nodal_value(j,u_nodal_index)*d_dpsidx_dX(p,q,j,
i);
252 d_dudx_dX(p,q,
i) = aux;
258 get_source_poisson(ipt,interpolated_x,source);
261 get_source_gradient_poisson(ipt,interpolated_x,d_source_dx);
267 for(
unsigned l=0;l<n_node;l++)
270 local_eqn = nodal_local_eqn(l,u_nodal_index);
276 for(
unsigned p=0;p<DIM;p++)
279 for(
unsigned q=0;q<n_node;q++)
281 double sum = source*test(l)*dJ_dX(p,q)
282 + d_source_dx[p]*test(l)*psi(q)*J;
284 for(
unsigned i=0;
i<DIM;
i++)
286 sum += interpolated_dudx[
i]*(dtestdx(l,
i)*dJ_dX(p,q) +
287 d_dtestdx_dX(p,q,l,
i)*J)
288 + d_dudx_dX(p,q,
i)*dtestdx(l,
i)*J;
292 dresidual_dnodal_coordinates(local_eqn,p,q) += sum*w;
305 template <
unsigned DIM>
338 template <
unsigned DIM>
340 const unsigned &nplot)
347 outfile << tecplot_zone_string(nplot);
350 unsigned num_plot_points=nplot_points(nplot);
351 for (
unsigned iplot=0;iplot<num_plot_points;iplot++)
355 get_s_plot(iplot,nplot,s);
357 for(
unsigned i=0;
i<DIM;
i++)
359 outfile << interpolated_x(s,
i) <<
" ";
361 outfile << interpolated_u_poisson(s) << std::endl;
366 write_tecplot_zone_footer(outfile,nplot);
378 template <
unsigned DIM>
380 const unsigned &nplot)
386 fprintf(file_pt,
"%s",tecplot_zone_string(nplot).c_str());
389 unsigned num_plot_points=nplot_points(nplot);
390 for (
unsigned iplot=0;iplot<num_plot_points;iplot++)
393 get_s_plot(iplot,nplot,s);
395 for(
unsigned i=0;
i<DIM;
i++)
397 fprintf(file_pt,
"%g ",interpolated_x(s,
i));
399 fprintf(file_pt,
"%g \n",interpolated_u_poisson(s));
403 write_tecplot_zone_footer(file_pt,nplot);
416 template <
unsigned DIM>
418 const unsigned &nplot,
428 outfile << tecplot_zone_string(nplot);
434 unsigned num_plot_points=nplot_points(nplot);
435 for (
unsigned iplot=0;iplot<num_plot_points;iplot++)
439 get_s_plot(iplot,nplot,s);
445 (*exact_soln_pt)(x,exact_soln);
448 for(
unsigned i=0;
i<DIM;
i++)
450 outfile << x[
i] <<
" ";
452 outfile << exact_soln[0] << std::endl;
456 write_tecplot_zone_footer(outfile,nplot);
465 template <
unsigned DIM>
479 unsigned n_node = this->nnode();
484 unsigned n_intpt = this->integral_pt()->nweight();
487 for(
unsigned ipt=0;ipt<n_intpt;ipt++)
491 for(
unsigned i=0;
i<DIM;
i++)
493 s[
i] = this->integral_pt()->knot(ipt,
i);
497 double w = this->integral_pt()->weight(ipt);
500 double J=this->J_eulerian(s);
506 u=this->interpolated_u_poisson(s);
520 template <
unsigned DIM>
523 double& error,
double& norm)
537 unsigned n_node = nnode();
542 unsigned n_intpt = integral_pt()->nweight();
545 outfile <<
"ZONE" << std::endl;
551 for(
unsigned ipt=0;ipt<n_intpt;ipt++)
555 for(
unsigned i=0;
i<DIM;
i++)
557 s[
i] = integral_pt()->knot(ipt,
i);
561 double w = integral_pt()->weight(ipt);
564 double J=J_eulerian(s);
573 double u_fe=interpolated_u_poisson(s);
576 (*exact_soln_pt)(x,exact_soln);
579 for(
unsigned i=0;
i<DIM;
i++)
581 outfile << x[
i] <<
" ";
583 outfile << exact_soln[0] <<
" " << exact_soln[0]-u_fe << std::endl;
586 norm+=exact_soln[0]*exact_soln[0]*
W;
587 error+=(exact_soln[0]-u_fe)*(exact_soln[0]-u_fe)*
W;
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK...
virtual void fill_in_generic_residual_contribution_poisson(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Compute element residual Vector only (if flag=and/or element Jacobian matrix.
void output(std::ostream &outfile)
Output with default number of plot points.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as ...
virtual void get_dresidual_dnodal_coordinates(RankThreeTensor< double > &dresidual_dnodal_coordinates)
Compute derivatives of elemental residual vector with respect to nodal coordinates. Overwrites default implementation in FiniteElement base class. dresidual_dnodal_coordinates(l,i,j) = d res(l) / dX_{ij}.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
void compute_norm(double &norm)
Compute norm of solution: square of the L2 norm.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: x,y,u_exact or x,y,z,u_exact at n_plot^DIM plot points.
unsigned self_test()
Self-test: Return 0 for OK.