libMesh
curl_curl_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 "curl_curl_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 #include "libmesh/zero_function.h"
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("vector_fe_ex3.in");
45 
46  // Add the solution variable
47  u_var = this->add_variable ("u", FIRST, NEDELEC_ONE);
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->extra_quadrature_order = infile("extra_quadrature_order", 0);
59 
60  this->init_dirichlet_bcs();
61 
62  // Do the parent's initialization after variables and boundary constraints are defined
63  FEMSystem::init_data();
64 }
65 
67 {
68  const boundary_id_type all_ids[4] = {0, 1, 2, 3};
69  std::set<boundary_id_type> boundary_ids(all_ids, all_ids+4);
70 
71  std::vector<unsigned int> vars;
72  vars.push_back(u_var);
73 
74  ZeroFunction<Real> func;
75 
76  // this->get_dof_map().add_dirichlet_boundary(libMesh::DirichletBoundary(boundary_ids, vars, &func));
77 }
78 
80 {
81  FEMContext & c = cast_ref<FEMContext &>(context);
82 
83  // Get finite element object
86 
87  // We should prerequest all the data
88  // we will need to build the linear system.
89  fe->get_JxW();
90  fe->get_phi();
91  fe->get_curl_phi();
92  fe->get_xyz();
93 
94  // Get finite element object
96  c.get_side_fe<RealGradient>(u_var, side_fe);
97 
98  side_fe->get_JxW();
99  side_fe->get_phi();
100  side_fe->get_xyz();
101 }
102 
103 
104 bool CurlCurlSystem::element_time_derivative (bool request_jacobian,
105  DiffContext & context)
106 {
107  FEMContext & c = cast_ref<FEMContext &>(context);
108 
109  // Get finite element object
112 
113  // First we get some references to cell-specific data that
114  // will be used to assemble the linear system.
115 
116  // Element Jacobian * quadrature weights for interior integration
117  const std::vector<Real> & JxW = fe->get_JxW();
118 
119  // The velocity shape functions at interior quadrature points.
120  const std::vector<std::vector<RealGradient>> & phi = fe->get_phi();
121 
122  // The velocity shape function gradients at interior
123  // quadrature points.
124  const std::vector<std::vector<RealGradient>> & curl_phi = fe->get_curl_phi();
125 
126  const std::vector<Point> & qpoint = fe->get_xyz();
127 
128  // The number of local degrees of freedom in each variable
129  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
130 
131  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
132 
134 
135  // Now we will build the element Jacobian and residual.
136  // Constructing the residual requires the solution and its
137  // gradient from the previous timestep. This must be
138  // calculated at each quadrature point by summing the
139  // solution degree-of-freedom values by the appropriate
140  // weight functions.
141  const unsigned int n_qpoints = c.get_element_qrule().n_points();
142 
143  // Loop over quadrature points
144  for (unsigned int qp=0; qp != n_qpoints; qp++)
145  {
146  Gradient u;
147  Gradient curl_u;
148 
149  c.interior_value(u_var, qp, u);
150 
151  c.interior_curl(u_var, qp, curl_u);
152 
153  // Value of the forcing function at this quadrature point
154  RealGradient f = this->forcing(qpoint[qp]);
155 
156  // First, an i-loop over the degrees of freedom.
157  for (unsigned int i=0; i != n_u_dofs; i++)
158  {
159  Fu(i) += (curl_u*curl_phi[i][qp] + u*phi[i][qp] - f*phi[i][qp])*JxW[qp];
160 
161  // Matrix contributions for the uu and vv couplings.
162  if (request_jacobian)
163  for (unsigned int j=0; j != n_u_dofs; j++)
164  Kuu(i,j) += (curl_phi[j][qp]*curl_phi[i][qp] +
165  phi[j][qp]*phi[i][qp])*JxW[qp];
166  }
167  } // end of the quadrature point qp-loop
168 
169  return request_jacobian;
170 }
171 
172 
173 bool CurlCurlSystem::side_time_derivative (bool request_jacobian,
174  DiffContext & context)
175 {
176  FEMContext & c = cast_ref<FEMContext &>(context);
177 
178  // Get finite element object
180  c.get_side_fe<RealGradient>(u_var, side_fe);
181 
182  // First we get some references to cell-specific data that
183  // will be used to assemble the linear system.
184 
185  // Element Jacobian * quadrature weights for interior integration
186  const std::vector<Real> & JxW = side_fe->get_JxW();
187 
188  // The velocity shape functions at interior quadrature points.
189  const std::vector<std::vector<RealGradient>> & phi = side_fe->get_phi();
190 
191  // The number of local degrees of freedom in each variable
192  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
193 
194  const std::vector<Point> & normals = side_fe->get_normals();
195 
196  const std::vector<Point> & qpoint = side_fe->get_xyz();
197 
198  // The penalty value. \frac{1}{\epsilon}
199  // in the discussion above.
200  const Real penalty = 1.e10;
201 
202  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
204 
205  const unsigned int n_qpoints = c.get_side_qrule().n_points();
206 
207  for (unsigned int qp=0; qp != n_qpoints; qp++)
208  {
209  Gradient u;
210  c.side_value(u_var, qp, u);
211 
212  RealGradient N(normals[qp](0), normals[qp](1));
213 
214  Gradient u_exact = this->exact_solution(qpoint[qp]);
215 
216  Gradient Ncu = (u - u_exact).cross(N);
217 
218  for (unsigned int i=0; i != n_u_dofs; i++)
219  {
220  Fu(i) += penalty*Ncu*(phi[i][qp].cross(N))*JxW[qp];
221 
222  if (request_jacobian)
223  for (unsigned int j=0; j != n_u_dofs; j++)
224  Kuu(i,j) += penalty*(phi[j][qp].cross(N))*(phi[i][qp].cross(N))*JxW[qp];
225  }
226  }
227 
228  return request_jacobian;
229 }
230 
232 {
233  Real x = p(0); Real y = p(1);
234 
235  return soln.forcing(x, y);
236 }
237 
239 {
240  Real x = p(0); Real y = p(1);
241 
242  return soln(x, y);
243 }
This is the EquationSystems class.
void interior_curl(unsigned int var, unsigned int qp, OutputType &curl_u) const
Definition: fem_context.C:512
virtual bool side_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on side of elem to elem_residual.
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:203
CurlCurlSystem(EquationSystems &es, const std::string &name, const unsigned int number)
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
virtual void init_context(DiffContext &context)
int extra_quadrature_order
A member int that can be employed to indicate increased or reduced quadrature order.
Definition: system.h:1508
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:772
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
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
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
bool print_jacobians
Set print_jacobians to true to print J whenever it is assembled.
Definition: diff_system.h:334
This class provides a specific system class.
Definition: fem_system.h:53
Number interior_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:326
RealGradient exact_solution(const Point &p)
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
unsigned int u_var
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
virtual void time_evolving(unsigned int var)
Tells the DiffSystem that variable var is evolving with respect to time.
Definition: diff_physics.h:248
CurlCurlExactSolution soln
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Number side_value(unsigned int var, unsigned int qp) const
Definition: fem_context.C:575
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:765
virtual void init_data()
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
bool print_element_jacobians
Set print_element_jacobians to true to print each J_elem contribution.
Definition: diff_system.h:349
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:299
virtual bool element_time_derivative(bool request_jacobian, DiffContext &context)
Adds the time derivative contribution on elem to elem_residual.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
unsigned int n_points() const
Definition: quadrature.h:113
RealGradient forcing(Real x, Real y)
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.
void init_dirichlet_bcs()