libMesh
laplace_system.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 "laplace_system.h"
21 
22 #include "libmesh/dof_map.h"
23 #include "libmesh/fe_base.h"
24 #include "libmesh/fe_interface.h"
25 #include "libmesh/fem_context.h"
26 #include "libmesh/mesh.h"
27 #include "libmesh/quadrature.h"
28 #include "libmesh/string_to_enum.h"
29 
30 
31 // Bring in everything from the libMesh namespace
32 using namespace libMesh;
33 
35  const std::string & name_in,
36  const unsigned int number_in) :
37  FEMSystem(es, name_in, number_in)
38 {}
39 
41 {
42  // Check the input file for Reynolds number, application type,
43  // approximation type
44  GetPot infile("laplace.in");
45 
46  // Add the solution variable
47  u_var = this->add_variable ("u", FIRST, LAGRANGE_VEC);
48 
49  // The solution is evolving, with a first order time derivative
50  this->time_evolving(u_var, 1);
51 
52  // Useful debugging options
53  // Set verify_analytic_jacobians to 1e-6 to use
54  this->verify_analytic_jacobians = infile("verify_analytic_jacobians", 0.);
55  this->print_jacobians = infile("print_jacobians", false);
56  this->print_element_jacobians = infile("print_element_jacobians", false);
57 
58  this->init_dirichlet_bcs();
59 
60  // Do the parent's initialization after variables and boundary constraints are defined
61  FEMSystem::init_data();
62 }
63 
65 {
66  const boundary_id_type all_ids[6] = {0, 1, 2, 3, 4, 5};
67  std::set<boundary_id_type> boundary_ids(all_ids, all_ids+6);
68 
69  std::vector<unsigned int> vars;
70  vars.push_back(u_var);
71 
72  // Note that for vector-valued variables, it is assumed each
73  // component is stored contiguously. For 2-D elements in 3-D space,
74  // only two components should be returned.
75  SolutionFunction func(u_var);
76 
77  // In general, when reusing a system-indexed exact solution, we want
78  // to use the default system-ordering constructor for
79  // DirichletBoundary, so we demonstrate that here. In this case,
80  // though, we have only one (vector-valued!) variable, so system-
81  // and local- orderings are the same.
82  this->get_dof_map().add_dirichlet_boundary
83  (libMesh::DirichletBoundary(boundary_ids, vars, func));
84 }
85 
87 {
88  FEMContext & c = cast_ref<FEMContext &>(context);
89 
90  // Get finite element object
92  c.get_element_fe<RealGradient>(u_var, fe);
93 
94  // We should prerequest all the data
95  // we will need to build the linear system.
96  fe->get_JxW();
97  fe->get_phi();
98  fe->get_dphi();
99  fe->get_xyz();
100 
101  FEGenericBase<RealGradient> * side_fe;
102  c.get_side_fe<RealGradient>(u_var, side_fe);
103 
104  side_fe->get_JxW();
105  side_fe->get_phi();
106 }
107 
108 
109 bool LaplaceSystem::element_time_derivative (bool request_jacobian,
110  DiffContext & context)
111 {
112  FEMContext & c = cast_ref<FEMContext &>(context);
113 
114  // Get finite element object
116  c.get_element_fe<RealGradient>(u_var, fe);
117 
118  // First we get some references to cell-specific data that
119  // will be used to assemble the linear system.
120 
121  // Element Jacobian * quadrature weights for interior integration
122  const std::vector<Real> & JxW = fe->get_JxW();
123 
124  // The velocity shape functions at interior quadrature points.
125  const std::vector<std::vector<RealGradient>> & phi = fe->get_phi();
126 
127  // The velocity shape function gradients at interior
128  // quadrature points.
129  const std::vector<std::vector<RealTensor>> & grad_phi = fe->get_dphi();
130 
131  const std::vector<Point> & qpoint = fe->get_xyz();
132 
133  // The number of local degrees of freedom in each variable
134  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
135 
136  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
137 
139 
140  // Now we will build the element Jacobian and residual.
141  // Constructing the residual requires the solution and its
142  // gradient from the previous timestep. This must be
143  // calculated at each quadrature point by summing the
144  // solution degree-of-freedom values by the appropriate
145  // weight functions.
146  const unsigned int n_qpoints = c.get_element_qrule().n_points();
147 
148  for (unsigned int qp=0; qp != n_qpoints; qp++)
149  {
150  Tensor grad_u;
151 
152  c.interior_gradient(u_var, qp, grad_u);
153 
154  // Value of the forcing function at this quadrature point
155  RealGradient f = this->forcing(qpoint[qp]);
156 
157  // First, an i-loop over the velocity degrees of freedom.
158  // We know that n_u_dofs == n_v_dofs so we can compute contributions
159  // for both at the same time.
160  for (unsigned int i=0; i != n_u_dofs; i++)
161  {
162  Fu(i) += (grad_u.contract(grad_phi[i][qp]) - f*phi[i][qp])*JxW[qp];
163 
164  // Matrix contributions for the uu and vv couplings.
165  if (request_jacobian)
166  for (unsigned int j=0; j != n_u_dofs; j++)
167  Kuu(i,j) += grad_phi[j][qp].contract(grad_phi[i][qp])*JxW[qp];
168  }
169  } // end of the quadrature point qp-loop
170 
171  return request_jacobian;
172 }
173 
174 
175 
176 // bool LaplaceSystem::side_constraint (bool request_jacobian,
177 // DiffContext & context)
178 // {
179 // FEMContext & c = cast_ref<FEMContext &>(context);
180 //
181 // // Get finite element object
182 // FEGenericBase<RealGradient> * side_fe = libmesh_nullptr;
183 // c.get_side_fe<RealGradient>(u_var, side_fe);
184 //
185 // // First we get some references to cell-specific data that
186 // // will be used to assemble the linear system.
187 //
188 // // Element Jacobian * quadrature weights for interior integration
189 // const std::vector<Real> & JxW = side_fe->get_JxW();
190 //
191 // // The velocity shape functions at interior quadrature points.
192 // const std::vector<std::vector<RealGradient>> & phi = side_fe->get_phi();
193 //
194 // // The number of local degrees of freedom in each variable
195 // const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
196 //
197 // const std::vector<Point> & qpoint = side_fe->get_xyz();
198 //
199 // // The penalty value. \frac{1}{\epsilon}
200 // // in the discussion above.
201 // const Real penalty = 1.e10;
202 //
203 // DenseSubMatrix<Number> & Kuu = *c.elem_subjacobians[u_var][u_var];
204 // DenseSubVector<Number> & Fu = *c.elem_subresiduals[u_var];
205 //
206 // const unsigned int n_qpoints = (c.get_side_qrule())->n_points();
207 //
208 // for (unsigned int qp=0; qp != n_qpoints; qp++)
209 // {
210 // Gradient u;
211 // c.side_value(u_var, qp, u);
212 //
213 // Gradient u_exact(this->exact_solution(0, qpoint[qp](0), qpoint[qp](1)),
214 // this->exact_solution(1, qpoint[qp](0), qpoint[qp](1)));
215 //
216 // for (unsigned int i=0; i != n_u_dofs; i++)
217 // {
218 // Fu(i) += penalty*(u - u_exact)*phi[i][qp]*JxW[qp];
219 //
220 // if (request_jacobian)
221 // for (unsigned int j=0; j != n_u_dofs; j++)
222 // Kuu(i,j) += penalty*phi[j][qp]*phi[i][qp]*JxW[qp];
223 // }
224 // }
225 //
226 // return request_jacobian;
227 // }
228 
229 
230 
232 {
233  Real x = p(0);
234  Real y = p(1);
235  Real z = p(2);
236 
237  const Real eps = 1.e-3;
238 
239  const Real fx = -(exact_solution(0, x, y, z-eps) +
240  exact_solution(0, x, y, z+eps) +
241  exact_solution(0, x, y-eps, z) +
242  exact_solution(0, x, y+eps, z) +
243  exact_solution(0, x-eps, y, z) +
244  exact_solution(0, x+eps, y, z) -
245  6.*exact_solution(0, x, y, z))/eps/eps;
246 
247  const Real fy = -(exact_solution(1, x, y, z-eps) +
248  exact_solution(1, x, y, z+eps) +
249  exact_solution(1, x, y-eps, z) +
250  exact_solution(1, x, y+eps, z) +
251  exact_solution(1, x-eps, y, z) +
252  exact_solution(1, x+eps, y, z) -
253  6.*exact_solution(1, x, y, z))/eps/eps;
254 
255  const Real fz = -(exact_solution(2, x, y, z-eps) +
256  exact_solution(2, x, y, z+eps) +
257  exact_solution(2, x, y-eps, z) +
258  exact_solution(2, x, y+eps, z) +
259  exact_solution(2, x-eps, y, z) +
260  exact_solution(2, x+eps, y, z) -
261  6.*exact_solution(2, x, y, z))/eps/eps;
262 
263  return RealGradient(fx, fy, fz);
264 }
virtual void init_context(DiffContext &context)
Definition: L-shaped.C:34
RealVectorValue RealGradient
This is the EquationSystems class.
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
RealGradient forcing(const Point &p)
Number(* exact_solution)(const Point &p, const Parameters &, const std::string &, const std::string &)
void init_dirichlet_bcs()
LaplaceSystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
Definition: L-shaped.h:15
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:366
const class libmesh_nullptr_t libmesh_nullptr
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
Definition: L-shaped.C:58
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 DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:248
This class provides a specific system class.
Definition: fem_system.h:53
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
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
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
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:381
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together, i.e.
Definition: type_tensor.h:1175
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
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
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: L-shaped.C:17
unsigned int n_points() const
Definition: quadrature.h:113
sys get_dof_map()
This class forms the foundation from which generic finite elements may be derived.
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.