refineable_unsteady_heat_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 /// Add element's contribution to the elemental
39 /// residual vector and/or Jacobian matrix.
40 /// flag=1: compute both
41 /// flag=0: compute only residual vector
42 //========================================================================
43 template<unsigned DIM>
46  DenseMatrix<double> &jacobian,
47  unsigned flag)
48 {
49 
50 //Find out how many nodes there are in the element
51 unsigned n_node = nnode();
52 
53 // Get continuous time from timestepper of first node
54 double time=node_pt(0)->time_stepper_pt()->time_pt()->time();
55 
56 //Find the index at which the unknown is stored
57  unsigned u_nodal_index = this->u_index_ust_heat();
58 
59 //Set up memory for the shape and test functions
60  Shape psi(n_node), test(n_node);
61  DShape dpsidx(n_node,DIM), dtestdx(n_node,DIM);
62 
63 //Set the value of n_intpt
64 unsigned n_intpt = integral_pt()->nweight();
65 
66 //Set the Vector to hold local coordinates
67 Vector<double> s(DIM);
68 
69  //Get Alpha and beta parameters number
70  double alpha_local = this->alpha();
71  double beta_local = this->beta();
72 
73 //Integers used to store the local equation number and local unknown
74 //indices for the residuals and jacobians
75 int local_eqn=0, local_unknown=0;
76 
77 //Local storage for pointers to hang info objects
78  HangInfo *hang_info_pt=0, *hang_info2_pt=0;
79 
80 //Local variable to determine the ALE stuff
81  bool ALE_is_disabled_flag = this->ALE_is_disabled;
82 
83 //Loop over the integration points
84 for(unsigned ipt=0;ipt<n_intpt;ipt++)
85 {
86 
87  //Assign values of s
88  for(unsigned i=0;i<DIM;i++) s[i] = integral_pt()->knot(ipt,i);
89 
90  //Get the integral weight
91  double w = integral_pt()->weight(ipt);
92 
93  //Call the derivatives of the shape and test functions
94  double J = this->
95  dshape_and_dtest_eulerian_at_knot_ust_heat(ipt,psi,dpsidx,test,dtestdx);
96 
97  //Premultiply the weights and the Jacobian
98  double W = w*J;
99 
100  //Calculate local values of the function
101  double dudt=0.0;
102  double interpolated_u=0.0;
103 
104  //This needs to be a Vector to be ANSI C++, initialise to zero
105  Vector<double> interpolated_x(DIM,0.0);
106  Vector<double> interpolated_dudx(DIM,0.0);
107  Vector<double> mesh_velocity(DIM,0.0);
108 
109 
110  //Calculate function value and derivatives:
111  //-----------------------------------------
112 
113  // Loop over nodes
114  for(unsigned l=0;l<n_node;l++)
115  {
116  //Get the value of the unknown at the node
117  double u_value = this->nodal_value(l,u_nodal_index);
118  interpolated_u += u_value*psi(l);
119  dudt += this->du_dt_ust_heat(l)*psi(l);
120  // Loop over directions
121  for(unsigned j=0;j<DIM;j++)
122  {
123  interpolated_x[j] += nodal_position(l,j)*psi(l);
124  interpolated_dudx[j] += u_value*dpsidx(l,j);
125  }
126  }
127 
128  if (!ALE_is_disabled_flag)
129  {
130  for(unsigned l=0;l<n_node;l++)
131  {
132  // Loop over directions
133  for(unsigned j=0;j<DIM;j++)
134  {
135  mesh_velocity[j] += dnodal_position_dt(l,j)*psi(l);
136  }
137  }
138  }
139 
140  //Get body force
141  double source;
142  this->get_source_ust_heat(time,ipt,interpolated_x,source);
143 
144 
145  // Assemble residuals and Jacobian
146  //================================
147 
148  // Loop over the nodes for the test functions
149  //-------------------------------------------
150  for(unsigned l=0;l<n_node;l++)
151  {
152  //Local variables to store the number of master nodes and
153  //the weight associated with the shape function if the node is hanging
154  unsigned n_master=1; double hang_weight=1.0;
155  //Local bool (is the node hanging)
156  bool is_node_hanging = this->node_pt(l)->is_hanging();
157 
158 
159  //If the node is hanging, get the number of master nodes
160  if(is_node_hanging)
161  {
162  hang_info_pt = this->node_pt(l)->hanging_pt();
163  n_master = hang_info_pt->nmaster();
164  }
165  //Otherwise there is just one master node, the node itself
166  else
167  {
168  n_master = 1;
169  }
170 
171  //Loop over the number of master nodes
172  for(unsigned m=0;m<n_master;m++)
173  {
174  //Get the local equation number and hang_weight
175  //If the node is hanging
176  if(is_node_hanging)
177  {
178  //Read out the local equation from the master node
179  local_eqn = this->local_hang_eqn(hang_info_pt->master_node_pt(m),
180  u_nodal_index);
181  //Read out the weight from the master node
182  hang_weight = hang_info_pt->master_weight(m);
183  }
184  //If the node is not hanging
185  else
186  {
187  //The local equation number comes from the node itself
188  local_eqn = this->nodal_local_eqn(l,u_nodal_index);
189  //The hang weight is one
190  hang_weight = 1.0;
191  }
192 
193  //If the nodal equation is not a boundary condition
194  if(local_eqn >= 0)
195  {
196  // Add body force/source term and time derivative
197  residuals[local_eqn]+=(alpha_local*dudt + source)*test(l)*W*hang_weight;
198 
199  // Mesh velocity and Laplace operator itself
200  for(unsigned k=0;k<DIM;k++)
201  {
202  double tmp=dtestdx(l,k)*beta_local;
203  if (!ALE_is_disabled_flag) tmp -= alpha_local*mesh_velocity[k]*test(l);
204  residuals[local_eqn] += interpolated_dudx[k]*tmp*W*hang_weight;
205  }
206 
207  // Calculate the Jacobian
208  if(flag)
209  {
210  //Local variables to store the number of master nodes
211  //and the weights associated with each hanging node
212  unsigned n_master2=1; double hang_weight2=1.0;
213  //Loop over the nodes for the variables
214  for(unsigned l2=0;l2<n_node;l2++)
215  {
216  //Local bool (is the node hanging)
217  bool is_node2_hanging = this->node_pt(l2)->is_hanging();
218  //If the node is hanging, get the number of master nodes
219  if(is_node2_hanging)
220  {
221  hang_info2_pt = this->node_pt(l2)->hanging_pt();
222  n_master2 = hang_info2_pt->nmaster();
223  }
224  //Otherwise there is one master node, the node itself
225  else
226  {
227  n_master2 = 1;
228  }
229 
230  //Loop over the master nodes
231  for(unsigned m2=0;m2<n_master2;m2++)
232  {
233  //Get the local unknown and weight
234  //If the node is hanging
235  if(is_node2_hanging)
236  {
237  //Read out the local unknown from the master node
238  local_unknown =
239  this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
240  u_nodal_index);
241  //Read out the hanging weight from the master node
242  hang_weight2 = hang_info2_pt->master_weight(m2);
243  }
244  //If the node is not hanging
245  else
246  {
247  //The local unknown number comes from the node
248  local_unknown = this->nodal_local_eqn(l2,u_nodal_index);
249  //The hang weight is one
250  hang_weight2 = 1.0;
251  }
252 
253  //If the unknown is not pinned
254  if(local_unknown >= 0)
255  {
256  //Add contribution to Elemental Matrix
257  // Mass matrix
258  jacobian(local_eqn,local_unknown)
259  += alpha_local*test(l)*psi(l2)*
260  this->node_pt(l2)->time_stepper_pt()->weight(1,0)
261  *W*hang_weight*hang_weight2;
262 
263  // Laplace operator and mesh veloc bit
264  for(unsigned k=0;k<DIM;k++)
265  {
266  double tmp=dtestdx(l,k)*beta_local;
267  if(!ALE_is_disabled_flag)
268  {
269  tmp -=
270  alpha_local*mesh_velocity[k]*test(l);
271  }
272  jacobian(local_eqn,local_unknown) +=
273  dpsidx(l2,k)*tmp*W*hang_weight*hang_weight2;
274  }
275  }
276  } //End of loop over master nodes
277  }
278  } //End of Jacobian calculation
279  }
280  } //End of loop over master nodes for residuals
281  } //End of loop over nodes
282 
283 } // End of loop over integration points
284 
285 
286 }
287 
288 //====================================================================
289 // Force build of templates
290 //====================================================================
294 
298 
299 }
bool ALE_is_disabled
Boolean flag to indicate if ALE formulation is disabled when time-derivatives are computed...
cstr elem_len * i
Definition: cfortran.h:607
void fill_in_generic_residual_contribution_ust_heat(Vector< double > &residuals, DenseMatrix< double > &jacobian, unsigned flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
static char t char * s
Definition: cfortran.h:572
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
Class that contains data for hanging nodes.
Definition: nodes.h:684
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736