libMesh
discontinuity_measure.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/discontinuity_measure.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 
38 namespace libMesh
39 {
40 
41 void
43 {
44  const unsigned int n_vars = c.n_vars();
45  for (unsigned int v=0; v<n_vars; v++)
46  {
47  // Possibly skip this variable
48  if (error_norm.weight(v) == 0.0) continue;
49 
50  // FIXME: Need to generalize this to vector-valued elements. [PB]
51  FEBase * side_fe = libmesh_nullptr;
52 
53  const std::set<unsigned char> & elem_dims =
54  c.elem_dimensions();
55 
56  for (std::set<unsigned char>::const_iterator dim_it =
57  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
58  {
59  const unsigned char dim = *dim_it;
60 
61  fine_context->get_side_fe( v, side_fe, dim );
62 
63  // We'll need values on both sides for discontinuity computation
64  side_fe->get_phi();
65  }
66  }
67 }
68 
69 
70 
71 void
73 {
74  const Elem & coarse_elem = coarse_context->get_elem();
75  const Elem & fine_elem = fine_context->get_elem();
76 
77  FEBase * fe_fine = libmesh_nullptr;
78  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
79 
80  FEBase * fe_coarse = libmesh_nullptr;
81  coarse_context->get_side_fe( var, fe_coarse, fine_elem.dim() );
82 
83  Real error = 1.e-30;
84  unsigned int n_qp = fe_fine->n_quadrature_points();
85 
86  std::vector<std::vector<Real>> phi_coarse = fe_coarse->get_phi();
87  std::vector<std::vector<Real>> phi_fine = fe_fine->get_phi();
88  std::vector<Real> JxW_face = fe_fine->get_JxW();
89 
90  for (unsigned int qp=0; qp != n_qp; ++qp)
91  {
92  // Calculate solution values on fine and coarse elements
93  // at this quadrature point
94  Number
95  u_fine = fine_context->side_value(var, qp),
96  u_coarse = coarse_context->side_value(var, qp);
97 
98  // Find the jump in the value
99  // at this quadrature point
100  const Number jump = u_fine - u_coarse;
101  const Real jump2 = TensorTools::norm_sq(jump);
102  // Accumulate the jump integral
103  error += JxW_face[qp] * jump2;
104  }
105 
106  // Add the h-weighted jump integral to each error term
107  fine_error =
108  error * fine_elem.hmax() * error_norm.weight(var);
109  coarse_error =
110  error * coarse_elem.hmax() * error_norm.weight(var);
111 }
112 
113 
114 bool
116 {
117  const Elem & fine_elem = fine_context->get_elem();
118 
119  FEBase * fe_fine = libmesh_nullptr;
120  fine_context->get_side_fe( var, fe_fine, fine_elem.dim() );
121 
122  const std::string & var_name =
123  fine_context->get_system().variable_name(var);
124 
125  std::vector<std::vector<Real>> phi_fine = fe_fine->get_phi();
126  std::vector<Real> JxW_face = fe_fine->get_JxW();
127  std::vector<Point> qface_point = fe_fine->get_xyz();
128 
129  // The reinitialization also recomputes the locations of
130  // the quadrature points on the side. By checking if the
131  // first quadrature point on the side is on an essential boundary
132  // for a particular variable, we will determine if the whole
133  // element is on an essential boundary (assuming quadrature points
134  // are strictly contained in the side).
135  if (this->_bc_function(fine_context->get_system(),
136  qface_point[0], var_name).first)
137  {
138  const Real h = fine_elem.hmax();
139 
140  // The number of quadrature points
141  const unsigned int n_qp = fe_fine->n_quadrature_points();
142 
143  // The error contribution from this face
144  Real error = 1.e-30;
145 
146  // loop over the integration points on the face.
147  for (unsigned int qp=0; qp<n_qp; qp++)
148  {
149  // Value of the imposed essential BC at this quadrature point.
150  const std::pair<bool,Real> essential_bc =
151  this->_bc_function(fine_context->get_system(), qface_point[qp],
152  var_name);
153 
154  // Be sure the BC function still thinks we're on the
155  // essential boundary.
156  libmesh_assert_equal_to (essential_bc.first, true);
157 
158  // The solution value on each point
159  Number u_fine = fine_context->side_value(var, qp);
160 
161  // The difference between the desired BC and the approximate solution.
162  const Number jump = essential_bc.second - u_fine;
163 
164  // The flux jump squared. If using complex numbers,
165  // norm_sq(z) returns |z|^2, where |z| is the modulus of z.
166  const Real jump2 = TensorTools::norm_sq(jump);
167 
168  // Integrate the error on the face. The error is
169  // scaled by an additional power of h, where h is
170  // the maximum side length for the element. This
171  // arises in the definition of the indicator.
172  error += JxW_face[qp]*jump2;
173 
174  } // End quadrature point loop
175 
176  fine_error = error*h*error_norm.weight(var);
177 
178  return true;
179  } // end if side on flux boundary
180  return false;
181 }
182 
183 
184 
185 void
187  const Point & p,
188  const std::string & var_name))
189 {
190  _bc_function = fptr;
191 
192  // We may be turning boundary side integration on or off
193  if (fptr)
195  else
196  integrate_boundary_sides = false;
197 }
198 
199 } // namespace libMesh
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
void attach_essential_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 essential BCs.
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.
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::set< unsigned char > & elem_dimensions() const
Definition: fem_context.h:913
const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:237
std::pair< bool, Real >(* _bc_function)(const System &system, const Point &p, const std::string &var_name)
Pointer to function that provides BC information.
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
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 internal_side_integration() libmesh_override
The function which calculates a normal derivative jump based error term on an internal side...
virtual void init_context(FEMContext &c) libmesh_override
An initialization function, for requesting specific data from the FE objects.
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
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.
unsigned int n_vars() const
Number of variables in solution.
Definition: diff_context.h:98