libMesh
laplace_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 #ifdef LIBMESH_ENABLE_DIRICHLET
67  // Note that for vector-valued variables, it is assumed each
68  // component is stored contiguously. For 2-D elements in 3-D space,
69  // only two components should be returned.
70  SolutionFunction func(u_var);
71 
72  // In general, when reusing a system-indexed exact solution, we want
73  // to use the default system-ordering constructor for
74  // DirichletBoundary, so we demonstrate that here. In this case,
75  // though, we have only one (vector-valued!) variable, so system-
76  // and local- orderings are the same.
78  (libMesh::DirichletBoundary({0,1,2,3,4,5}, {u_var}, func));
79 #endif // LIBMESH_ENABLE_DIRICHLET
80 }
81 
83 {
84  FEMContext & c = cast_ref<FEMContext &>(context);
85 
86  // Get finite element object
89 
90  // We should prerequest all the data
91  // we will need to build the linear system.
92  fe->get_JxW();
93  fe->get_phi();
94  fe->get_dphi();
95  fe->get_xyz();
96 
98  c.get_side_fe<RealGradient>(u_var, side_fe);
99 
100  side_fe->get_JxW();
101  side_fe->get_phi();
102 }
103 
104 
105 bool LaplaceSystem::element_time_derivative (bool request_jacobian,
106  DiffContext & context)
107 {
108  FEMContext & c = cast_ref<FEMContext &>(context);
109 
110  // Get finite element object
111  FEGenericBase<RealGradient> * fe = nullptr;
113 
114  // First we get some references to cell-specific data that
115  // will be used to assemble the linear system.
116 
117  // Element Jacobian * quadrature weights for interior integration
118  const std::vector<Real> & JxW = fe->get_JxW();
119 
120  // The velocity shape functions at interior quadrature points.
121  const std::vector<std::vector<RealGradient>> & phi = fe->get_phi();
122 
123  // The velocity shape function gradients at interior
124  // quadrature points.
125  const std::vector<std::vector<RealTensor>> & grad_phi = fe->get_dphi();
126 
127  const std::vector<Point> & qpoint = fe->get_xyz();
128 
129  // The number of local degrees of freedom in each variable
130  const unsigned int n_u_dofs = c.n_dof_indices(u_var);
131 
133 
135 
136  // Now we will build the element Jacobian and residual.
137  // Constructing the residual requires the solution and its
138  // gradient from the previous timestep. This must be
139  // calculated at each quadrature point by summing the
140  // solution degree-of-freedom values by the appropriate
141  // weight functions.
142  const unsigned int n_qpoints = c.get_element_qrule().n_points();
143 
144  for (unsigned int qp=0; qp != n_qpoints; qp++)
145  {
146  Tensor grad_u;
147 
148  c.interior_gradient(u_var, qp, grad_u);
149 
150  // Value of the forcing function at this quadrature point
151  RealGradient f = this->forcing(qpoint[qp]);
152 
153  // First, an i-loop over the velocity degrees of freedom.
154  // We know that n_u_dofs == n_v_dofs so we can compute contributions
155  // for both at the same time.
156  for (unsigned int i=0; i != n_u_dofs; i++)
157  {
158  Fu(i) += (grad_u.contract(grad_phi[i][qp]) - f*phi[i][qp])*JxW[qp];
159 
160  // Matrix contributions for the uu and vv couplings.
161  if (request_jacobian)
162  for (unsigned int j=0; j != n_u_dofs; j++)
163  Kuu(i,j) += grad_phi[j][qp].contract(grad_phi[i][qp])*JxW[qp];
164  }
165  } // end of the quadrature point qp-loop
166 
167  return request_jacobian;
168 }
169 
170 
171 
172 // bool LaplaceSystem::side_constraint (bool request_jacobian,
173 // DiffContext & context)
174 // {
175 // FEMContext & c = cast_ref<FEMContext &>(context);
176 //
177 // // Get finite element object
178 // FEGenericBase<RealGradient> * side_fe = nullptr;
179 // c.get_side_fe<RealGradient>(u_var, side_fe);
180 //
181 // // First we get some references to cell-specific data that
182 // // will be used to assemble the linear system.
183 //
184 // // Element Jacobian * quadrature weights for interior integration
185 // const std::vector<Real> & JxW = side_fe->get_JxW();
186 //
187 // // The velocity shape functions at interior quadrature points.
188 // const std::vector<std::vector<RealGradient>> & phi = side_fe->get_phi();
189 //
190 // // The number of local degrees of freedom in each variable
191 // const unsigned int n_u_dofs = c.dof_indices_var[u_var].size();
192 //
193 // const std::vector<Point> & qpoint = side_fe->get_xyz();
194 //
195 // // The penalty value. \frac{1}{\epsilon}
196 // // in the discussion above.
197 // const Real penalty = 1.e10;
198 //
199 // DenseSubMatrix<Number> & Kuu = *c.elem_subjacobians[u_var][u_var];
200 // DenseSubVector<Number> & Fu = *c.elem_subresiduals[u_var];
201 //
202 // const unsigned int n_qpoints = (c.get_side_qrule())->n_points();
203 //
204 // for (unsigned int qp=0; qp != n_qpoints; qp++)
205 // {
206 // Gradient u;
207 // c.side_value(u_var, qp, u);
208 //
209 // Gradient u_exact(this->exact_solution(0, qpoint[qp](0), qpoint[qp](1)),
210 // this->exact_solution(1, qpoint[qp](0), qpoint[qp](1)));
211 //
212 // for (unsigned int i=0; i != n_u_dofs; i++)
213 // {
214 // Fu(i) += penalty*(u - u_exact)*phi[i][qp]*JxW[qp];
215 //
216 // if (request_jacobian)
217 // for (unsigned int j=0; j != n_u_dofs; j++)
218 // Kuu(i,j) += penalty*phi[j][qp]*phi[i][qp]*JxW[qp];
219 // }
220 // }
221 //
222 // return request_jacobian;
223 // }
224 
225 
226 
228 {
229  Real x = p(0);
230  Real y = p(1);
231  Real z = p(2);
232 
233  const Real eps = 1.e-3;
234 
235  const Real fx = -(exact_solution(0, x, y, z-eps) +
236  exact_solution(0, x, y, z+eps) +
237  exact_solution(0, x, y-eps, z) +
238  exact_solution(0, x, y+eps, z) +
239  exact_solution(0, x-eps, y, z) +
240  exact_solution(0, x+eps, y, z) -
241  6.*exact_solution(0, x, y, z))/eps/eps;
242 
243  const Real fy = -(exact_solution(1, x, y, z-eps) +
244  exact_solution(1, x, y, z+eps) +
245  exact_solution(1, x, y-eps, z) +
246  exact_solution(1, x, y+eps, z) +
247  exact_solution(1, x-eps, y, z) +
248  exact_solution(1, x+eps, y, z) -
249  6.*exact_solution(1, x, y, z))/eps/eps;
250 
251  const Real fz = -(exact_solution(2, x, y, z-eps) +
252  exact_solution(2, x, y, z+eps) +
253  exact_solution(2, x, y-eps, z) +
254  exact_solution(2, x, y+eps, z) +
255  exact_solution(2, x-eps, y, z) +
256  exact_solution(2, x+eps, y, z) -
257  6.*exact_solution(2, x, y, z))/eps/eps;
258 
259  return RealGradient(fx, fy, fz);
260 }
virtual void init_context(DiffContext &context)
Definition: L-shaped.C:34
RealVectorValue RealGradient
This is the EquationSystems class.
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:274
Real verify_analytic_jacobians
If verify_analytic_jacobian is equal to zero (as it is by default), no numeric jacobians will be calc...
Definition: fem_system.h:215
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:55
RealGradient forcing(const Point &p)
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:317
unsigned int n_dof_indices() const
Total number of dof indices on the element.
Definition: diff_context.h:395
void init_dirichlet_bcs()
LaplaceSystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
Definition: L-shaped.h:15
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
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
unsigned int u_var
The libMesh namespace provides an interface to certain functionality in the library.
virtual void time_evolving(unsigned int var, unsigned int order)
Tells the DiffSystem that variable var is evolving with respect to time.
Definition: diff_physics.C:41
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:364
This class provides a specific system class.
Definition: fem_system.h:53
Defines a dense subvector for use in finite element computations.
Gradient interior_gradient(unsigned int var, unsigned int qp) const
Definition: fem_context.C:459
Defines a dense submatrix for use in Finite Element-type computations.
virtual_for_inffe const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:295
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
unsigned int add_variable(std::string_view var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1305
unsigned int n_points() const
Definition: quadrature.h:123
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together to return a scalar, i.e.
Definition: type_tensor.h:1268
const DenseVector< Number > & get_elem_residual() const
Const accessor for element residual.
Definition: diff_context.h:242
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool print_element_jacobians
Set print_element_jacobians to true to print each J_elem contribution.
Definition: diff_system.h:379
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
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:277
LaplaceExactSolution exact_solution
const DofMap & get_dof_map() const
Definition: system.h:2293
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
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
This class forms the foundation from which generic finite elements may be derived.
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:802
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207