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
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// 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)
68 add_internal_data(new Data(1));
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,
85 Data* p_traded_data_pt,
86 const unsigned& index_of_traded_pressure)
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...
int ptraded_local_eqn()
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.
unsigned Index_of_traded_pressure_value
Index of the value in traded pressure data that corresponds to the traded pressure.
unsigned External_or_internal_data_index_of_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...
bool Traded_pressure_stored_as_internal_data
The Data that contains the traded pressure is stored as external or internal Data for the element....
int ptraded_local_eqn()
The local eqn number for the traded pressure.
unsigned index_of_traded_pressure()
Return the index of Data object at which the traded pressure is stored.
double * Prescribed_volume_pt
Pointer to the desired value of the volume.
Data * p_traded_data_pt()
Access to Data that contains the traded pressure.
double dot(const Vector< double > &a, const Vector< double > &b)
Probably not always best/fastest because not optimised for dimension but useful...
Definition: Vector.h:291
//////////////////////////////////////////////////////////////////// ////////////////////////////////...