35 template<
unsigned NNODE_1D>
51 Vector<double>& residuals)
54 const unsigned n_node =
nnode();
57 Shape psi(n_node), test(n_node);
58 DShape dpsidr(n_node, 1), dtestdr(n_node, 1);
75 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
82 ipt, psi, dpsidr, test, dtestdr);
85 double interpolated_r = 0;
87 double interpolated_w = 0;
88 double interpolated_laplacian_w = 0;
89 double interpolated_phi = 0;
90 double interpolated_laplacian_phi = 0;
92 double interpolated_dwdr = 0;
93 double interpolated_dlaplacian_wdr = 0;
94 double interpolated_dphidr = 0;
95 double interpolated_dlaplacian_phidr = 0;
97 double interpolated_smooth_dwdr = 0;
98 double interpolated_smooth_dphidr = 0;
99 double interpolated_continuous_d2wdr2 = 0;
100 double interpolated_continuous_d2phidr2 = 0;
106 for (
unsigned l = 0; l < n_node; l++)
122 interpolated_laplacian_w +=
nodal_value[1] * psi(l);
127 interpolated_laplacian_phi +=
nodal_value[3] * psi(l);
128 interpolated_smooth_dwdr +=
nodal_value[4] * psi(l);
129 interpolated_smooth_dphidr +=
nodal_value[5] * psi(l);
131 interpolated_continuous_d2wdr2 +=
nodal_value[4] * dpsidr(l, 0);
132 interpolated_continuous_d2phidr2 +=
nodal_value[5] * dpsidr(l, 0);
136 interpolated_dwdr +=
nodal_value[0] * dpsidr(l, 0);
137 interpolated_dlaplacian_wdr +=
nodal_value[1] * dpsidr(l, 0);
141 interpolated_dphidr +=
nodal_value[2] * dpsidr(l, 0);
142 interpolated_dlaplacian_phidr +=
nodal_value[3] * dpsidr(l, 0);
148 double W = w * interpolated_r * J;
163 for (
unsigned l = 0; l < n_node; l++)
171 residuals[local_eqn] += pressure * test(l) *
W;
174 residuals[local_eqn] +=
175 interpolated_dlaplacian_wdr * dtestdr(l, 0) *
W;
180 residuals[local_eqn] +=
182 (interpolated_continuous_d2wdr2 * interpolated_dphidr +
183 interpolated_dwdr * interpolated_continuous_d2phidr2) /
184 interpolated_r * test(l) *
W;
195 residuals[local_eqn] += interpolated_laplacian_w * test(l) *
W;
196 residuals[local_eqn] += interpolated_dwdr * dtestdr(l, 0) *
W;
205 residuals[local_eqn] += airy_forcing * test(l) *
W;
208 residuals[local_eqn] +=
209 interpolated_dlaplacian_phidr * dtestdr(l, 0) *
W;
212 residuals[local_eqn] -= interpolated_continuous_d2wdr2 *
213 interpolated_dwdr / interpolated_r * test(l) *
224 residuals[local_eqn] += interpolated_laplacian_phi * test(l) *
W;
225 residuals[local_eqn] += interpolated_dphidr * dtestdr(l, 0) *
W;
232 residuals[local_eqn] +=
233 (interpolated_dwdr - interpolated_smooth_dwdr) * test(l) *
W;
239 residuals[local_eqn] +=
240 (interpolated_dphidr - interpolated_smooth_dphidr) * test(l) *
W;
277 const Vector<double>&
s,
double& sigma_r_r,
double& sigma_phi_phi)
const
288 unsigned n_dim = this->
dim();
289 unsigned n_node = this->
nnode();
291 DShape dpsi_dr(n_node, n_dim);
300 double interpolated_r = 0;
301 double interpolated_dphi_dr = 0;
302 double interpolated_continuous_d2phi_dr2 = 0;
308 for (
unsigned l = 0; l < n_node; l++)
312 interpolated_dphi_dr +=
313 this->
nodal_value(l, smooth_dphi_dr_nodal_index) * psi(l);
314 interpolated_continuous_d2phi_dr2 +=
315 this->
nodal_value(l, smooth_dphi_dr_nodal_index) * dpsi_dr(l, 0);
319 if (interpolated_r == 0.0)
323 sigma_r_r = interpolated_continuous_d2phi_dr2;
327 sigma_r_r = interpolated_dphi_dr / interpolated_r;
329 sigma_phi_phi = interpolated_continuous_d2phi_dr2;
342 const unsigned& nplot)
352 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
358 double sigma_r_r = 0.0;
359 double sigma_phi_phi = 0.0;
364 << sigma_r_r <<
" " << sigma_phi_phi << std::endl;
375 const unsigned& nplot)
385 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
408 std::ostream& outfile,
409 const unsigned& nplot,
423 Vector<double> exact_soln(1);
427 for (
unsigned iplot = 0; iplot < num_plot_points; iplot++)
436 (*exact_soln_pt)(r, exact_soln);
439 outfile << r[0] <<
" ";
440 outfile << exact_soln[0] << std::endl;
456 std::ostream& outfile,
472 unsigned n_node =
nnode();
480 outfile <<
"ZONE" << std::endl;
484 Vector<double> exact_soln(1);
487 for (
unsigned ipt = 0; ipt < n_intpt; ipt++)
508 (*exact_soln_pt)(r, exact_soln);
511 outfile << r[0] <<
" ";
512 outfile << exact_soln[0] <<
" " << exact_soln[0] - w_fe << std::endl;
515 norm += exact_soln[0] * exact_soln[0] *
W;
516 error += (exact_soln[0] - w_fe) * (exact_soln[0] - w_fe) *
W;
524 template class AxisymFoepplvonKarmanElement<2>;
525 template class AxisymFoepplvonKarmanElement<3>;
526 template class AxisymFoepplvonKarmanElement<4>;
static const unsigned Initial_Nvalue
Static int that holds the number of variables at nodes: always the same.
static double Default_Physical_Constant_Value
Default value for physical constants.
void output(std::ostream &outfile)
Output with default number of plot points.
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Fill in the residuals with this element's contribution.
unsigned self_test()
Self-test: Return 0 for OK.
virtual double dshape_and_dtest_eulerian_at_knot_axisym_fvk(const unsigned &ipt, Shape &psi, DShape &dpsidr, Shape &test, DShape &dtestdr) const =0
Shape/test functions and derivs w.r.t. to global coords at integration point ipt; return Jacobian of ...
virtual unsigned nodal_index_fvk(const unsigned &i=0) const
Return the index at which the i-th unknown value is stored. The default value, i, is appropriate for ...
void interpolated_stress(const Vector< double > &s, double &sigma_r_r, double &sigma_phi_phi) const
Compute in-plane stresses.
bool Linear_bending_model
Flag which stores whether we are using a linear, pure bending model instead of the full non-linear Fo...
const double & eta() const
FvK parameter.
void output_fct(std::ostream &outfile, const unsigned &n_plot, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt)
Output exact soln: r,w_exact at n_plot plot points.
void compute_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_soln_pt, double &error, double &norm)
Get error against and norm of exact solution.
double interpolated_w_fvk(const Vector< double > &s) const
Return FE representation of transverse displacement.
virtual void get_pressure_fvk(const unsigned &ipt, const double &r, double &pressure) const
Get pressure term at (Eulerian) position r. This function is virtual to allow overloading in multi-ph...
virtual void get_airy_forcing_fvk(const unsigned &ipt, const double &r, double &airy_forcing) const
Get Airy forcing term at (Eulerian) position r. This function is virtual to allow overloading in mult...
virtual double J_eulerian(const Vector< double > &s) const
Return the Jacobian of mapping from local to global coordinates at local position s.
double nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n. Produces suitably interpolated values for hanging nodes...
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction")
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
int nodal_local_eqn(const unsigned &n, const unsigned &i) const
Return the local equation number corresponding to the i-th value at the n-th local node.
unsigned dim() const
Return the spatial dimension of the element, i.e. the number of local coordinates required to paramet...
unsigned nnode() const
Return the number of nodes.
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
virtual void get_s_plot(const unsigned &i, const unsigned &nplot, Vector< double > &s, const bool &shifted_to_interior=false) const
Get cector of local coordinates of plot point i (when plotting nplot points in each "coordinate direc...
virtual unsigned nplot_points(const unsigned &nplot) const
Return total number of plot points (when plotting nplot points in each "coordinate direction")
double dshape_eulerian(const Vector< double > &s, Shape &psi, DShape &dpsidx) const
Compute the geometric shape functions and also first derivatives w.r.t. global coordinates at local c...
double raw_nodal_value(const unsigned &n, const unsigned &i) const
Return the i-th value stored at local node n but do NOT take hanging nodes into account.
double raw_nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. Do not use the hanging node representation....
virtual void write_tecplot_zone_footer(std::ostream &outfile, const unsigned &nplot) const
Add tecplot zone "footer" to output stream (when plotting nplot points in each "coordinate direction"...
virtual unsigned self_test()
Self-test: Check inversion of element & do self-test for GeneralisedElement. Return 0 if OK.
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...