refineable_young_laplace_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//====================================================================
31 
32 
33 namespace oomph
34 {
35 
36 
37 //======================================================================
38 /// Compute element residual vector taking hanging nodes into account
39 //======================================================================
42 {
43 
44  //Find out how many nodes there are
45  unsigned n_node = nnode();
46 
47  //Set up memory for the shape (test) functions
48  Shape psi(n_node);
49  DShape dpsidzeta(n_node,2);
50 
51  //Set the value of n_intpt
52  unsigned n_intpt = integral_pt()->nweight();
53 
54  //Integers to store the local equation numbers
55  int local_eqn=0;
56 
57  //Loop over the integration points
58  for(unsigned ipt=0;ipt<n_intpt;ipt++)
59  {
60  //Get the integral weight
61  double w = integral_pt()->weight(ipt);
62 
63  //Call the derivatives of the shape (test) functions
64  double J = dshape_eulerian_at_knot(ipt,psi,dpsidzeta);
65 
66  //Premultiply the weights and the Jacobian
67  double W = w*J;
68 
69  //Calculate local values of the pressure and velocity components
70  //Allocate and initialise to zero
71  double interpolated_u=0.0;
73  Vector<double> interpolated_du_dzeta(2,0.0);
74 
75  //Calculate function value and derivatives:
76  //-----------------------------------------
77  // Loop over nodes
78  for(unsigned l=0;l<n_node;l++)
79  {
80  interpolated_u += u(l)*psi(l);
81  // Loop over directions
82  for(unsigned j=0;j<2;j++)
83  {
84  interpolated_zeta[j] += nodal_position(l,j)*psi(l);
85  interpolated_du_dzeta[j] += u(l)*dpsidzeta(l,j);
86  }
87  }
88 
89  // Allocation and definition of variables necessary for
90  // further calculations
91 
92  /// "Simple" case
93  ///--------------
94  double nonlinearterm=1.0;
95  double sqnorm=0.0;
96 
97  /// Spine case
98  ///-----------
99 
100  // Derivs of position vector w.r.t. global intrinsic coords
101  Vector<Vector<double> > dRdzeta;
102  allocate_vector_of_vectors(2,3,dRdzeta);
103 
104  // Unnormalised normal
105  Vector<double> N_unnormalised(3,0.0);
106 
107  // Spine and spine basis vectors, entries initialised to zero
108  Vector<double> spine_base(3,0.0), spine(3,0.0);
109 
110  // Derivative of spine basis vector w.r.t to the intrinsic
111  // coordinates: dspine_base[i,j] = j-th component of the deriv.
112  // of the spine basis vector w.r.t. to the i-th global intrinsic
113  // coordinate
114  Vector< Vector<double> > dspine_base;
115  allocate_vector_of_vectors(2,3,dspine_base);
116 
117  // Derivative of spine vector w.r.t to the intrinsic
118  // coordinates: dspine[i,j] = j-th component of the deriv.
119  // of the spine vector w.r.t. to the i-th global intrinsic
120  // coordinate
121  Vector< Vector<double> > dspine;
122  allocate_vector_of_vectors(2,3,dspine);
123 
124  // Vector v_\alpha contains the numerator of the variations of the
125  // area element {\cal A}^{1/2} w.r.t. the components of dR/d\zeta_\alpha
126  Vector<double> area_variation_numerator_0(3,0.0);
127  Vector<double> area_variation_numerator_1(3,0.0);
128 
129  // Vector position
130  Vector<double> r(3,0.0);
131 
132  //No spines
133  //---------
134  if (!use_spines())
135  {
136 
137  for (unsigned j=0; j<2; j++)
138  sqnorm += interpolated_du_dzeta[j]*interpolated_du_dzeta[j];
139 
140  nonlinearterm=1.0/sqrt(1.0+sqnorm);
141  }
142 
143  //Spines
144  //------
145  else
146  {
147  // Get the spines
148  get_spine_base(interpolated_zeta, spine_base, dspine_base);
149  get_spine(interpolated_zeta, spine, dspine);
150 
151  // calculation of dR/d\zeta_\alpha
152  for (unsigned alpha=0;alpha<2;alpha++)
153  {
154  // Product rule for d(u {\bf S} ) / d \zeta_\alpha
155  Vector<double> dudzeta_times_spine(3,0.0);
156  scalar_times_vector(interpolated_du_dzeta[alpha],
157  spine,dudzeta_times_spine);
158 
159  Vector<double> u_times_dspinedzeta(3,0.0);
160  scalar_times_vector(interpolated_u,dspine[alpha],u_times_dspinedzeta);
161 
162  Vector<double> d_u_times_spine_dzeta(3,0.0);
163  vector_sum(dudzeta_times_spine,
164  u_times_dspinedzeta,
165  d_u_times_spine_dzeta);
166 
167  // Add derivative of spine base
168  vector_sum(d_u_times_spine_dzeta,dspine_base[alpha],dRdzeta[alpha]);
169  }
170 
171  /// Get the unnormalized normal
172  cross_product(dRdzeta[0],dRdzeta[1],N_unnormalised);
173 
174  /// Tmp storage
175  Vector<double> v_tmp_1(3,0.0);
176  Vector<double> v_tmp_2(3,0.0);
177 
178  // Calculation of
179  // |dR/d\zeta_1|^2 dR/d\zeta_0 - <dR/d\zeta_0,dR/d\zeta_1>dR/d\zeta_1
180  scalar_times_vector(pow(two_norm(dRdzeta[1]),2), dRdzeta[0], v_tmp_1);
181  scalar_times_vector(-1*scalar_product(dRdzeta[0],dRdzeta[1]),
182  dRdzeta[1], v_tmp_2);
183  vector_sum(v_tmp_1,v_tmp_2,area_variation_numerator_0);
184 
185  // Calculation of
186  // |dR/d\zeta_0|^2 dR/d\zeta_1 - <dR/d\zeta_0,dR/d\zeta_1>dR/d\zeta_0
187  scalar_times_vector(pow(two_norm(dRdzeta[0]),2), dRdzeta[1], v_tmp_1);
188  scalar_times_vector(-1*scalar_product(dRdzeta[0],dRdzeta[1]),
189  dRdzeta[0], v_tmp_2);
190  vector_sum(v_tmp_1,v_tmp_2,area_variation_numerator_1);
191 
192  // Global Eulerian cooordinates
193  for (unsigned j=0; j<3; j++)
194  {
195  r[j]=spine_base[j]+interpolated_u*spine[j];
196  }
197 
198  }
199 
200 
201 
202  // Assemble residuals
203  //-------------------
204 
205  // Loop over the nodes for the test functions
206  for(unsigned l=0;l<n_node;l++)
207  {
208  //Local variables used to store the number of master nodes and the
209  //weight associated with the shape function if the node is hanging
210  unsigned n_master=1; double hang_weight=1.0;
211 
212  //If the node is hanging, get the number of master nodes
213  if(node_pt(l)->is_hanging())
214  {
215  n_master = node_pt(l)->hanging_pt()->nmaster();
216  }
217  //Otherwise there is just one master node, the node itself
218  else
219  {
220  n_master = 1;
221  }
222 
223  //Loop over the master nodes
224  for(unsigned m=0;m<n_master;m++)
225  {
226  //Get the local equation number and hang_weight
227  //If the node is hanging
228  if(node_pt(l)->is_hanging())
229  {
230  //Read out the local equation from the master node
231  local_eqn =
232  local_hang_eqn(node_pt(l)->hanging_pt()->master_node_pt(m),0);
233  //Read out the weight from the master node
234  hang_weight = node_pt(l)->hanging_pt()->master_weight(m);
235  }
236  //If the node is not hanging
237  else
238  {
239  //The local equation number comes from the node itself
240  local_eqn = this->u_local_eqn(l);
241  //The hang weight is one
242  hang_weight = 1.0;
243  }
244 
245  /*IF it's not a boundary condition*/
246  if(local_eqn >= 0)
247  {
248 
249  // "simple" calculation case
250  if (!use_spines())
251  {
252  // Add source term: The curvature
253  residuals[local_eqn] += get_kappa()*psi(l)*W*hang_weight;
254 
255  // The YoungLaplace bit itself
256  for(unsigned k=0;k<2;k++)
257  {
258  residuals[local_eqn] += nonlinearterm*
259  interpolated_du_dzeta[k]*dpsidzeta(l,k)*W*hang_weight;
260  }
261  }
262 
263  // Spine calculation case
264  else
265  {
266 
267  // Calculation of d(u S)/d\zeta_0
268  //-------------------------------
269  Vector<double> v_tmp_1(3,0.0);
270  scalar_times_vector(dpsidzeta(l,0), spine, v_tmp_1);
271 
272  Vector<double> v_tmp_2(3,0.0);
273  scalar_times_vector(psi(l),dspine[0], v_tmp_2);
274 
275  Vector<double> d_uS_dzeta0(3,0.0);
276  vector_sum(v_tmp_1,v_tmp_2,d_uS_dzeta0);
277 
278  // Add contribution to residual
279  residuals[local_eqn] +=
280  W*hang_weight*
281  scalar_product(area_variation_numerator_0,d_uS_dzeta0)
282  /two_norm(N_unnormalised);
283 
284  // Calculation of d(u S)/d\zeta_1
285  scalar_times_vector(dpsidzeta(l,1), spine, v_tmp_1);
286  scalar_times_vector(psi(l),dspine[1], v_tmp_2);
287  Vector<double> d_uS_dzeta1(3,0.0);
288  vector_sum(v_tmp_1,v_tmp_2,d_uS_dzeta1);
289 
290  // Add contribution to residual
291  residuals[local_eqn] +=
292  W*hang_weight*
293  scalar_product(area_variation_numerator_1,d_uS_dzeta1)
294  /two_norm(N_unnormalised);
295 
296  // Curvature contribution to the residual : kappa N S test
297  residuals[local_eqn] += W*hang_weight*(get_kappa())*
298  scalar_product(N_unnormalised,spine)*psi(l);
299  }
300  }
301 
302  } //End of loop over master nodes for residual
303 
304  } //End of loop over nodes
305 
306  } // End of loop over integration points
307 
308 }
309 
310 //====================================================================
311 // Force build of templates
312 //====================================================================
313 template class RefineableQYoungLaplaceElement<2>;
314 template class RefineableQYoungLaplaceElement<3>;
315 template class RefineableQYoungLaplaceElement<4>;
316 
317 }
int local_hang_eqn(Node *const &node_pt, const unsigned &i)
Access function that returns the local equation number for the hanging node variables (values stored ...
void get_spine(const Vector< double > &x, Vector< double > &spine, Vector< Vector< double > > &dspine) const
Get spine vector field: Defaults to standard cartesian representation if no spine base fct pointers h...
virtual double u(const unsigned &n) const
virtual void get_spine_base(const Vector< double > &x, Vector< double > &spine_base, Vector< Vector< double > > &dspine_base) const
Get spine base vector field: Defaults to standard cartesian representation if no spine base fct point...
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
void fill_in_contribution_to_residuals(Vector< double > &residuals)
Compute element residual vector taking hanging nodes into account.
virtual double weight(const unsigned &i) const =0
Return weight of i-th integration point.
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
bool use_spines() const
Use spines or not? (Based on availability of function pointers to to spine and spine base vector fiel...
double interpolated_u(const Vector< double > &s) const
Return FE representation of function value u(s) at local coordinate s.
static void allocate_vector_of_vectors(unsigned n_rows, unsigned n_cols, Vector< Vector< double > > &v)
Helper fct: Allocate storage for a vector of vectors of doubles to v(n_rows,n_cols) and initialise ea...
static void scalar_times_vector(const double &lambda, const Vector< double > &v, Vector< double > &lambda_times_v)
Multiply a vector by a scalar.
Integral *const & integral_pt() const
Return the pointer to the integration scheme (const version)
Definition: elements.h:1900
static void cross_product(const Vector< double > &v1, const Vector< double > &v2, Vector< double > &v_cross)
Cross-product: v_cross= v1 x v2.
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
static double scalar_product(const Vector< double > &v1, const Vector< double > &v2)
Scalar product between two vectors.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
static void vector_sum(const Vector< double > &v1, const Vector< double > &v2, Vector< double > &vs)
Vectorial sum of two vectors.
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
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
void interpolated_zeta(const Vector< double > &s, Vector< double > &zeta) const
Calculate the interpolated value of zeta, the intrinsic coordinate of the element when viewed as a co...
Definition: elements.cc:4480
virtual int u_local_eqn(const unsigned &n)
Get the local equation number of the (one and only) unknown stored at local node n (returns -1 if val...
static double two_norm(const Vector< double > &v)
2-norm of a vector
double get_kappa() const
Get curvature.
virtual double dshape_eulerian_at_knot(const unsigned &ipt, Shape &psi, DShape &dpsidx) const
Return the geometric shape functions and also first derivatives w.r.t. global coordinates at the ipt-...
Definition: elements.cc:3252