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_ex4.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[6] = {0, 1, 2, 3, 4, 5};
69  std::set<boundary_id_type> boundary_ids(all_ids, all_ids+6);
70 
71  std::vector<unsigned int> vars;
72  vars.push_back(u_var);
73 
74  ZeroFunction<Real> func;
75 
76  // HCurl DirichletBoundary not supported yet
77  // this->get_dof_map().add_dirichlet_boundary(libMesh::DirichletBoundary(boundary_ids, vars, &func));
78 }
79 
81 {
82  FEMContext & c = cast_ref<FEMContext &>(context);
83 
84  // Get finite element object
87 
88  // We should prerequest all the data
89  // we will need to build the linear system.
90  fe->get_JxW();
91  fe->get_phi();
92  fe->get_curl_phi();
93  fe->get_xyz();
94 
95  // Get finite element object
97  c.get_side_fe<RealGradient>(u_var, side_fe);
98 
99  side_fe->get_JxW();
100  side_fe->get_phi();
101  side_fe->get_xyz();
102 }
103 
104 
105 bool CurlCurlSystem::element_time_derivative (bool request_jacobian,
106  DiffContext & context)
107 {
108  FEMContext & c = cast_ref<FEMContext &>(context);
109 
110  // Get finite element object
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<RealGradient>> & curl_phi = fe->get_curl_phi();
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.get_dof_indices(u_var).size();
131 
132  DenseSubMatrix<Number> & Kuu = c.get_elem_jacobian(u_var, u_var);
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  // Loop over quadrature points
145  for (unsigned int qp=0; qp != n_qpoints; qp++)
146  {
147  Gradient u;
148  Gradient curl_u;
149 
150  c.interior_value(u_var, qp, u);
151 
152  c.interior_curl(u_var, qp, curl_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 degrees of freedom.
158  for (unsigned int i=0; i != n_u_dofs; i++)
159  {
160  Fu(i) += (curl_u*curl_phi[i][qp] + u*phi[i][qp] - f*phi[i][qp])*JxW[qp];
161 
162  // Matrix contributions for the uu and vv couplings.
163  if (request_jacobian)
164  for (unsigned int j=0; j != n_u_dofs; j++)
165  Kuu(i,j) += (curl_phi[j][qp]*curl_phi[i][qp] +
166  phi[j][qp]*phi[i][qp])*JxW[qp];
167  }
168  } // end of the quadrature point qp-loop
169 
170  return request_jacobian;
171 }
172 
173 
174 bool CurlCurlSystem::side_time_derivative (bool request_jacobian,
175  DiffContext & context)
176 {
177  FEMContext & c = cast_ref<FEMContext &>(context);
178 
179  // Get finite element object
181  c.get_side_fe<RealGradient>(u_var, side_fe);
182 
183  // First we get some references to cell-specific data that
184  // will be used to assemble the linear system.
185 
186  // Element Jacobian * quadrature weights for interior integration
187  const std::vector<Real> & JxW = side_fe->get_JxW();
188 
189  // The velocity shape functions at interior quadrature points.
190  const std::vector<std::vector<RealGradient>> & phi = side_fe->get_phi();
191 
192  // The number of local degrees of freedom in each variable
193  const unsigned int n_u_dofs = c.get_dof_indices(u_var).size();
194 
195  const std::vector<Point> & normals = side_fe->get_normals();
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.get_elem_jacobian(u_var, 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  RealGradient N(normals[qp](0), normals[qp](1), normals[qp](2));
214 
215  Gradient u_exact = this->exact_solution(qpoint[qp]);
216 
217  Gradient Ncu = (u - u_exact).cross(N);
218 
219  for (unsigned int i=0; i != n_u_dofs; i++)
220  {
221  Fu(i) += penalty*Ncu*(phi[i][qp].cross(N))*JxW[qp];
222 
223  if (request_jacobian)
224  for (unsigned int j=0; j != n_u_dofs; j++)
225  Kuu(i,j) += penalty*(phi[j][qp].cross(N))*(phi[i][qp].cross(N))*JxW[qp];
226  }
227  }
228 
229  return request_jacobian;
230 }
231 
233 {
234  Real x = p(0);
235  Real y = p(1);
236  Real z = p(2);
237 
238  return soln.forcing(x, y, z);
239 }
240 
242 {
243  Real x = p(0);
244  Real y = p(1);
245  Real z = p(2);
246 
247  return soln(x, y, z);
248 }
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()