libMesh
kelly_error_estimator.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 
19 // C++ includes
20 #include <algorithm> // for std::fill
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for sqrt
23 
24 
25 // Local Includes
26 #include "libmesh/libmesh_common.h"
27 #include "libmesh/kelly_error_estimator.h"
28 #include "libmesh/error_vector.h"
29 #include "libmesh/fe_base.h"
30 #include "libmesh/libmesh_logging.h"
31 #include "libmesh/elem.h"
32 #include "libmesh/system.h"
33 
34 #include "libmesh/dense_vector.h"
35 #include "libmesh/tensor_tools.h"
36 
37 namespace libMesh
38 {
39 
40 
41 
42 void
44 {
45  const unsigned int n_vars = c.n_vars();
46  for (unsigned int v=0; v<n_vars; v++)
47  {
48  // Possibly skip this variable
49  if (error_norm.weight(v) == 0.0) continue;
50 
51  // FIXME: Need to generalize this to vector-valued elements. [PB]
52  FEBase * side_fe = libmesh_nullptr;
53 
54  const std::set<unsigned char> & elem_dims =
55  c.elem_dimensions();
56 
57  for (std::set<unsigned char>::const_iterator dim_it =
58  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
59  {
60  const unsigned char dim = *dim_it;
61 
62  fine_context->get_side_fe( v, side_fe, dim );
63 
64  // We'll need gradients on both sides for flux jump computation
65  side_fe->get_dphi();
66 
67  // But we only need normal vectors from one side
68  if (&c != coarse_context.get())
69  side_fe->get_normals();
70  }
71  }
72 }
73 
74 
75 
76 void
78 {
79  const Elem & coarse_elem = coarse_context->get_elem();
80  const Elem & fine_elem = fine_context->get_elem();
81 
82  FEBase * fe_fine = libmesh_nullptr;
83  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
84 
85  FEBase * fe_coarse = libmesh_nullptr;
86  coarse_context->get_side_fe( var, fe_coarse, fine_elem.dim() );
87 
88  Real error = 1.e-30;
89  unsigned int n_qp = fe_fine->n_quadrature_points();
90 
91  std::vector<std::vector<RealGradient>> dphi_coarse = fe_coarse->get_dphi();
92  std::vector<std::vector<RealGradient>> dphi_fine = fe_fine->get_dphi();
93  std::vector<Point> face_normals = fe_fine->get_normals();
94  std::vector<Real> JxW_face = fe_fine->get_JxW();
95 
96  for (unsigned int qp=0; qp != n_qp; ++qp)
97  {
98  // Calculate solution gradients on fine and coarse elements
99  // at this quadrature point
100  Gradient
101  grad_fine = fine_context->side_gradient(var, qp),
102  grad_coarse = coarse_context->side_gradient(var, qp);
103 
104  // Find the jump in the normal derivative
105  // at this quadrature point
106  const Number jump = (grad_fine - grad_coarse)*face_normals[qp];
107  const Real jump2 = TensorTools::norm_sq(jump);
108 
109  // Accumulate the jump integral
110  error += JxW_face[qp] * jump2;
111  }
112 
113  // Add the h-weighted jump integral to each error term
114  fine_error =
115  error * fine_elem.hmax() * error_norm.weight(var);
116  coarse_error =
117  error * coarse_elem.hmax() * error_norm.weight(var);
118 }
119 
120 
121 bool
123 {
124  const Elem & fine_elem = fine_context->get_elem();
125 
126  FEBase * fe_fine = libmesh_nullptr;
127  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
128 
129  const std::string & var_name =
130  fine_context->get_system().variable_name(var);
131 
132  std::vector<std::vector<RealGradient>> dphi_fine = fe_fine->get_dphi();
133  std::vector<Point> face_normals = fe_fine->get_normals();
134  std::vector<Real> JxW_face = fe_fine->get_JxW();
135  std::vector<Point> qface_point = fe_fine->get_xyz();
136 
137  // The reinitialization also recomputes the locations of
138  // the quadrature points on the side. By checking if the
139  // first quadrature point on the side is on a flux boundary
140  // for a particular variable, we will determine if the whole
141  // element is on a flux boundary (assuming quadrature points
142  // are strictly contained in the side).
143  if (this->_bc_function(fine_context->get_system(),
144  qface_point[0], var_name).first)
145  {
146  const Real h = fine_elem.hmax();
147 
148  // The number of quadrature points
149  const unsigned int n_qp = fe_fine->n_quadrature_points();
150 
151  // The error contribution from this face
152  Real error = 1.e-30;
153 
154  // loop over the integration points on the face.
155  for (unsigned int qp=0; qp<n_qp; qp++)
156  {
157  // Value of the imposed flux BC at this quadrature point.
158  const std::pair<bool,Real> flux_bc =
159  this->_bc_function(fine_context->get_system(),
160  qface_point[qp], var_name);
161 
162  // Be sure the BC function still thinks we're on the
163  // flux boundary.
164  libmesh_assert_equal_to (flux_bc.first, true);
165 
166  // The solution gradient from each element
167  Gradient grad_fine = fine_context->side_gradient(var, qp);
168 
169  // The difference between the desired BC and the approximate solution.
170  const Number jump = flux_bc.second - grad_fine*face_normals[qp];
171 
172  // The flux jump squared. If using complex numbers,
173  // TensorTools::norm_sq(z) returns |z|^2, where |z| is the modulus of z.
174  const Real jump2 = TensorTools::norm_sq(jump);
175 
176  // Integrate the error on the face. The error is
177  // scaled by an additional power of h, where h is
178  // the maximum side length for the element. This
179  // arises in the definition of the indicator.
180  error += JxW_face[qp]*jump2;
181 
182  } // End quadrature point loop
183 
184  fine_error = error*h*error_norm.weight(var);
185 
186  return true;
187  } // end if side on flux boundary
188  return false;
189 }
190 
191 
192 
193 void
194 KellyErrorEstimator::attach_flux_bc_function (std::pair<bool,Real> fptr(const System & system,
195  const Point & p,
196  const std::string & var_name))
197 {
198  _bc_function = fptr;
199 
200  // We may be turning boundary side integration on or off
201  if (fptr)
203  else
204  integrate_boundary_sides = false;
205 }
206 
207 } // namespace libMesh
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)
Pointer to function that provides BC information.
virtual void internal_side_integration() libmesh_override
The function which calculates a normal derivative jump based error term on an internal side...
SystemNorm error_norm
When estimating the error in a single system, the error_norm is used to control the scaling and norm ...
unsigned int dim
virtual unsigned int n_quadrature_points() const =0
Real weight(unsigned int var) const
Definition: system_norm.h:292
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
UniquePtr< FEMContext > fine_context
Context objects for integrating on the fine and coarse elements sharing a face.
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
void attach_flux_bc_function(std::pair< bool, Real > fptr(const System &system, const Point &p, const std::string &var_name))
Register a user function to use in computing the flux BCs.
Number fptr(const Point &p, const Parameters &, const std::string &libmesh_dbg_var(sys_name), const std::string &unknown_name)
Definition: projection.C:78
const std::vector< Point > & get_normals() const
Definition: fe_abstract.h:377
const std::set< unsigned char > & elem_dimensions() const
Definition: fem_context.h:913
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:61
Real fine_error
The fine and coarse error values to be set by each side_integration();.
virtual bool boundary_side_integration() libmesh_override
The function which calculates a normal derivative jump based error term on a boundary side...
virtual Real hmax() const
Definition: elem.C:475
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:216
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned int dim() const =0
unsigned int var
The variable number currently being evaluated.
bool integrate_boundary_sides
A boolean flag, by default false, to be set to true if integrations with boundary_side_integration() ...
UniquePtr< FEMContext > coarse_context
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
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 void init_context(FEMContext &c) libmesh_override
An initialization function, for requesting specific data from the FE objects.
unsigned int n_vars() const
Number of variables in solution.
Definition: diff_context.h:98