biharmonic_flux_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 // Config header generated by autoconfig
31 #ifdef HAVE_CONFIG_H
32  #include <oomph-lib-config.h>
33 #endif
34 
35 //oomph-lib includes
37 
38 
39 namespace oomph
40 {
41 //=============================================================================
42 /// Constructor, takes the pointer to the "bulk" element, the
43 /// index of the fixed local coordinate and its value represented
44 /// by an integer (+/- 1), indicating that the face is located
45 /// at the max. or min. value of the "fixed" local coordinate
46 /// in the bulk element. And the macro element position of the flux element
47 /// along the edge of the problem
48 //=============================================================================
49 template<>
52  const int &face_index, const unsigned& b) :
54  {
55 
56  // Let the bulk element build the FaceElement, i.e. setup the pointers
57  // to its nodes (by referring to the appropriate nodes in the bulk
58  // element), etc.
59  bulk_el_pt->build_face_element(face_index,this);
60 
61  // Initialise the prescribed-flux function pointer to zero
62  Flux0_fct_pt = 0;
63  Flux1_fct_pt = 0;
64 
65  // set the number of nodal degrees of freedom for the face element
66  Nface_nodal_dof = 2;
67 
68  //
69  Boundary = b;
70 
71  }
72 
73 
74 //=============================================================================
75 /// \short Calculate the Jacobian of the mapping between local and global
76 /// coordinates at the position s for face elements for two dimensional
77 /// problems
78 /// The jacobian of the 1D face element is computed which is dt/ds_t
79 //=============================================================================
80 template<>
82 {
83 
84  //Find the number of nodes
85  unsigned n_node = this->nnode();
86 
87  //Set up dummy memory for the shape functions
88  Shape psi(n_node,Nface_nodal_dof);
89  DShape dpsids(n_node,Nface_nodal_dof,1);
90 
91  //Get the shape functions and local derivatives
92  this->dshape_local(s,psi,dpsids);
93 
94  // computed dx_i/ds_t along edge element
95  Vector<double> interpolated_dxds_t(2);
96 
97  //Get the bulk position type corresponding to the slope
98  const unsigned slope_index = this->bulk_position_type(1);
99  for(unsigned l=0;l<n_node;l++)
100  {
101  for (unsigned d = 0; d < 2; d++)
102  {
103  interpolated_dxds_t[d] += this->node_pt(l)->x_gen(0,d) * dpsids(l,0,0)
104  + this->node_pt(l)->x_gen(slope_index,d)*dpsids(l,1,0);
105  }
106  }
107 
108  // dt/ds_t = dx_0/ds_t*t_0 + dx_1/ds_t*t_1
109  //
110  // given : t_0 = dx_0/ds_t / ( (dx_0/ds_t)^2 + (dx_1/ds_t)^2 )^0.5
111  // t_1 = dx_1/ds_t / ( (dx_0/ds_t)^2 + (dx_1/ds_t)^2 )^0.5
112  //
113  // then : dt/ds_t = ( (dx_0/ds_t)^2 + (dx_1/ds_t)^2 )^0.5
114  //
115  // (note : it is gaurantee that dt/ds_t is +ve, therefore do not need PARANOID
116  // check of jacobian)
117  double dtds_t = sqrt(interpolated_dxds_t[0]*interpolated_dxds_t[0] +
118  interpolated_dxds_t[1]*interpolated_dxds_t[1]);
119 
120  //Return the Jacobian
121  return dtds_t;
122 }
123 
124 
125 
126 //=============================================================================
127 /// \short Compute the elemental contribution to the residual vector for 2D
128 /// problem type 0 biharmonic flux elements
129 //=============================================================================
130 template<>
133 {
134  //Find out how many nodes there are in the face element
135  unsigned n_node = this->nnode();
136 
137  //set up memory for shape functions
138  Shape psi(n_node,Nface_nodal_dof);
139 
140  //setup memory for derivative of shape functions
141  DShape dpsi(n_node,Nface_nodal_dof,1);
142 
143  //Set the value of Nintpt
144  unsigned n_intpt = this->integral_pt()->nweight();
145 
146  // local coordinate position of integration point in face element
147  Vector<double> s(1);
148 
149  // edge sign adjust
150  int edge_sign = //-int(this->s_fixed_value())*int((s_fixed_index-0.5)*2);
151  - this->normal_sign();
152 
153  //Flip for the different normal conventions (TIDY THIS UP)
154  if((this->face_index() == 2) || (this->face_index() == -2))
155  {edge_sign *= -1;}
156 
157  // Get the bulk position type corresponding to the slope
158  const unsigned slope_index = this->bulk_position_type(1);
159 
160  // Ge the position type corresponding to the outer slope
161  const unsigned outer_slope_index = 3 - slope_index;
162 
163 //Loop over the integration points
164  //--------------------------------
165  for(unsigned ipt=0;ipt<n_intpt;ipt++)
166  {
167 
168  //Get the integral weight
169  double w = this->integral_pt()->weight(ipt);
170 
171  // value of local coordinate s at integration point
172  s[0] = this->integral_pt()->knot(ipt,0);
173 
174  // get the (1D) shape fn
175  this->dshape_local(s,psi,dpsi);
176 
177  // compute dx_i/ds_t and dx_i/ds_n at ipt
178  Vector<double> dxds_t(2,0.0);
179  Vector<double> dxds_n(2,0.0);
180  for (unsigned n = 0; n < n_node; n++)
181  {
182  for (unsigned d = 0; d < 2; d++)
183  {
184  for (unsigned k = 0; k < Nface_nodal_dof; k++)
185  {
186  dxds_t[d] += dpsi(n,k,0)*node_pt(n)->x_gen(slope_index*k,d);
187  dxds_n[d] += psi(n,k)*node_pt(n)->x_gen(outer_slope_index+
188  slope_index*k,d);
189  }
190  }
191  }
192 
193  // compute ds_n/dn
194  double ds_ndn = -edge_sign*sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1])/
195  (-dxds_n[0]*dxds_t[1]+dxds_t[0]*dxds_n[1]);
196 
197  // compute ds_t/dn
198  double ds_tdn = edge_sign*(dxds_t[1]*dxds_n[1]+dxds_t[0]*dxds_n[0])/
199  (sqrt(dxds_t[0]*dxds_t[0]+dxds_t[1]*dxds_t[1])*
200  (-dxds_n[0]*dxds_t[1]+dxds_t[0]*dxds_n[1]));
201 
202  // compute interpolated_m at s
203  double interpolated_m = 0.0;
204  Vector<double> m(2);
205  for (unsigned n = 0; n < n_node; n++)
206  {
207  this->node_pt(n)->get_coordinates_on_boundary(Boundary,m);
208  for (unsigned k = 0; k < Nface_nodal_dof; k++)
209  {
210  interpolated_m += psi(n,k)*m[k];
211  }
212  }
213 
214  // get the fluxes
215  double flux0 = 0.0;
216  get_flux0(interpolated_m,flux0);
217  double flux1 = 0.0;
218  get_flux1(interpolated_m,flux1);
219 
220  // get J
221  double J = J_eulerian(s);
222 
223  //Premultiply the weights and the Jacobian
224  double W = w*J;
225 
226  //Now add to the appropriate equations
227 
228 
229  //Loop over the test function nodes
230  for(unsigned n=0;n<n_node;n++)
231  {
232  // loop over the face element position types
233  for (unsigned k = 0;k < Nface_nodal_dof;k++)
234  {
235 
236  // apply contribution for flux0
237 
238  // determine bulk position type
239  unsigned bulk_p_type = slope_index*k;
240 
241  // local equation number
242  int local_eqn = this->nodal_local_eqn(n,bulk_p_type);
243 
244  // if node not pinned
245  if (local_eqn >= 0)
246  {
247  // add contribution to residual
248  residual[local_eqn] += flux1*psi(n,k)*W;
249  }
250 
251  // apply first contribution for flux1
252 
253  // if node not pinned
254  if (local_eqn >= 0)
255  {
256  // add contribution to residual
257  residual[local_eqn] -= flux0*dpsi(n,k,0)*ds_tdn*W;
258  }
259 
260  // apply second contribution for flux1
261 
262  // determine bulk position type
263  bulk_p_type = outer_slope_index + slope_index*k;
264 
265  // local equation number
266  local_eqn = this->nodal_local_eqn(n,bulk_p_type);
267 
268  // if node not pinned
269  if (local_eqn >= 0)
270  {
271  // add contribution to residual
272  residual[local_eqn] -= flux0*psi(n,k)*ds_ndn*W;
273  }
274  }
275  }
276  }
277 }
278 
279 // Ensure Flux elements are build
280 template class BiharmonicFluxElement<2>;
281 
282 }
283 
284 
285 
286 
287 
288 
289 
290 
void fill_in_generic_residual_contribution_biharmonic_flux(Vector< double > &residuals)
Add the element's contribution to its residual vector. Flux elements only make contribution to the re...
A general Finite Element class.
Definition: elements.h:1271
FluxFctPt Flux1_fct_pt
Function pointer to the prescribed flux.
unsigned Nface_nodal_dof
the number of nodal degrees of freedom for the face element basis functions
BiharmonicFluxElement()
Broken empty constructor.
static char t char * s
Definition: cfortran.h:572
FluxFctPt Flux0_fct_pt
Function pointer to the prescribed flux.
biharmonic element class
virtual void build_face_element(const int &face_index, FaceElement *face_element_pt)
Function for building a lower dimensional FaceElement on the specified face of the FiniteElement...
Definition: elements.cc:4922
double J_eulerian(const Vector< double > &s) const
Calculate the Jacobian of the mapping between local and global coordinates at the position s for face...