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// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 // The element-independent guts for imposition of "constant volume"
31 // constraints in free surface/interface problems.
32 
33 
35 
36 namespace oomph
37 {
38 
39 //=====================================================================
40 /// \short Fill in the residuals for the volume constraint
41 //====================================================================
44  Vector<double> &residuals)
45  {
46  // Note: This element can only be used with the associated
47  // VolumeConstraintBoundingElement elements which compute the actual
48  // enclosed volume; here we only add the contribution to the
49  // residual; everything else, incl. the derivatives of this
50  // residual w.r.t. the nodal positions of the
51  // VolumeConstraintBoundingElements
52  // is handled by them
53  const int local_eqn = this->ptraded_local_eqn();
54  if(local_eqn >= 0)
55  {
56  residuals[local_eqn] -= *Prescribed_volume_pt;
57  }
58  }
59 
60 //===========================================================================
61 /// \short Constructor: Pass pointer to target volume. "Pressure" value that
62 /// "traded" for the volume contraint is created internally (as a Data
63 /// item with a single pressure value)
64 //===========================================================================
66  {
67  // Store pointer to prescribed volume
68  Prescribed_volume_pt = prescribed_volume_pt;
69 
70  // Create data, add as internal data and record the index
71  // (gets deleted automatically in destructor of GeneralisedElement)
73  add_internal_data(new Data(1));
74 
75  // Record that we've created it as internal data
77 
78  // ...and stored the "traded pressure" value as first value
80  }
81 
82 //======================================================================
83 /// \short Constructor: Pass pointer to target volume, pointer to Data
84 /// item whose value specified by index_of_traded_pressure represents
85 /// the "Pressure" value that "traded" for the volume contraint.
86 /// The Data is stored as external Data for this element.
87 //======================================================================
89  double* prescribed_volume_pt,
90  Data* p_traded_data_pt,
91  const unsigned& index_of_traded_pressure)
92  {
93  // Store pointer to prescribed volume
94  Prescribed_volume_pt = prescribed_volume_pt;
95 
96  // Add as external data and record the index
98  add_external_data(p_traded_data_pt);
99 
100  // Record that it is external data
102 
103  // Record index
105  }
106 
107 
108 /////////////////////////////////////////////////////////////////////////
109 /////////////////////////////////////////////////////////////////////////
110 /////////////////////////////////////////////////////////////////////////
111 
112 
113 //====================================================================
114 /// \short Helper function to fill in contributions to residuals
115 /// (remember that part of the residual is added by the
116 /// the associated VolumeConstraintElement). This is specific for
117 /// 1D line elements that bound 2D cartesian fluid elements.
118 //====================================================================
121  Vector<double> &residuals)
122  {
123  //Add in the volume constraint term if required
124  const int local_eqn=this->ptraded_local_eqn();
125  if(local_eqn >=0)
126  {
127  //Find out how many nodes there are
128  const unsigned n_node = this->nnode();
129 
130  //Set up memeory for the shape functions
131  Shape psif(n_node);
132  DShape dpsifds(n_node,1);
133 
134  //Set the value of n_intpt
135  const unsigned n_intpt = this->integral_pt()->nweight();
136 
137  //Storage for the local coordinate
138  Vector<double> s(1);
139 
140  //Loop over the integration points
141  for(unsigned ipt=0;ipt<n_intpt;ipt++)
142  {
143  //Get the local coordinate at the integration point
144  s[0] = this->integral_pt()->knot(ipt,0);
145 
146  //Get the integral weight
147  double W = this->integral_pt()->weight(ipt);
148 
149  //Call the derivatives of the shape function at the knot point
150  this->dshape_local_at_knot(ipt,psif,dpsifds);
151 
152  // Get position and tangent vector
153  Vector<double> interpolated_t1(2,0.0);
155  for(unsigned l=0;l<n_node;l++)
156  {
157  //Loop over directional components
158  for(unsigned i=0;i<2;i++)
159  {
160  interpolated_x[i] += this->nodal_position(l,i)*psif(l);
161  interpolated_t1[i] += this->nodal_position(l,i)*dpsifds(l,0);
162  }
163  }
164 
165  //Calculate the length of the tangent Vector
166  double tlength = interpolated_t1[0]*interpolated_t1[0] +
167  interpolated_t1[1]*interpolated_t1[1];
168 
169  //Set the Jacobian of the line element
170  double J = sqrt(tlength);
171 
172  //Now calculate the normal Vector
173  Vector<double> interpolated_n(2);
174  this->outer_unit_normal(ipt,interpolated_n);
175 
176  // Assemble dot product
177  double dot = 0.0;
178  for(unsigned k=0;k<2;k++)
179  {
180  dot += interpolated_x[k]*interpolated_n[k];
181  }
182 
183  // Add to residual with sign chosen so that the volume is
184  // positive when the elements bound the fluid
185  residuals[local_eqn] += 0.5*dot*W*J;
186  }
187  }
188  }
189 
190 
191 
192 //====================================================================
193 /// \short Return this element's contribution to the total volume enclosed
194 //====================================================================
196  {
197 
198  // Initialise
199  double vol=0.0;
200 
201  //Find out how many nodes there are
202  const unsigned n_node = this->nnode();
203 
204  //Set up memeory for the shape functions
205  Shape psif(n_node);
206  DShape dpsifds(n_node,1);
207 
208  //Set the value of n_intpt
209  const unsigned n_intpt = this->integral_pt()->nweight();
210 
211  //Storage for the local cooridinate
212  Vector<double> s(1);
213 
214  //Loop over the integration points
215  for(unsigned ipt=0;ipt<n_intpt;ipt++)
216  {
217  //Get the local coordinate at the integration point
218  s[0] = this->integral_pt()->knot(ipt,0);
219 
220  //Get the integral weight
221  double W = this->integral_pt()->weight(ipt);
222 
223  //Call the derivatives of the shape function at the knot point
224  this->dshape_local_at_knot(ipt,psif,dpsifds);
225 
226  // Get position and tangent vector
227  Vector<double> interpolated_t1(2,0.0);
229  for(unsigned l=0;l<n_node;l++)
230  {
231  //Loop over directional components
232  for(unsigned i=0;i<2;i++)
233  {
234  interpolated_x[i] += this->nodal_position(l,i)*psif(l);
235  interpolated_t1[i] += this->nodal_position(l,i)*dpsifds(l,0);
236  }
237  }
238 
239  //Calculate the length of the tangent Vector
240  double tlength = interpolated_t1[0]*interpolated_t1[0] +
241  interpolated_t1[1]*interpolated_t1[1];
242 
243  //Set the Jacobian of the line element
244  double J = sqrt(tlength);
245 
246  //Now calculate the normal Vector
247  Vector<double> interpolated_n(2);
248  this->outer_unit_normal(ipt,interpolated_n);
249 
250  // Assemble dot product
251  double dot = 0.0;
252  for(unsigned k=0;k<2;k++)
253  {
254  dot += interpolated_x[k]*interpolated_n[k];
255  }
256 
257  // Add to volume with sign chosen so that the volume is
258  // positive when the elements bound the fluid
259  vol += 0.5*dot*W*J;
260  }
261 
262  return vol;
263  }
264 
265 
266 //////////////////////////////////////////////////////////////////////
267 //////////////////////////////////////////////////////////////////////
268 //////////////////////////////////////////////////////////////////////
269 
270 //====================================================================
271 /// \short Return this element's contribution to the total volume enclosed
272 //====================================================================
275  {
276  // Initialise
277  double vol=0.0;
278 
279  //Find out how many nodes there are
280  const unsigned n_node = this->nnode();
281 
282  //Set up memeory for the shape functions
283  Shape psif(n_node);
284  DShape dpsifds(n_node,1);
285 
286  //Set the value of n_intpt
287  const unsigned n_intpt = this->integral_pt()->nweight();
288 
289  //Storage for the local cooridinate
290  Vector<double> s(1);
291 
292  //Loop over the integration points
293  for(unsigned ipt=0;ipt<n_intpt;ipt++)
294  {
295  //Get the local coordinate at the integration point
296  s[0] = this->integral_pt()->knot(ipt,0);
297 
298  //Get the integral weight
299  double W = this->integral_pt()->weight(ipt);
300 
301  //Call the derivatives of the shape function at the knot point
302  this->dshape_local_at_knot(ipt,psif,dpsifds);
303 
304  // Get position and tangent vector
305  Vector<double> interpolated_t1(2,0.0);
307  for(unsigned l=0;l<n_node;l++)
308  {
309  //Loop over directional components
310  for(unsigned i=0;i<2;i++)
311  {
312  interpolated_x[i] += this->nodal_position(l,i)*psif(l);
313  interpolated_t1[i] += this->nodal_position(l,i)*dpsifds(l,0);
314  }
315  }
316 
317  //Calculate the length of the tangent Vector
318  double tlength = interpolated_t1[0]*interpolated_t1[0] +
319  interpolated_t1[1]*interpolated_t1[1];
320 
321  //Set the Jacobian of the line element
322  double J = sqrt(tlength)*interpolated_x[0];
323 
324  //Now calculate the normal Vector
325  Vector<double> interpolated_n(2);
326  this->outer_unit_normal(ipt,interpolated_n);
327 
328  // Assemble dot product
329  double dot = 0.0;
330  for(unsigned k=0;k<2;k++)
331  {
332  dot += interpolated_x[k]*interpolated_n[k];
333  }
334 
335  // Add to volume with sign chosen so that the volume is
336  // positive when the elements bound the fluid
337  vol += dot*W*J/3.0;
338  }
339 
340  return vol;
341  }
342 
343 
344 //////////////////////////////////////////////////////////////////////
345 //////////////////////////////////////////////////////////////////////
346 //////////////////////////////////////////////////////////////////////
347 
348 
349 //====================================================================
350 /// \short Helper function to fill in contributions to residuals
351 /// (remember that part of the residual is added by the
352 /// the associated VolumeConstraintElement). This is specific for
353 /// axisymmetric line elements that bound 2D axisymmetric fluid elements.
354 /// The only difference from the 1D case is the addition of the
355 /// weighting factor in the integrand and division of the dot product
356 /// by three.
357 //====================================================================
360  Vector<double> &residuals)
361  {
362  //Add in the volume constraint term if required
363  const int local_eqn=this->ptraded_local_eqn();
364  if(local_eqn >=0)
365  {
366  //Find out how many nodes there are
367  const unsigned n_node = this->nnode();
368 
369  //Set up memeory for the shape functions
370  Shape psif(n_node);
371  DShape dpsifds(n_node,1);
372 
373  //Set the value of n_intpt
374  const unsigned n_intpt = this->integral_pt()->nweight();
375 
376  //Storage for the local cooridinate
377  Vector<double> s(1);
378 
379  //Loop over the integration points
380  for(unsigned ipt=0;ipt<n_intpt;ipt++)
381  {
382  //Get the local coordinate at the integration point
383  s[0] = this->integral_pt()->knot(ipt,0);
384 
385  //Get the integral weight
386  double W = this->integral_pt()->weight(ipt);
387 
388  //Call the derivatives of the shape function at the knot point
389  this->dshape_local_at_knot(ipt,psif,dpsifds);
390 
391  // Get position and tangent vector
392  Vector<double> interpolated_t1(2,0.0);
394  for(unsigned l=0;l<n_node;l++)
395  {
396  //Loop over directional components
397  for(unsigned i=0;i<2;i++)
398  {
399  interpolated_x[i] += this->nodal_position(l,i)*psif(l);
400  interpolated_t1[i] += this->nodal_position(l,i)*dpsifds(l,0);
401  }
402  }
403 
404  //Calculate the length of the tangent Vector
405  double tlength = interpolated_t1[0]*interpolated_t1[0] +
406  interpolated_t1[1]*interpolated_t1[1];
407 
408  //Set the Jacobian of the line element
409  //multiplied by r (x[0])
410  double J = sqrt(tlength)*interpolated_x[0];
411 
412  //Now calculate the normal Vector
413  Vector<double> interpolated_n(2);
414  this->outer_unit_normal(ipt,interpolated_n);
415 
416  // Assemble dot product
417  double dot = 0.0;
418  for(unsigned k=0;k<2;k++)
419  {
420  dot += interpolated_x[k]*interpolated_n[k];
421  }
422 
423  // Add to residual with sign chosen so that the volume is
424  // positive when the elements bound the fluid
425  residuals[local_eqn] += dot*W*J/3.0;
426  }
427  }
428  }
429 
430 //////////////////////////////////////////////////////////////////////
431 //////////////////////////////////////////////////////////////////////
432 //////////////////////////////////////////////////////////////////////
433 
434 
435  //=================================================================
436  /// \short Helper function to fill in contributions to residuals
437  /// (remember that part of the residual is added by the
438  /// the associated VolumeConstraintElement). This is specific for
439  /// 2D surface elements that bound 3D cartesian fluid elements.
440  //=================================================================
443  Vector<double> &residuals)
444  {
445  //Add in the volume constraint term if required
446  const int local_eqn=this->ptraded_local_eqn();
447  if(local_eqn >=0)
448  {
449  //Find out how many nodes there are
450  const unsigned n_node = this->nnode();
451 
452  //Set up memeory for the shape functions
453  Shape psif(n_node);
454  DShape dpsifds(n_node,2);
455 
456  //Set the value of n_intpt
457  const unsigned n_intpt = this->integral_pt()->nweight();
458 
459  //Storage for the local coordinate
460  Vector<double> s(2);
461 
462  //Loop over the integration points
463  for(unsigned ipt=0;ipt<n_intpt;ipt++)
464  {
465  //Get the local coordinate at the integration point
466  for(unsigned i=0;i<2;i++) {s[i] = this->integral_pt()->knot(ipt,i);}
467 
468  //Get the integral weight
469  double W = this->integral_pt()->weight(ipt);
470 
471  //Call the derivatives of the shape function at the knot point
472  this->dshape_local_at_knot(ipt,psif,dpsifds);
473 
474  // Get position and tangent vector
475  DenseMatrix<double> interpolated_g(2,3,0.0);
477  for(unsigned l=0;l<n_node;l++)
478  {
479  //Loop over directional components
480  for(unsigned i=0;i<3;i++)
481  {
482  //Cache the nodal position
483  const double x_ = this->nodal_position(l,i);
484  //Calculate the position
485  interpolated_x[i] += x_*psif(l);
486  //Calculate the two tangent vecors
487  for(unsigned j=0;j<2;j++) {interpolated_g(j,i) += x_*dpsifds(l,j);}
488  }
489  }
490 
491  // Define the (non-unit) normal vector (cross product of tangent vectors)
492  Vector<double> interpolated_n(3);
493  interpolated_n[0] =
494  interpolated_g(0,1)*interpolated_g(1,2) -
495  interpolated_g(0,2)*interpolated_g(1,1);
496  interpolated_n[1] =
497  interpolated_g(0,2)*interpolated_g(1,0) -
498  interpolated_g(0,0)*interpolated_g(1,2);
499  interpolated_n[2] =
500  interpolated_g(0,0)*interpolated_g(1,1) -
501  interpolated_g(0,1)*interpolated_g(1,0);
502 
503  //The determinant of the local metric tensor is
504  //equal to the length of the normal vector, so
505  //rather than dividing and multipling, we just leave it out
506  //completely, which is why there is no J in the sum below
507 
508  //We can now find the sign to get the OUTER UNIT normal
509  for(unsigned i=0;i<3;i++) {interpolated_n[i] *= this->normal_sign();}
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 }
528 
529 
530 
531 
532 
533 
double contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
cstr elem_len * i
Definition: cfortran.h:607
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...
void fill_in_generic_contribution_to_residuals_volume_constraint(Vector< double > &residuals)
Fill in the residuals for the volume constraint.
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:5695
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned add_external_data(Data *const &data_pt, const bool &fd=true)
Definition: elements.cc:293
int & normal_sign()
Sign of outer unit normal (relative to cross-products of tangent vectors in the corresponding "bulk" ...
Definition: elements.h:4359
bool Traded_pressure_stored_as_internal_data
The Data that contains the traded pressure is stored as external or internal Data for the element...
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...
VolumeConstraintElement(double *prescribed_volume_pt)
Constructor: Pass pointer to target volume. "Pressure" value that "traded" for the volume contraint i...
static char t char * s
Definition: cfortran.h:572
virtual double knot(const unsigned &i, const unsigned &j) const =0
Return local coordinate s[j] of i-th integration point.
double contribution_to_enclosed_volume()
Return this element's contribution to the total volume enclosed.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
int ptraded_local_eqn()
The local eqn number for the traded pressure.
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
double * Prescribed_volume_pt
Pointer to the desired value of the volume.
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:3174
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:2215
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 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:4277
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:289
virtual unsigned nweight() const =0
Return the number of integration points of the scheme.
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
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:66
int ptraded_local_eqn()
The local eqn number for 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...
unsigned index_of_traded_pressure()
Return the index of Data object at which the traded pressure is stored.