refineable_linear_elasticity_elements.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Non-inline member functions and static member data for refineable
27// linear elasticity elements
28
29
31
32namespace oomph
33{
34 //====================================================================
35 /// Residuals for Refineable QLinearElasticityElements
36 //====================================================================
37 template<unsigned DIM>
40 Vector<double>& residuals, DenseMatrix<double>& jacobian, unsigned flag)
41 {
42#ifdef PARANOID
43
44 // Find out how many positional dofs there are
45 unsigned n_position_type = this->nnodal_position_type();
46 if (n_position_type != 1)
47 {
48 throw OomphLibError("LinearElasticity is not yet implemented for more "
49 "than one position type",
50 OOMPH_CURRENT_FUNCTION,
51 OOMPH_EXCEPTION_LOCATION);
52 }
53
54 // Throw and error if an elasticity tensor has not been set
55 if (this->Elasticity_tensor_pt == 0)
56 {
57 throw OomphLibError("No elasticity tensor set",
58 OOMPH_CURRENT_FUNCTION,
59 OOMPH_EXCEPTION_LOCATION);
60 }
61
62#endif
63
64
65 // Find out how many nodes there are
66 unsigned n_node = this->nnode();
67
68 // Find the indices at which the local displacements are stored
69 unsigned u_nodal_index[DIM];
70 for (unsigned i = 0; i < DIM; i++)
71 {
72 u_nodal_index[i] = this->u_index_linear_elasticity(i);
73 }
74
75 // Timescale ratio (non-dim density)
76 // hierher double Lambda_sq = this->lambda_sq();
77
78 // Set up memory for the shape functions
79 Shape psi(n_node);
80 DShape dpsidx(n_node, DIM);
81
82 // Set the value of Nintpt -- the number of integration points
83 unsigned n_intpt = this->integral_pt()->nweight();
84
85 // Set the vector to hold the local coordinates in the element
86 Vector<double> s(DIM);
87
88 // Integer to store the local equation number
89 int local_eqn = 0, local_unknown = 0;
90
91 // Local storage for pointers to hang_info objects
92 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
93
94 // Loop over the integration points
95 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
96 {
97 // Assign the values of s
98 for (unsigned i = 0; i < DIM; ++i)
99 {
100 s[i] = this->integral_pt()->knot(ipt, i);
101 }
102
103 // Get the integral weight
104 double w = this->integral_pt()->weight(ipt);
105
106 // Call the derivatives of the shape functions (and get Jacobian)
107 double J = this->dshape_eulerian_at_knot(ipt, psi, dpsidx);
108
109 // Storage for Eulerian coordinates (initialised to zero)
110 Vector<double> interpolated_x(DIM, 0.0);
111
112 // Calculate interpolated values of the derivative of displacements
113 DenseMatrix<double> interpolated_dudx(DIM, DIM, 0.0);
114
115 // Setup memory for accelerations (initialised to zero)
116 // hierher Vector<double> accel(DIM,0.0);
117
118
119 // Calculate displacements and derivatives
120 for (unsigned l = 0; l < n_node; l++)
121 {
122 // Loop over displacement components
123 for (unsigned i = 0; i < DIM; i++)
124 {
125 // Calculate the coordinates and the accelerations
126 interpolated_x[i] += this->nodal_position(l, i) * psi(l);
127
128 // Only compute accelerations if inertia is switched on
129 // otherwise the timestepper might not be able to
130 // work out dx_gen_dt(2,...)
131 // if (this->Unsteady)
132 // {
133 // // hierher accel[i] += this->dnodal_position_dt(2,l,i)*psi(l);
134 // }
135
136 // Get the nodal displacements
137 const double u_value = nodal_value(l, u_nodal_index[i]);
138
139 // Loop over derivative directions
140 for (unsigned j = 0; j < DIM; j++)
141 {
142 interpolated_dudx(i, j) += u_value * dpsidx(l, j);
143 }
144 }
145 }
146
147 // Get body force at current time
148 Vector<double> b(DIM);
149 this->body_force(interpolated_x, b);
150
151 // Premultiply the weights and the Jacobian
152 double W = w * J;
153
154 // Number of master nodes and storage for the weight of the shape function
155 unsigned n_master = 1;
156 double hang_weight = 1.0;
157
158 // Loop over the test functions, nodes of the element
159 for (unsigned l = 0; l < n_node; l++)
160 {
161 // Local boolean to indicate whether the node is hanging
162 bool is_node_hanging = node_pt(l)->is_hanging();
163
164 // If the node is hanging
165 if (is_node_hanging)
166 {
167 hang_info_pt = node_pt(l)->hanging_pt();
168 // Read out number of master nodes from hanging data
169 n_master = hang_info_pt->nmaster();
170 }
171 // Otherwise the node is its own master
172 else
173 {
174 n_master = 1;
175 }
176
177 // Loop over the master nodes
178 for (unsigned m = 0; m < n_master; m++)
179 {
180 // Loop over the displacement components
181 for (unsigned a = 0; a < DIM; a++)
182 {
183 // Get the equation number
184 if (is_node_hanging)
185 {
186 // Get the equation number from the master node
187 local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
188 u_nodal_index[a]);
189 // Get the hang weight from the master node
190 hang_weight = hang_info_pt->master_weight(m);
191 }
192 // Otherwise the node is not hanging
193 else
194 {
195 local_eqn = this->nodal_local_eqn(l, u_nodal_index[a]);
196 hang_weight = 1.0;
197 }
198
199 /*IF it's not a boundary condition*/
200 if (local_eqn >= 0)
201 {
202 // Acceleration and body force
203 residuals[local_eqn] += (/*hierher Lambda_sq*accel[a]*/ -b[a]) *
204 psi(l) * W * hang_weight;
205
206 // Stress term
207 for (unsigned b = 0; b < DIM; b++)
208 {
209 for (unsigned c = 0; c < DIM; c++)
210 {
211 for (unsigned d = 0; d < DIM; d++)
212 {
213 // Add the stress terms to the residuals
214 residuals[local_eqn] += this->E(a, b, c, d) *
215 interpolated_dudx(c, d) *
216 dpsidx(l, b) * W * hang_weight;
217 }
218 }
219 }
220
221 // Jacobian entries
222 if (flag)
223 {
224 // Number of master nodes and weights
225 unsigned n_master2 = 1;
226 double hang_weight2 = 1.0;
227 // Loop over the displacement basis functions again
228 for (unsigned l2 = 0; l2 < n_node; l2++)
229 {
230 // Local boolean to indicate whether the node is hanging
231 bool is_node2_hanging = node_pt(l2)->is_hanging();
232
233 // If the node is hanging
234 if (is_node2_hanging)
235 {
236 hang_info2_pt = node_pt(l2)->hanging_pt();
237 // Read out number of master nodes from hanging data
238 n_master2 = hang_info2_pt->nmaster();
239 }
240 // Otherwise the node is its own master
241 else
242 {
243 n_master2 = 1;
244 }
245
246 // Loop over the master nodes
247 for (unsigned m2 = 0; m2 < n_master2; m2++)
248 {
249 // Loop over the displacement components again
250 for (unsigned c = 0; c < DIM; c++)
251 {
252 // Get the number of the unknown
253 // If the node is hanging
254 if (is_node2_hanging)
255 {
256 // Get the equation number from the master node
257 local_unknown = this->local_hang_eqn(
258 hang_info2_pt->master_node_pt(m2), u_nodal_index[c]);
259 // Get the hang weights
260 hang_weight2 = hang_info2_pt->master_weight(m2);
261 }
262 else
263 {
264 local_unknown =
265 this->nodal_local_eqn(l2, u_nodal_index[c]);
266 hang_weight2 = 1.0;
267 }
268
269 // If it's not pinned
270 if (local_unknown >= 0)
271 {
272 for (unsigned b = 0; b < DIM; b++)
273 {
274 for (unsigned d = 0; d < DIM; d++)
275 {
276 // Add the contribution to the Jacobian matrix
277 jacobian(local_eqn, local_unknown) +=
278 this->E(a, b, c, d) * dpsidx(l2, d) *
279 dpsidx(l, b) * W * hang_weight * hang_weight2;
280 }
281 }
282 } // End of if not boundary condition
283 }
284 }
285 }
286 } // End of jacobian calculation
287
288 } // End of if not boundary condition
289 } // End of loop over coordinate directions
290 }
291 } // End of loop over shape functions
292 } // End of loop over integration points
293 }
294
295 /// Get error against and norm of exact solution
296 template<unsigned DIM>
298 std::ostream& outfile,
300 double& error,
301 double& norm)
302 {
303 // Initialise
304 error = 0.0;
305 norm = 0.0;
306
307 // Vector of local coordinates
308 Vector<double> s(DIM);
309
310 // Vector for coordintes
311 Vector<double> x(DIM);
312
313 // Set the value of n_intpt
314 unsigned n_intpt = this->integral_pt()->nweight();
315
316 // Setup output structure: Conversion is fishy but it's only output...
317 unsigned nplot;
318 if (DIM == 1)
319 {
320 nplot = n_intpt;
321 }
322 else
323 {
324 nplot = unsigned(pow(n_intpt, 1.0 / double(DIM)));
325 }
326
327 // Tecplot header info
328 outfile << this->tecplot_zone_string(nplot);
329
330 // Exact gradient Vector
331 Vector<double> exact_grad(DIM);
332
333 // Loop over the integration points
334 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
335 {
336 // Assign values of s
337 for (unsigned i = 0; i < DIM; i++)
338 {
339 s[i] = this->integral_pt()->knot(ipt, i);
340 }
341
342 // Get the integral weight
343 double w = this->integral_pt()->weight(ipt);
344
345 // Get jacobian of mapping
346 double J = this->J_eulerian(s);
347
348 // Premultiply the weights and the Jacobian
349 double W = w * J;
350
351 // Get x position as Vector
352 this->interpolated_x(s, x);
353
354 // Get FE du/dx
355 Vector<double> dudx_fe(DIM);
356
357 // Get exact gradient at this point
358 (*exact_grad_pt)(x, exact_grad);
359
360 // Output x,y,...,error
361 for (unsigned i = 0; i < DIM; i++)
362 {
363 outfile << x[i] << " ";
364 }
365 for (unsigned i = 0; i < DIM; i++)
366 {
367 outfile << exact_grad[i] << " " << exact_grad[i] - dudx_fe[i]
368 << std::endl;
369 }
370
371 // Add to error and norm
372 for (unsigned i = 0; i < DIM; i++)
373 {
374 // LinearElasticityEquations<DIM>::get_flux(s,i,dudx_fe);
375 norm += exact_grad[i] * exact_grad[i] * W;
376 error +=
377 (exact_grad[i] - dudx_fe[i]) * (exact_grad[i] - dudx_fe[i]) * W;
378 }
379 }
380 }
381
382 template<unsigned DIM>
384 {
385 if (this->tree_pt()->father_pt() != 0)
386 {
387 // Needed to set the source function pointer (if there is a father)
389 }
390 // Now do the PRefineableQElement further_build()
391 // PRefineableQElement<DIM>::further_build();
392 }
393
394
395 //====================================================================
396 /// Force building of required templates
397 //====================================================================
400
402 // template class PRefineableQLinearElasticityElement<3>;
403
404} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
void(* SteadyExactSolutionFctPt)(const Vector< double > &, Vector< double > &)
Function pointer for function that computes vector-valued steady "exact solution" as .
Definition: elements.h:1759
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
An OomphLibError object which should be thrown when an run-time error is encountered....
p-refineable version of 2D QLinearElasticityElement elements
void compute_energy_error(std::ostream &outfile, FiniteElement::SteadyExactSolutionFctPt exact_grad_pt, double &error, double &norm)
Get error against and norm of exact solution.
Class for Refineable LinearElasticity equations.
void further_build()
Further build function, pass the pointers down to the sons.
void fill_in_generic_contribution_to_residuals_linear_elasticity(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Overloaded helper function to take hanging nodes into account.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...