refineable_helmholtz_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// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
27
28
29namespace oomph
30{
31 //========================================================================
32 /// Add element's contribution to the elemental
33 /// residual vector and/or Jacobian matrix.
34 /// flag=1: compute both
35 /// flag=0: compute only residual vector
36 //========================================================================
37 template<unsigned DIM>
40 Vector<double>& residuals,
41 DenseMatrix<double>& jacobian,
42 const unsigned& flag)
43 {
44 // Find out how many nodes there are in the element
45 unsigned n_node = nnode();
46
47 // Set up memory for the shape and test functions
48 Shape psi(n_node), test(n_node);
49 DShape dpsidx(n_node, DIM), dtestdx(n_node, DIM);
50
51 // Set the value of n_intpt
52 unsigned n_intpt = integral_pt()->nweight();
53
54 // Integers to store the local equation and unknown numbers
55 int local_eqn_real = 0, local_unknown_real = 0;
56 int local_eqn_imag = 0, local_unknown_imag = 0;
57
58 // Local storage for pointers to hang_info objects
59 HangInfo *hang_info_pt = 0, *hang_info2_pt = 0;
60
61 // Loop over the integration points
62 for (unsigned ipt = 0; ipt < n_intpt; ipt++)
63 {
64 // Get the integral weight
65 double w = integral_pt()->weight(ipt);
66
67 // Call the derivatives of the shape and test functions
68 double J = this->dshape_and_dtest_eulerian_at_knot_helmholtz(
69 ipt, psi, dpsidx, test, dtestdx);
70
71 // Premultiply the weights and the Jacobian
72 double W = w * J;
73
74 // Position and gradient
75 std::complex<double> interpolated_u(0.0, 0.0);
76 Vector<double> interpolated_x(DIM, 0.0);
77 Vector<std::complex<double>> interpolated_dudx(DIM);
78
79 // Calculate function value and derivatives:
80 //-----------------------------------------
81
82 // Loop over nodes
83 for (unsigned l = 0; l < n_node; l++)
84 {
85 // Loop over directions
86 for (unsigned j = 0; j < DIM; j++)
87 {
88 interpolated_x[j] += nodal_position(l, j) * psi(l);
89 }
90 // Get the nodal value of the helmholtz unknown
91 const std::complex<double> u_value(
92 nodal_value(l, this->u_index_helmholtz().real()),
93 nodal_value(l, this->u_index_helmholtz().imag()));
94 // Add to the interpolated value
95 interpolated_u += u_value * psi(l);
96
97 // Loop over directions
98 for (unsigned j = 0; j < DIM; j++)
99 {
100 interpolated_dudx[j] += u_value * dpsidx(l, j);
101 }
102 }
103
104 // Get body force
105 std::complex<double> source(0.0, 0.0);
106 this->get_source_helmholtz(ipt, interpolated_x, source);
107
108
109 // Assemble residuals and Jacobian
110
111 // Loop over the nodes for the test functions
112 for (unsigned l = 0; l < n_node; l++)
113 {
114 // Local variables used to store the number of master nodes and the
115 // weight associated with the shape function if the node is hanging
116 unsigned n_master = 1;
117 double hang_weight = 1.0;
118
119 // Local bool (is the node hanging)
120 bool is_node_hanging = this->node_pt(l)->is_hanging();
121
122 // If the node is hanging, get the number of master nodes
123 if (is_node_hanging)
124 {
125 hang_info_pt = this->node_pt(l)->hanging_pt();
126 n_master = hang_info_pt->nmaster();
127 }
128 // Otherwise there is just one master node, the node itself
129 else
130 {
131 n_master = 1;
132 }
133
134 // Loop over the master nodes
135 for (unsigned m = 0; m < n_master; m++)
136 {
137 // Get the local equation number and hang_weight
138 // If the node is hanging
139 if (is_node_hanging)
140 {
141 // Read out the local equation number from the m-th master node
142 local_eqn_real =
143 this->local_hang_eqn(hang_info_pt->master_node_pt(m),
144 this->u_index_helmholtz().real());
145 local_eqn_imag =
146 this->local_hang_eqn(hang_info_pt->master_node_pt(m),
147 this->u_index_helmholtz().imag());
148
149 // Read out the weight from the master node
150 hang_weight = hang_info_pt->master_weight(m);
151 }
152 // If the node is not hanging
153 else
154 {
155 // The local equation number comes from the node itself
156 local_eqn_real =
157 this->nodal_local_eqn(l, this->u_index_helmholtz().real());
158 local_eqn_imag =
159 this->nodal_local_eqn(l, this->u_index_helmholtz().imag());
160
161 // The hang weight is one
162 hang_weight = 1.0;
163 }
164
165 // If the nodal equation is not a boundary condition
166 if (local_eqn_real >= 0)
167 {
168 // Add body force/source term here and Helmholtz bit
169 residuals[local_eqn_real] +=
170 (source.real() - this->k_squared() * interpolated_u.real()) *
171 test(l) * W * hang_weight;
172
173 // The Helmholtz bit itself
174 for (unsigned k = 0; k < DIM; k++)
175 {
176 residuals[local_eqn_real] +=
177 interpolated_dudx[k].real() * dtestdx(l, k) * W * hang_weight;
178 }
179
180 // Calculate the Jacobian
181 if (flag)
182 {
183 // Local variables to store the number of master nodes
184 // and the weights associated with each hanging node
185 unsigned n_master2 = 1;
186 double hang_weight2 = 1.0;
187
188 // Loop over the nodes for the variables
189 for (unsigned l2 = 0; l2 < n_node; l2++)
190 {
191 // Local bool (is the node hanging)
192 bool is_node2_hanging = this->node_pt(l2)->is_hanging();
193
194 // If the node is hanging, get the number of master nodes
195 if (is_node2_hanging)
196 {
197 hang_info2_pt = this->node_pt(l2)->hanging_pt();
198 n_master2 = hang_info2_pt->nmaster();
199 }
200 // Otherwise there is one master node, the node itself
201 else
202 {
203 n_master2 = 1;
204 }
205
206 // Loop over the master nodes
207 for (unsigned m2 = 0; m2 < n_master2; m2++)
208 {
209 // Get the local unknown and weight
210 // If the node is hanging
211 if (is_node2_hanging)
212 {
213 // Read out the local unknown from the master node
214 local_unknown_real =
215 this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
216 this->u_index_helmholtz().real());
217
218 // Read out the hanging weight from the master node
219 hang_weight2 = hang_info2_pt->master_weight(m2);
220 }
221 // If the node is not hanging
222 else
223 {
224 // The local unknown number comes from the node
225 local_unknown_real = this->nodal_local_eqn(
226 l2, this->u_index_helmholtz().real());
227
228 // The hang weight is one
229 hang_weight2 = 1.0;
230 }
231
232 // If the unknown is not pinned
233 if (local_unknown_real >= 0)
234 {
235 // Add contribution to Elemental Matrix
236 for (unsigned i = 0; i < DIM; i++)
237 {
238 jacobian(local_eqn_real, local_unknown_real) +=
239 dpsidx(l2, i) * dtestdx(l, i) * W * hang_weight *
240 hang_weight2;
241 }
242
243 // Add the helmholtz contribution
244 jacobian(local_eqn_real, local_unknown_real) -=
245 this->k_squared() * psi(l2) * test(l) * W * hang_weight *
246 hang_weight2;
247 }
248 } // End of loop over master nodes
249 } // End of loop over nodes
250 } // End of Jacobian calculation
251 } // End of case when residual equation is not pinned
252
253
254 if (local_eqn_imag >= 0)
255 {
256 // Add body force/source term here and Helmholtz bit
257 residuals[local_eqn_imag] +=
258 (source.imag() - this->k_squared() * interpolated_u.imag()) *
259 test(l) * W * hang_weight;
260
261 // The Helmholtz bit itself
262 for (unsigned k = 0; k < DIM; k++)
263 {
264 residuals[local_eqn_imag] +=
265 interpolated_dudx[k].imag() * dtestdx(l, k) * W * hang_weight;
266 }
267
268 // Calculate the Jacobian
269 if (flag)
270 {
271 // Local variables to store the number of master nodes
272 // and the weights associated with each hanging node
273 unsigned n_master2 = 1;
274 double hang_weight2 = 1.0;
275
276 // Loop over the nodes for the variables
277 for (unsigned l2 = 0; l2 < n_node; l2++)
278 {
279 // Local bool (is the node hanging)
280 bool is_node2_hanging = this->node_pt(l2)->is_hanging();
281
282 // If the node is hanging, get the number of master nodes
283 if (is_node2_hanging)
284 {
285 hang_info2_pt = this->node_pt(l2)->hanging_pt();
286 n_master2 = hang_info2_pt->nmaster();
287 }
288 // Otherwise there is one master node, the node itself
289 else
290 {
291 n_master2 = 1;
292 }
293
294 // Loop over the master nodes
295 for (unsigned m2 = 0; m2 < n_master2; m2++)
296 {
297 // Get the local unknown and weight
298 // If the node is hanging
299 if (is_node2_hanging)
300 {
301 // Read out the local unknown from the master node
302 local_unknown_imag =
303 this->local_hang_eqn(hang_info2_pt->master_node_pt(m2),
304 this->u_index_helmholtz().imag());
305
306 // Read out the hanging weight from the master node
307 hang_weight2 = hang_info2_pt->master_weight(m2);
308 }
309 // If the node is not hanging
310 else
311 {
312 // The local unknown number comes from the node
313
314 local_unknown_imag = this->nodal_local_eqn(
315 l2, this->u_index_helmholtz().imag());
316
317 // The hang weight is one
318 hang_weight2 = 1.0;
319 }
320
321 if (local_unknown_imag >= 0)
322 {
323 // Add contribution to Elemental Matrix
324 for (unsigned i = 0; i < DIM; i++)
325 {
326 jacobian(local_eqn_imag, local_unknown_imag) +=
327 dpsidx(l2, i) * dtestdx(l, i) * W * hang_weight *
328 hang_weight2;
329 }
330
331 // Add the helmholtz contribution
332 jacobian(local_eqn_imag, local_unknown_imag) -=
333 this->k_squared() * psi(l2) * test(l) * W * hang_weight *
334 hang_weight2;
335 }
336 } // End of loop over master nodes
337 } // End of loop over nodes
338 } // End of Jacobian calculation
339 }
340 } // End of loop over master nodes for residual
341 } // End of loop over nodes
342
343 } // End of loop over integration points
344 }
345
346
347 //====================================================================
348 // Force build of templates
349 //====================================================================
353
357
361
362} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
A Class for the derivatives of shape functions The class design is essentially the same as Shape,...
Definition: shape.h:278
Class that contains data for hanging nodes.
Definition: nodes.h:742
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:791
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:785
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:808
void fill_in_generic_residual_contribution_helmholtz(Vector< double > &residuals, DenseMatrix< double > &jacobian, const unsigned &flag)
Add element's contribution to elemental residual vector and/or Jacobian matrix flag=1: compute both f...
Refineable version of 2D QHelmholtzElement elements.
A Class for shape functions. In simple cases, the shape functions have only one index that can be tho...
Definition: shape.h:76
//////////////////////////////////////////////////////////////////// ////////////////////////////////...