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