libMesh
naviersystem.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 #include "libmesh/getpot.h"
19 
20 #include "naviersystem.h"
21 
22 #include "libmesh/dirichlet_boundaries.h"
23 #include "libmesh/dof_map.h"
24 #include "libmesh/fe_base.h"
25 #include "libmesh/fe_interface.h"
26 #include "libmesh/fem_context.h"
27 #include "libmesh/mesh.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/string_to_enum.h"
30 #include "libmesh/zero_function.h"
31 #include "libmesh/elem.h"
32 
33 // Bring in everything from the libMesh namespace
34 using namespace libMesh;
35 
36 
37 // Boundary conditions for the 3D test case
38 class BdyFunction : public FunctionBase<Number>
39 {
40 public:
41  BdyFunction (unsigned int u_var,
42  unsigned int v_var,
43  unsigned int w_var,
44  Real Reynolds)
45  : _u_var(u_var), _v_var(v_var), _w_var(w_var), _Re(Reynolds)
46  { this->_initialized = true; }
47 
48  virtual Number operator() (const Point &, const Real = 0)
49  { libmesh_not_implemented(); }
50 
51  virtual void operator() (const Point & p,
52  const Real,
53  DenseVector<Number> & output)
54  {
55  output.zero();
56  const Real x=p(0), y=p(1), z=p(2);
57  output(_u_var) = (_Re+1)*(y*y + z*z);
58  output(_v_var) = (_Re+1)*(x*x + z*z);
59  output(_w_var) = (_Re+1)*(x*x + y*y);
60  }
61 
63  { return UniquePtr<FunctionBase<Number>> (new BdyFunction(_u_var, _v_var, _w_var, _Re)); }
64 
65 private:
66  const unsigned int _u_var, _v_var, _w_var;
67  const Real _Re;
68 };
69 
70 
72 {
73  const unsigned int dim = this->get_mesh().mesh_dimension();
74 
75  // Check the input file for Reynolds number, application type,
76  // approximation type
77  GetPot infile("navier.in");
78  Reynolds = infile("Reynolds", 1.);
79  application = infile("application", 0);
80  unsigned int pressure_p = infile("pressure_p", 1);
81  std::string fe_family = infile("fe_family", std::string("LAGRANGE"));
82 
83  // LBB needs better-than-quadratic velocities for better-than-linear
84  // pressures, and libMesh needs non-Lagrange elements for
85  // better-than-quadratic velocities.
86  libmesh_assert((pressure_p == 1) || (fe_family != "LAGRANGE"));
87 
88  FEFamily fefamily = Utility::string_to_enum<FEFamily>(fe_family);
89 
90  // Add the velocity components "u" & "v". They
91  // will be approximated using second-order approximation.
92  u_var = this->add_variable ("u", static_cast<Order>(pressure_p+1),
93  fefamily);
94  v_var = this->add_variable ("v",
95  static_cast<Order>(pressure_p+1),
96  fefamily);
97 
98  if (dim == 3)
99  w_var = this->add_variable ("w",
100  static_cast<Order>(pressure_p+1),
101  fefamily);
102  else
103  w_var = u_var;
104 
105  // Add the pressure variable "p". This will
106  // be approximated with a first-order basis,
107  // providing an LBB-stable pressure-velocity pair.
108  p_var = this->add_variable ("p",
109  static_cast<Order>(pressure_p),
110  fefamily);
111 
112  // Tell the system to march velocity forward in time, but
113  // leave p as a constraint only
114  this->time_evolving(u_var, 1);
115  this->time_evolving(v_var, 1);
116  if (dim == 3)
117  this->time_evolving(w_var, 1);
118 
119  // Useful debugging options
120  // Set verify_analytic_jacobians to 1e-6 to use
121  this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
122  this->print_jacobians = infile("print_jacobians", false);
123  this->print_element_jacobians = infile("print_element_jacobians", false);
124 
125  // Set Dirichlet boundary conditions
126  const boundary_id_type top_id = (dim==3) ? 5 : 2;
127 
128  std::set<boundary_id_type> top_bdys;
129  top_bdys.insert(top_id);
130 
131  const boundary_id_type all_ids[6] = {0, 1, 2, 3, 4, 5};
132  std::set<boundary_id_type> all_bdys(all_ids, all_ids+(dim*2));
133 
134  std::set<boundary_id_type> nontop_bdys = all_bdys;
135  nontop_bdys.erase(top_id);
136 
137  std::vector<unsigned int> u_only(1, u_var);
138  std::vector<unsigned int> vw(1, v_var), uvw(1, u_var);
139  uvw.push_back(v_var);
140  if (dim == 3)
141  {
142  vw.push_back(w_var);
143  uvw.push_back(w_var);
144  }
145 
147  ConstFunction<Number> one(1);
148  // For lid-driven cavity, set u=1,v=w=0 on the lid and u=v=w=0 elsewhere
149  if (application == 0)
150  {
151  this->get_dof_map().add_dirichlet_boundary
152  (DirichletBoundary (top_bdys, u_only, &one));
153  this->get_dof_map().add_dirichlet_boundary
154  (DirichletBoundary (top_bdys, vw, &zero));
155  this->get_dof_map().add_dirichlet_boundary
156  (DirichletBoundary (nontop_bdys, uvw, &zero));
157  }
158  // For forcing with zero wall velocity, set homogeneous Dirichlet BCs
159  else if (application == 1)
160  {
161  this->get_dof_map().add_dirichlet_boundary
162  (DirichletBoundary (all_bdys, uvw, &zero));
163  }
164  // For 3D test case with quadratic velocity field, set that field on walls
165  else if (application == 2)
166  {
167  BdyFunction bdy(u_var, v_var, w_var, Reynolds);
168  this->get_dof_map().add_dirichlet_boundary
169  (DirichletBoundary (all_bdys, uvw, &bdy));
170  }
171 
172  // Do the parent's initialization after variables and boundary constraints are defined
174 }
175 
176 
177 
179 {
180  FEMContext & c = cast_ref<FEMContext &>(context);
181 
182  FEBase * u_elem_fe;
183  FEBase * p_elem_fe;
184  FEBase * u_side_fe;
185 
186  c.get_element_fe(u_var, u_elem_fe);
187  c.get_element_fe(p_var, p_elem_fe);
188  c.get_side_fe(u_var, u_side_fe);
189 
190  // We should prerequest all the data
191  // we will need to build the linear system.
192  u_elem_fe->get_JxW();
193  u_elem_fe->get_phi();
194  u_elem_fe->get_dphi();
195  u_elem_fe->get_xyz();
196 
197  p_elem_fe->get_phi();
198 
199  u_side_fe->get_JxW();
200  u_side_fe->get_phi();
201  u_side_fe->get_xyz();
202 }
203 
204 
205 bool NavierSystem::element_time_derivative (bool request_jacobian,
206  DiffContext & context)
207 {
208  FEMContext & c = cast_ref<FEMContext &>(context);
209 
210  FEBase * u_elem_fe;
211  FEBase * p_elem_fe;
212 
213  c.get_element_fe(u_var, u_elem_fe);
214  c.get_element_fe(p_var, p_elem_fe);
215 
216  // First we get some references to cell-specific data that
217  // will be used to assemble the linear system.
218 
219  // Element Jacobian * quadrature weights for interior integration
220  const std::vector<Real> & JxW = u_elem_fe->get_JxW();
221 
222  // The velocity shape functions at interior quadrature points.
223  const std::vector<std::vector<Real>> & phi = u_elem_fe->get_phi();
224 
225  // The velocity shape function gradients at interior
226  // quadrature points.
227  const std::vector<std::vector<RealGradient>> & dphi = u_elem_fe->get_dphi();
228 
229  // The pressure shape functions at interior
230  // quadrature points.
231  const std::vector<std::vector<Real>> & psi = p_elem_fe->get_phi();
232 
233  // Physical location of the quadrature points
234  const std::vector<Point> & qpoint = u_elem_fe->get_xyz();
235 
236  // The number of local degrees of freedom in each variable
237  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
238  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
239  libmesh_assert_equal_to (n_u_dofs, c.get_dof_indices(v_var).size());
240 
241  // The subvectors and submatrices we need to fill:
242  const unsigned int dim = this->get_mesh().mesh_dimension();
243  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
244  DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
245  DenseSubMatrix<Number> & Kww = c.get_elem_jacobian(w_var, w_var);
246  DenseSubMatrix<Number> & Kuv = c.get_elem_jacobian(u_var, v_var);
247  DenseSubMatrix<Number> & Kuw = c.get_elem_jacobian(u_var, w_var);
248  DenseSubMatrix<Number> & Kvu = c.get_elem_jacobian(v_var, u_var);
249  DenseSubMatrix<Number> & Kvw = c.get_elem_jacobian(v_var, w_var);
250  DenseSubMatrix<Number> & Kwu = c.get_elem_jacobian(w_var, u_var);
251  DenseSubMatrix<Number> & Kwv = c.get_elem_jacobian(w_var, v_var);
252  DenseSubMatrix<Number> & Kup = c.get_elem_jacobian(u_var, p_var);
253  DenseSubMatrix<Number> & Kvp = c.get_elem_jacobian(v_var, p_var);
254  DenseSubMatrix<Number> & Kwp = c.get_elem_jacobian(w_var, p_var);
258 
259  // Now we will build the element Jacobian and residual.
260  // Constructing the residual requires the solution and its
261  // gradient from the previous timestep. This must be
262  // calculated at each quadrature point by summing the
263  // solution degree-of-freedom values by the appropriate
264  // weight functions.
265  unsigned int n_qpoints = c.get_element_qrule().n_points();
266 
267  // Variables to store solutions & its gradient at old Newton iterate
268  Number p, u, v, w;
269  Gradient grad_u, grad_v, grad_w;
270 
271  for (unsigned int qp=0; qp != n_qpoints; qp++)
272  {
273  // Compute the solution & its gradient at the old Newton iterate
274  c.interior_value(p_var, qp, p),
275  c.interior_value(u_var, qp, u),
276  c.interior_value(v_var, qp, v),
277  c.interior_value(w_var, qp, w);
278  c.interior_gradient(u_var, qp, grad_u),
279  c.interior_gradient(v_var, qp, grad_v),
280  c.interior_gradient(w_var, qp, grad_w);
281 
282  // Definitions for convenience. It is sometimes simpler to do a
283  // dot product if you have the full vector at your disposal.
284  NumberVectorValue U (u, v);
285  if (dim == 3)
286  U(2) = w;
287  const Number u_x = grad_u(0);
288  const Number u_y = grad_u(1);
289  const Number u_z = (dim == 3) ? grad_u(2) : 0;
290  const Number v_x = grad_v(0);
291  const Number v_y = grad_v(1);
292  const Number v_z = (dim == 3) ? grad_v(2) : 0;
293  const Number w_x = (dim == 3) ? grad_w(0) : 0;
294  const Number w_y = (dim == 3) ? grad_w(1) : 0;
295  const Number w_z = (dim == 3) ? grad_w(2) : 0;
296 
297  // Value of the forcing function at this quadrature point
298  Point f = this->forcing(qpoint[qp]);
299 
300  // First, an i-loop over the velocity degrees of freedom.
301  // We know that n_u_dofs == n_v_dofs so we can compute contributions
302  // for both at the same time.
303  for (unsigned int i=0; i != n_u_dofs; i++)
304  {
305  Fu(i) += JxW[qp] *
306  (-Reynolds*(U*grad_u)*phi[i][qp] + // convection term
307  p*dphi[i][qp](0) - // pressure term
308  (grad_u*dphi[i][qp]) + // diffusion term
309  f(0)*phi[i][qp] // forcing function
310  );
311 
312 
313  Fv(i) += JxW[qp] *
314  (-Reynolds*(U*grad_v)*phi[i][qp] + // convection term
315  p*dphi[i][qp](1) - // pressure term
316  (grad_v*dphi[i][qp]) + // diffusion term
317  f(1)*phi[i][qp] // forcing function
318  );
319 
320 
321  if (dim == 3)
322  Fw(i) += JxW[qp] *
323  (-Reynolds*(U*grad_w)*phi[i][qp] + // convection term
324  p*dphi[i][qp](2) - // pressure term
325  (grad_w*dphi[i][qp]) + // diffusion term
326  f(2)*phi[i][qp] // forcing function
327  );
328 
329 
330  // Note that the Fp block is identically zero unless we are using
331  // some kind of artificial compressibility scheme...
332 
333  if (request_jacobian && c.elem_solution_derivative)
334  {
335  libmesh_assert_equal_to (c.elem_solution_derivative, 1.0);
336 
337  // Matrix contributions for the uu and vv couplings.
338  for (unsigned int j=0; j != n_u_dofs; j++)
339  {
340  Kuu(i,j) += JxW[qp] *
341  /* convection term */ (-Reynolds*(U*dphi[j][qp])*phi[i][qp] -
342  /* diffusion term */ (dphi[i][qp]*dphi[j][qp]) -
343  /* Newton term */ Reynolds*u_x*phi[i][qp]*phi[j][qp]);
344 
345  Kuv(i,j) += JxW[qp] *
346  /* Newton term */ -Reynolds*u_y*phi[i][qp]*phi[j][qp];
347 
348  Kvv(i,j) += JxW[qp] *
349  /* convection term */ (-Reynolds*(U*dphi[j][qp])*phi[i][qp] -
350  /* diffusion term */ (dphi[i][qp]*dphi[j][qp]) -
351  /* Newton term */ Reynolds*v_y*phi[i][qp]*phi[j][qp]);
352 
353  Kvu(i,j) += JxW[qp] *
354  /* Newton term */ -Reynolds*v_x*phi[i][qp]*phi[j][qp];
355  if (dim == 3)
356  {
357  Kww(i,j) += JxW[qp] *
358  /* convection term */ (-Reynolds*(U*dphi[j][qp])*phi[i][qp] -
359  /* diffusion term */ (dphi[i][qp]*dphi[j][qp]) -
360  /* Newton term */ Reynolds*w_z*phi[i][qp]*phi[j][qp]);
361  Kuw(i,j) += JxW[qp] *
362  /* Newton term */ -Reynolds*u_z*phi[i][qp]*phi[j][qp];
363  Kvw(i,j) += JxW[qp] *
364  /* Newton term */ -Reynolds*v_z*phi[i][qp]*phi[j][qp];
365  Kwu(i,j) += JxW[qp] *
366  /* Newton term */ -Reynolds*w_x*phi[i][qp]*phi[j][qp];
367  Kwv(i,j) += JxW[qp] *
368  /* Newton term */ -Reynolds*w_y*phi[i][qp]*phi[j][qp];
369  }
370  }
371 
372  // Matrix contributions for the up and vp couplings.
373  for (unsigned int j=0; j != n_p_dofs; j++)
374  {
375  Kup(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](0);
376  Kvp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](1);
377  if (dim == 3)
378  Kwp(i,j) += JxW[qp]*psi[j][qp]*dphi[i][qp](2);
379  }
380  }
381  }
382  } // end of the quadrature point qp-loop
383 
384  return request_jacobian;
385 }
386 
387 
388 
389 bool NavierSystem::element_constraint (bool request_jacobian,
390  DiffContext & context)
391 {
392  FEMContext & c = cast_ref<FEMContext &>(context);
393 
394  FEBase * u_elem_fe;
395  FEBase * p_elem_fe;
396 
397  c.get_element_fe(u_var, u_elem_fe);
398  c.get_element_fe(p_var, p_elem_fe);
399 
400  // Here we define some references to cell-specific data that
401  // will be used to assemble the linear system.
402 
403  // Element Jacobian * quadrature weight for interior integration
404  const std::vector<Real> & JxW = u_elem_fe->get_JxW();
405 
406  // The velocity shape function gradients at interior
407  // quadrature points.
408  const std::vector<std::vector<RealGradient>> & dphi = u_elem_fe->get_dphi();
409 
410  // The pressure shape functions at interior
411  // quadrature points.
412  const std::vector<std::vector<Real>> & psi = p_elem_fe->get_phi();
413 
414  // The number of local degrees of freedom in each variable
415  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
416  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
417 
418  // The subvectors and submatrices we need to fill:
419  const unsigned int dim = this->get_mesh().mesh_dimension();
420  DenseSubMatrix<Number> & Kpu = c.get_elem_jacobian(p_var, u_var);
421  DenseSubMatrix<Number> & Kpv = c.get_elem_jacobian(p_var, v_var);
422  DenseSubMatrix<Number> & Kpw = c.get_elem_jacobian(p_var, w_var);
424 
425  // Add the constraint given by the continuity equation
426  unsigned int n_qpoints = c.get_element_qrule().n_points();
427 
428  // Variables to store solutions & its gradient at old Newton iterate
429  Gradient grad_u, grad_v, grad_w;
430 
431  for (unsigned int qp=0; qp != n_qpoints; qp++)
432  {
433  // Compute the velocity gradient at the old Newton iterate
434  c.interior_gradient(u_var, qp, grad_u),
435  c.interior_gradient(v_var, qp, grad_v),
436  c.interior_gradient(w_var, qp, grad_w);
437 
438  // Now a loop over the pressure degrees of freedom. This
439  // computes the contributions of the continuity equation.
440  for (unsigned int i=0; i != n_p_dofs; i++)
441  {
442  Fp(i) += JxW[qp] * psi[i][qp] *
443  (grad_u(0) + grad_v(1));
444  if (dim == 3)
445  Fp(i) += JxW[qp] * psi[i][qp] *
446  (grad_w(2));
447 
448  if (request_jacobian && c.get_elem_solution_derivative())
449  {
450  libmesh_assert_equal_to (c.get_elem_solution_derivative(), 1.0);
451 
452  for (unsigned int j=0; j != n_u_dofs; j++)
453  {
454  Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
455  Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
456  if (dim == 3)
457  Kpw(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](2);
458  }
459  }
460  }
461  } // end of the quadrature point qp-loop
462 
463  return request_jacobian;
464 }
465 
466 
467 
468 bool NavierSystem::side_constraint (bool request_jacobian,
469  DiffContext & context)
470 {
471  FEMContext & c = cast_ref<FEMContext &>(context);
472 
473  FEBase * p_elem_fe;
474 
475  c.get_element_fe(p_var, p_elem_fe);
476 
477  // Pin p = 0 at the origin
478  const Point zero(0., 0.);
479 
480  if (c.get_elem().contains_point(zero))
481  {
482  // The pressure penalty value. \f$ \frac{1}{\epsilon} \f$
483  const Real penalty = 1.e9;
484 
485  DenseSubMatrix<Number> & Kpp = c.get_elem_jacobian(p_var, p_var);
487  const unsigned int n_p_dofs = c.get_dof_indices(p_var).size();
488 
489  Number p;
490  c.point_value(p_var, zero, p);
491 
492  Number p_value = 0.;
493 
494  unsigned int dim = get_mesh().mesh_dimension();
495  FEType fe_type = p_elem_fe->get_fe_type();
496  Point p_master = FEInterface::inverse_map(dim, fe_type, &c.get_elem(), zero);
497 
498  std::vector<Real> point_phi(n_p_dofs);
499  for (unsigned int i=0; i != n_p_dofs; i++)
500  {
501  point_phi[i] = FEInterface::shape(dim, fe_type, &c.get_elem(), i, p_master);
502  }
503 
504  for (unsigned int i=0; i != n_p_dofs; i++)
505  {
506  Fp(i) += penalty * (p - p_value) * point_phi[i];
507  if (request_jacobian && c.get_elem_solution_derivative())
508  {
509  libmesh_assert_equal_to (c.get_elem_solution_derivative(), 1.0);
510 
511  for (unsigned int j=0; j != n_p_dofs; j++)
512  Kpp(i,j) += penalty * point_phi[i] * point_phi[j];
513  }
514  }
515  }
516 
517  return request_jacobian;
518 }
519 
520 
521 // We override the default mass_residual function,
522 // because in the non-dimensionalized Navier-Stokes equations
523 // the time derivative of velocity has a Reynolds number coefficient.
524 // Alternatively we could divide the whole equation by
525 // Reynolds number (or choose a more complicated non-dimensionalization
526 // of time), but this gives us an opportunity to demonstrate overriding
527 // FEMSystem::mass_residual()
528 bool NavierSystem::mass_residual (bool request_jacobian,
529  DiffContext & context)
530 {
531  FEMContext & c = cast_ref<FEMContext &>(context);
532 
533  FEBase * u_elem_fe;
534 
535  c.get_element_fe(u_var, u_elem_fe);
536 
537  // The subvectors and submatrices we need to fill:
538  const unsigned int dim = this->get_mesh().mesh_dimension();
539 
540  // Element Jacobian * quadrature weight for interior integration
541  const std::vector<Real> & JxW = u_elem_fe->get_JxW();
542 
543  // The velocity shape functions at interior quadrature points.
544  const std::vector<std::vector<Real>> & phi = u_elem_fe->get_phi();
545 
546  // The subvectors and submatrices we need to fill:
550  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
551  DenseSubMatrix<Number> & Kvv = c.get_elem_jacobian(v_var, v_var);
552  DenseSubMatrix<Number> & Kww = c.get_elem_jacobian(w_var, w_var);
553 
554  // The number of local degrees of freedom in velocity
555  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
556 
557  unsigned int n_qpoints = c.get_element_qrule().n_points();
558 
559  // Variables to store time derivatives at old Newton iterate
560  Number u_dot, v_dot, w_dot;
561 
562  for (unsigned int qp = 0; qp != n_qpoints; ++qp)
563  {
564  // Compute time derivatives
565  c.interior_rate(u_var, qp, u_dot),
566  c.interior_rate(v_var, qp, v_dot),
567  c.interior_rate(w_var, qp, w_dot);
568 
569  // We pull as many calculations as possible outside of loops
570  Number JxWxRe = JxW[qp] * Reynolds;
571  Number JxWxRexU = JxWxRe * u_dot;
572  Number JxWxRexV = JxWxRe * v_dot;
573  Number JxWxRexW = JxWxRe * w_dot;
574 
575  for (unsigned int i = 0; i != n_u_dofs; ++i)
576  {
577  Fu(i) -= JxWxRexU * phi[i][qp];
578  Fv(i) -= JxWxRexV * phi[i][qp];
579  if (dim == 3)
580  Fw(i) -= JxWxRexW * phi[i][qp];
581 
582  if (request_jacobian && c.get_elem_solution_derivative())
583  {
584  libmesh_assert_equal_to (c.get_elem_solution_derivative(), 1.0);
585 
586  Number JxWxRexPhiI = JxWxRe * phi[i][qp];
587  Number JxWxRexPhiII = -JxWxRexPhiI * phi[i][qp];
588  Kuu(i,i) += JxWxRexPhiII;
589  Kvv(i,i) += JxWxRexPhiII;
590  if (dim == 3)
591  Kww(i,i) += JxWxRexPhiII;
592 
593  // The mass matrix is symmetric, so we calculate
594  // one triangle and add it to both upper and lower
595  // triangles
596  for (unsigned int j = i+1; j != n_u_dofs; ++j)
597  {
598  Number Kij = -JxWxRexPhiI * phi[j][qp];
599  Kuu(i,j) += Kij;
600  Kuu(j,i) += Kij;
601  Kvv(i,j) += Kij;
602  Kvv(j,i) += Kij;
603  if (dim == 3)
604  {
605  Kww(i,j) += Kij;
606  Kww(j,i) += Kij;
607  }
608  }
609  }
610  }
611  }
612 
613  return request_jacobian;
614 }
615 
616 
617 
619 {
620  const unsigned int dim = this->get_mesh().mesh_dimension();
621 
622  Point p(1./3., 1./3.);
623  Number u = point_value(u_var, p),
624  v = point_value(v_var, p),
625  w = (dim == 3) ? point_value(w_var, p) : 0;
626 
627  libMesh::out << "u(1/3,1/3) = ("
628  << u << ", "
629  << v << ", "
630  << w << ")"
631  << std::endl;
632 }
633 
634 
635 
636 
638 {
639  switch (application)
640  {
641  // lid driven cavity
642  case 0:
643  {
644  // No forcing
645  return Point(0.,0.,0.);
646  }
647 
648  // Homogeneous Dirichlet BCs + sinusoidal forcing
649  case 1:
650  {
651  const unsigned int dim = this->get_mesh().mesh_dimension();
652 
653  // This assumes your domain is defined on [0,1]^dim.
654  Point f;
655 
656  // Counter-Clockwise vortex in x-y plane
657  if (dim==2)
658  {
659  f(0) = std::sin(2.*libMesh::pi*p(1));
660  f(1) = -std::sin(2.*libMesh::pi*p(0));
661  }
662 
663  // Counter-Clockwise vortex in x-z plane
664  else if (dim==3)
665  {
666  f(0) = std::sin(2.*libMesh::pi*p(1));
667  f(1) = 0.;
668  f(2) = -std::sin(2.*libMesh::pi*p(0));
669  }
670 
671  return f;
672  }
673 
674  // 3D test case with quadratic velocity and linear pressure field
675  case 2:
676  {
677  // This problem doesn't make sense in 1D...
678  libmesh_assert_not_equal_to (1, this->get_mesh().mesh_dimension());
679  const Real x=p(0), y=p(1), z=p(2);
680  const Real
681  u=(Reynolds+1)*(y*y + z*z),
682  v=(Reynolds+1)*(x*x + z*z),
683  w=(Reynolds+1)*(x*x + y*y);
684 
685  if (this->get_mesh().mesh_dimension() == 2)
686  return 2*Reynolds*(Reynolds+1)*Point(v*y,
687  u*x);
688  else
689  return 2*Reynolds*(Reynolds+1)*Point(v*y + w*z,
690  u*x + w*z,
691  u*x + v*y);
692  }
693 
694  default:
695  libmesh_error_msg("Invalid application id = " << application);
696  }
697 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
ConstFunction that simply returns 0.
Definition: zero_function.h:35
Number point_value(unsigned int var, const Point &p) const
Definition: fem_context.C:788
Point forcing(const Point &p)
Definition: naviersystem.C:637
unsigned int dim
virtual bool side_constraint(bool request_jacobian, DiffContext &context)
Adds the constraint contribution on side of elem to elem_residual.
Definition: naviersystem.C:468
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:262
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:178
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
const Elem & get_elem() const
Accessor for current Elem object.
Definition: fem_context.h:871
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:326
void interior_rate(unsigned int var, unsigned int qp, OutputType &u) const
Definition: fem_context.C:1273
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
Defines a dense subvector for use in finite element computations.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
FEFamily
defines an enum for finite element families.
PetscErrorCode Vec x
int8_t boundary_id_type
Definition: id_types.h:51
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: naviersystem.C:71
static Real shape(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int i, const Point &p)
Definition: fe_interface.C:641
Defines a dense submatrix for use in Finite Element-type computations.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:208
Real elem_solution_derivative
The derivative of elem_solution with respect to the current nonlinear solution.
Definition: diff_context.h:483
const unsigned int _w_var
Definition: naviersystem.C:66
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:569
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:381
virtual UniquePtr< FunctionBase< Number > > clone() const
Definition: naviersystem.C:62
BdyFunction(unsigned int u_var, unsigned int v_var, unsigned int w_var, Real Reynolds)
Definition: naviersystem.C:41
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual void init_context(DiffContext &context)
Definition: naviersystem.C:178
OStreamProxy out
Real get_elem_solution_derivative() const
The derivative of the current elem_solution w.r.t.
Definition: diff_context.h:419
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
Function that returns a single value that never changes.
virtual bool mass_residual(bool request_jacobian, DiffContext &context)
Subtracts a mass vector contribution on elem from elem_residual.
Definition: naviersystem.C:528
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: fem_system.C:845
This is the base class for functor-like classes.
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:299
const Real _Re
Definition: naviersystem.C:67
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113
sys get_dof_map()
virtual bool element_constraint(bool request_jacobian, DiffContext &context)
Adds the constraint contribution on elem to elem_residual.
Definition: naviersystem.C:389
const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:230
This class forms the foundation from which generic finite elements may be derived.
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
Definition: naviersystem.C:205
virtual bool contains_point(const Point &p, Real tol=TOLERANCE) const
Definition: elem.C:2448
virtual void postprocess()
Runs a postprocessing loop over all elements, and if postprocess_sides is true over all sides...
Definition: naviersystem.C:618
const Real pi
.
Definition: libmesh.h:172