constrained_volume_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
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// The element-independent guts for imposition of "constant volume"
27// constraints in free surface/interface problems.
28
29
31
32namespace oomph
33{
34 //=====================================================================
35 /// Fill in the residuals for the volume constraint
36 //====================================================================
39 Vector<double>& residuals)
40 {
41 // Note: This element can only be used with the associated
42 // VolumeConstraintBoundingElement elements which compute the actual
43 // enclosed volume; here we only add the contribution to the
44 // residual; everything else, incl. the derivatives of this
45 // residual w.r.t. the nodal positions of the
46 // VolumeConstraintBoundingElements
47 // is handled by them
48 const int local_eqn = this->ptraded_local_eqn();
49 if (local_eqn >= 0)
50 {
51 residuals[local_eqn] -= *Prescribed_volume_pt;
52 }
53 }
54
55 //===========================================================================
56 /// Constructor: Pass pointer to target volume. "Pressure" value that
57 /// "traded" for the volume contraint is created internally (as a Data
58 /// item with a single pressure value)
59 //===========================================================================
61 {
62 // Store pointer to prescribed volume
63 Prescribed_volume_pt = prescribed_volume_pt;
64
65 // Create data, add as internal data and record the index
66 // (gets deleted automatically in destructor of GeneralisedElement)
69
70 // Record that we've created it as internal data
72
73 // ...and stored the "traded pressure" value as first value
75 }
76
77 //======================================================================
78 /// Constructor: Pass pointer to target volume, pointer to Data
79 /// item whose value specified by index_of_traded_pressure represents
80 /// the "Pressure" value that "traded" for the volume contraint.
81 /// The Data is stored as external Data for this element.
82 //======================================================================
84 double* prescribed_volume_pt,
87 {
88 // Store pointer to prescribed volume
89 Prescribed_volume_pt = prescribed_volume_pt;
90
91 // Add as external data and record the index
94
95 // Record that it is external data
97
98 // Record index
100 }
101
102
103 /// //////////////////////////////////////////////////////////////////////
104 /// //////////////////////////////////////////////////////////////////////
105 /// //////////////////////////////////////////////////////////////////////
106
107
108 //====================================================================
109 /// Helper function to fill in contributions to residuals
110 /// (remember that part of the residual is added by the
111 /// the associated VolumeConstraintElement). This is specific for
112 /// 1D line elements that bound 2D cartesian fluid elements.
113 //====================================================================
116 Vector<double>& residuals)
117 {
118 // Add in the volume constraint term if required
119 const int local_eqn = this->ptraded_local_eqn();
120 if (local_eqn >= 0)
121 {
122 // Find out how many nodes there are
123 const unsigned n_node = this->nnode();
124
125 // Set up memeory for the shape functions
126 Shape psif(n_node);
127 DShape dpsifds(n_node, 1);
128
129 // Set the value of n_intpt
130 const unsigned n_intpt = this->integral_pt()->nweight();
131
132 // Storage for the local coordinate
133 Vector<double> s(1);
134
135 // Loop over the integration points
136 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
137 {
138 // Get the local coordinate at the integration point
139 s[0] = this->integral_pt()->knot(ipt, 0);
140
141 // Get the integral weight
142 double W = this->integral_pt()->weight(ipt);
143
144 // Call the derivatives of the shape function at the knot point
145 this->dshape_local_at_knot(ipt, psif, dpsifds);
146
147 // Get position and tangent vector
148 Vector<double> interpolated_t1(2, 0.0);
150 for (unsigned l = 0; l < n_node; l++)
151 {
152 // Loop over directional components
153 for (unsigned i = 0; i < 2; i++)
154 {
155 interpolated_x[i] += this->nodal_position(l, i) * psif(l);
156 interpolated_t1[i] += this->nodal_position(l, i) * dpsifds(l, 0);
157 }
158 }
159
160 // Calculate the length of the tangent Vector
161 double tlength = interpolated_t1[0] * interpolated_t1[0] +
162 interpolated_t1[1] * interpolated_t1[1];
163
164 // Set the Jacobian of the line element
165 double J = sqrt(tlength);
166
167 // Now calculate the normal Vector
168 Vector<double> interpolated_n(2);
169 this->outer_unit_normal(ipt, interpolated_n);
170
171 // Assemble dot product
172 double dot = 0.0;
173 for (unsigned k = 0; k < 2; k++)
174 {
175 dot += interpolated_x[k] * interpolated_n[k];
176 }
177
178 // Add to residual with sign chosen so that the volume is
179 // positive when the elements bound the fluid
180 residuals[local_eqn] += 0.5 * dot * W * J;
181 }
182 }
183 }
184
185
186 //====================================================================
187 /// Return this element's contribution to the total volume enclosed
188 //====================================================================
190 {
191 // Initialise
192 double vol = 0.0;
193
194 // Find out how many nodes there are
195 const unsigned n_node = this->nnode();
196
197 // Set up memeory for the shape functions
198 Shape psif(n_node);
199 DShape dpsifds(n_node, 1);
200
201 // Set the value of n_intpt
202 const unsigned n_intpt = this->integral_pt()->nweight();
203
204 // Storage for the local cooridinate
205 Vector<double> s(1);
206
207 // Loop over the integration points
208 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
209 {
210 // Get the local coordinate at the integration point
211 s[0] = this->integral_pt()->knot(ipt, 0);
212
213 // Get the integral weight
214 double W = this->integral_pt()->weight(ipt);
215
216 // Call the derivatives of the shape function at the knot point
217 this->dshape_local_at_knot(ipt, psif, dpsifds);
218
219 // Get position and tangent vector
220 Vector<double> interpolated_t1(2, 0.0);
222 for (unsigned l = 0; l < n_node; l++)
223 {
224 // Loop over directional components
225 for (unsigned i = 0; i < 2; i++)
226 {
227 interpolated_x[i] += this->nodal_position(l, i) * psif(l);
228 interpolated_t1[i] += this->nodal_position(l, i) * dpsifds(l, 0);
229 }
230 }
231
232 // Calculate the length of the tangent Vector
233 double tlength = interpolated_t1[0] * interpolated_t1[0] +
234 interpolated_t1[1] * interpolated_t1[1];
235
236 // Set the Jacobian of the line element
237 double J = sqrt(tlength);
238
239 // Now calculate the normal Vector
240 Vector<double> interpolated_n(2);
241 this->outer_unit_normal(ipt, interpolated_n);
242
243 // Assemble dot product
244 double dot = 0.0;
245 for (unsigned k = 0; k < 2; k++)
246 {
247 dot += interpolated_x[k] * interpolated_n[k];
248 }
249
250 // Add to volume with sign chosen so that the volume is
251 // positive when the elements bound the fluid
252 vol += 0.5 * dot * W * J;
253 }
254
255 return vol;
256 }
257
258
259 /// ///////////////////////////////////////////////////////////////////
260 /// ///////////////////////////////////////////////////////////////////
261 /// ///////////////////////////////////////////////////////////////////
262
263 //====================================================================
264 /// Return this element's contribution to the total volume enclosed
265 //====================================================================
268 {
269 // Initialise
270 double vol = 0.0;
271
272 // Find out how many nodes there are
273 const unsigned n_node = this->nnode();
274
275 // Set up memeory for the shape functions
276 Shape psif(n_node);
277 DShape dpsifds(n_node, 1);
278
279 // Set the value of n_intpt
280 const unsigned n_intpt = this->integral_pt()->nweight();
281
282 // Storage for the local cooridinate
283 Vector<double> s(1);
284
285 // Loop over the integration points
286 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
287 {
288 // Get the local coordinate at the integration point
289 s[0] = this->integral_pt()->knot(ipt, 0);
290
291 // Get the integral weight
292 double W = this->integral_pt()->weight(ipt);
293
294 // Call the derivatives of the shape function at the knot point
295 this->dshape_local_at_knot(ipt, psif, dpsifds);
296
297 // Get position and tangent vector
298 Vector<double> interpolated_t1(2, 0.0);
300 for (unsigned l = 0; l < n_node; l++)
301 {
302 // Loop over directional components
303 for (unsigned i = 0; i < 2; i++)
304 {
305 interpolated_x[i] += this->nodal_position(l, i) * psif(l);
306 interpolated_t1[i] += this->nodal_position(l, i) * dpsifds(l, 0);
307 }
308 }
309
310 // Calculate the length of the tangent Vector
311 double tlength = interpolated_t1[0] * interpolated_t1[0] +
312 interpolated_t1[1] * interpolated_t1[1];
313
314 // Set the Jacobian of the line element
315 double J = sqrt(tlength) * interpolated_x[0];
316
317 // Now calculate the normal Vector
318 Vector<double> interpolated_n(2);
319 this->outer_unit_normal(ipt, interpolated_n);
320
321 // Assemble dot product
322 double dot = 0.0;
323 for (unsigned k = 0; k < 2; k++)
324 {
325 dot += interpolated_x[k] * interpolated_n[k];
326 }
327
328 // Add to volume with sign chosen so that the volume is
329 // positive when the elements bound the fluid
330 vol += dot * W * J / 3.0;
331 }
332
333 return vol;
334 }
335
336
337 /// ///////////////////////////////////////////////////////////////////
338 /// ///////////////////////////////////////////////////////////////////
339 /// ///////////////////////////////////////////////////////////////////
340
341
342 //====================================================================
343 /// Helper function to fill in contributions to residuals
344 /// (remember that part of the residual is added by the
345 /// the associated VolumeConstraintElement). This is specific for
346 /// axisymmetric line elements that bound 2D axisymmetric fluid elements.
347 /// The only difference from the 1D case is the addition of the
348 /// weighting factor in the integrand and division of the dot product
349 /// by three.
350 //====================================================================
353 Vector<double>& residuals)
354 {
355 // Add in the volume constraint term if required
356 const int local_eqn = this->ptraded_local_eqn();
357 if (local_eqn >= 0)
358 {
359 // Find out how many nodes there are
360 const unsigned n_node = this->nnode();
361
362 // Set up memeory for the shape functions
363 Shape psif(n_node);
364 DShape dpsifds(n_node, 1);
365
366 // Set the value of n_intpt
367 const unsigned n_intpt = this->integral_pt()->nweight();
368
369 // Storage for the local cooridinate
370 Vector<double> s(1);
371
372 // Loop over the integration points
373 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
374 {
375 // Get the local coordinate at the integration point
376 s[0] = this->integral_pt()->knot(ipt, 0);
377
378 // Get the integral weight
379 double W = this->integral_pt()->weight(ipt);
380
381 // Call the derivatives of the shape function at the knot point
382 this->dshape_local_at_knot(ipt, psif, dpsifds);
383
384 // Get position and tangent vector
385 Vector<double> interpolated_t1(2, 0.0);
387 for (unsigned l = 0; l < n_node; l++)
388 {
389 // Loop over directional components
390 for (unsigned i = 0; i < 2; i++)
391 {
392 interpolated_x[i] += this->nodal_position(l, i) * psif(l);
393 interpolated_t1[i] += this->nodal_position(l, i) * dpsifds(l, 0);
394 }
395 }
396
397 // Calculate the length of the tangent Vector
398 double tlength = interpolated_t1[0] * interpolated_t1[0] +
399 interpolated_t1[1] * interpolated_t1[1];
400
401 // Set the Jacobian of the line element
402 // multiplied by r (x[0])
403 double J = sqrt(tlength) * interpolated_x[0];
404
405 // Now calculate the normal Vector
406 Vector<double> interpolated_n(2);
407 this->outer_unit_normal(ipt, interpolated_n);
408
409 // Assemble dot product
410 double dot = 0.0;
411 for (unsigned k = 0; k < 2; k++)
412 {
413 dot += interpolated_x[k] * interpolated_n[k];
414 }
415
416 // Add to residual with sign chosen so that the volume is
417 // positive when the elements bound the fluid
418 residuals[local_eqn] += dot * W * J / 3.0;
419 }
420 }
421 }
422
423 /// ///////////////////////////////////////////////////////////////////
424 /// ///////////////////////////////////////////////////////////////////
425 /// ///////////////////////////////////////////////////////////////////
426
427
428 //=================================================================
429 /// Helper function to fill in contributions to residuals
430 /// (remember that part of the residual is added by the
431 /// the associated VolumeConstraintElement). This is specific for
432 /// 2D surface elements that bound 3D cartesian fluid elements.
433 //=================================================================
436 Vector<double>& residuals)
437 {
438 // Add in the volume constraint term if required
439 const int local_eqn = this->ptraded_local_eqn();
440 if (local_eqn >= 0)
441 {
442 // Find out how many nodes there are
443 const unsigned n_node = this->nnode();
444
445 // Set up memeory for the shape functions
446 Shape psif(n_node);
447 DShape dpsifds(n_node, 2);
448
449 // Set the value of n_intpt
450 const unsigned n_intpt = this->integral_pt()->nweight();
451
452 // Storage for the local coordinate
453 Vector<double> s(2);
454
455 // Loop over the integration points
456 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
457 {
458 // Get the local coordinate at the integration point
459 for (unsigned i = 0; i < 2; i++)
460 {
461 s[i] = this->integral_pt()->knot(ipt, i);
462 }
463
464 // Get the integral weight
465 double W = this->integral_pt()->weight(ipt);
466
467 // Call the derivatives of the shape function at the knot point
468 this->dshape_local_at_knot(ipt, psif, dpsifds);
469
470 // Get position and tangent vector
471 DenseMatrix<double> interpolated_g(2, 3, 0.0);
473 for (unsigned l = 0; l < n_node; l++)
474 {
475 // Loop over directional components
476 for (unsigned i = 0; i < 3; i++)
477 {
478 // Cache the nodal position
479 const double x_ = this->nodal_position(l, i);
480 // Calculate the position
481 interpolated_x[i] += x_ * psif(l);
482 // Calculate the two tangent vecors
483 for (unsigned j = 0; j < 2; j++)
484 {
485 interpolated_g(j, i) += x_ * dpsifds(l, j);
486 }
487 }
488 }
489
490 // Define the (non-unit) normal vector (cross product of tangent
491 // vectors)
492 Vector<double> interpolated_n(3);
493 interpolated_n[0] = interpolated_g(0, 1) * interpolated_g(1, 2) -
494 interpolated_g(0, 2) * interpolated_g(1, 1);
495 interpolated_n[1] = interpolated_g(0, 2) * interpolated_g(1, 0) -
496 interpolated_g(0, 0) * interpolated_g(1, 2);
497 interpolated_n[2] = interpolated_g(0, 0) * interpolated_g(1, 1) -
498 interpolated_g(0, 1) * interpolated_g(1, 0);
499
500 // The determinant of the local metric tensor is
501 // equal to the length of the normal vector, so
502 // rather than dividing and multipling, we just leave it out
503 // completely, which is why there is no J in the sum below
504
505 // We can now find the sign to get the OUTER UNIT normal
506 for (unsigned i = 0; i < 3; i++)
507 {
508 interpolated_n[i] *= this->normal_sign();
509 }
510
511 // Assemble dot product
512 double dot = 0.0;
513 for (unsigned k = 0; k < 3; k++)
514 {
515 dot += interpolated_x[k] * interpolated_n[k];
516 }
517
518 // Add to residual with the sign chosen such that the volume is
519 // positive when the faces are assembled such that the outer unit
520 // normal is directed out of the bulk fluid
521 residuals[local_eqn] += dot * W / 3.0;
522 }
523 }
524 }
525
526
527} // namespace oomph
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
void fill_in_generic_residual_contribution_volume_constraint(Vector< double > &residuals)
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
double contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
int & normal_sign()
Sign of outer unit normal (relative to cross-products of tangent vectors in the corresponding "bulk" ...
Definition: elements.h:4612
void outer_unit_normal(const Vector< double > &s, Vector< double > &unit_normal) const
Compute outer unit normal at the specified local coordinate.
Definition: elements.cc:6006
double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s. Overloaded to get information from bulk...
Definition: elements.h:4528
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1963
virtual void dshape_local_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsids) const
Return the geometric shape function and its derivative w.r.t. the local coordinates at the ipt-th int...
Definition: elements.cc:3239
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
double nodal_position(const unsigned &n, const unsigned &i) const
Return the i-th coordinate at local node n. If the node is hanging, the appropriate interpolation is ...
Definition: elements.h:2317
unsigned add_internal_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) internal data object to the element and return the index required to obtain it ...
Definition: elements.cc:62
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Add a (pointer to an) external data object to the element and return its index (i....
Definition: elements.cc:307
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.
void fill_in_generic_residual_contribution_volume_constraint(Vector< double > &residuals)
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
double contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
void fill_in_generic_residual_contribution_volume_constraint(Vector< double > &residuals)
Helper function to fill in contributions to residuals (remember that part of the residual is added by...
The local eqn number for the traded pressure.
void fill_in_generic_contribution_to_residuals_volume_constraint(Vector< double > &residuals)
Fill in the residuals for the volume constraint.
Index of the value in traded pressure data that corresponds to the traded pressure.
The Data that contains the traded pressure is stored as external or internal Data for the element....
VolumeConstraintElement(double *prescribed_volume_pt)
Constructor: Pass pointer to target volume. "Pressure" value that "traded" for the volume contraint i...