libMesh
fourth_error_estimators.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 #include "libmesh/libmesh_config.h"
20 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
21 
22 
23 // C++ includes
24 #include <algorithm> // for std::fill
25 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
26 #include <cmath> // for sqrt
27 
28 
29 // Local Includes
30 #include "libmesh/libmesh_common.h"
31 #include "libmesh/fourth_error_estimators.h"
32 #include "libmesh/error_vector.h"
33 #include "libmesh/fe_base.h"
34 #include "libmesh/libmesh_logging.h"
35 #include "libmesh/elem.h"
36 #include "libmesh/system.h"
37 
38 #include "libmesh/dense_vector.h"
39 #include "libmesh/tensor_tools.h"
40 
41 
42 namespace libMesh
43 {
44 
45 
46 void
48 {
49  const unsigned int n_vars = c.n_vars();
50  for (unsigned int v=0; v<n_vars; v++)
51  {
52  // Possibly skip this variable
53  if (error_norm.weight(v) == 0.0) continue;
54 
55  // FIXME: Need to generalize this to vector-valued elements. [PB]
56  FEBase * side_fe = libmesh_nullptr;
57 
58  const std::set<unsigned char> & elem_dims =
59  c.elem_dimensions();
60 
61  for (std::set<unsigned char>::const_iterator dim_it =
62  elem_dims.begin(); dim_it != elem_dims.end(); ++dim_it)
63  {
64  const unsigned char dim = *dim_it;
65 
66  fine_context->get_side_fe( v, side_fe, dim );
67 
68  // We'll need hessians on both sides for flux jump computation
69  side_fe->get_d2phi();
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  const DenseVector<Number> & Ucoarse = coarse_context->get_elem_solution();
83  const DenseVector<Number> & Ufine = fine_context->get_elem_solution();
84 
85  unsigned int dim = fine_elem.dim();
86 
87  FEBase * fe_fine = libmesh_nullptr;
88  fine_context->get_side_fe( var, fe_fine, dim );
89 
90  FEBase * fe_coarse = libmesh_nullptr;
91  coarse_context->get_side_fe( var, fe_coarse, dim );
92 
93  Real error = 1.e-30;
94  unsigned int n_qp = fe_fine->n_quadrature_points();
95 
96  std::vector<std::vector<RealTensor>> d2phi_coarse = fe_coarse->get_d2phi();
97  std::vector<std::vector<RealTensor>> d2phi_fine = fe_fine->get_d2phi();
98  std::vector<Real> JxW_face = fe_fine->get_JxW();
99 
100  for (unsigned int qp=0; qp != n_qp; ++qp)
101  {
102  // Calculate solution gradients on fine and coarse elements
103  // at this quadrature point
104  Number laplacian_fine = 0., laplacian_coarse = 0.;
105 
106  const unsigned int n_coarse_dofs = Ucoarse.size();
107  for (unsigned int i=0; i != n_coarse_dofs; ++i)
108  {
109  laplacian_coarse += d2phi_coarse[i][qp](0,0) * Ucoarse(i);
110  if (dim > 1)
111  laplacian_coarse += d2phi_coarse[i][qp](1,1) * Ucoarse(i);
112  if (dim > 2)
113  laplacian_coarse += d2phi_coarse[i][qp](2,2) * Ucoarse(i);
114  }
115 
116  const unsigned int n_fine_dofs = Ufine.size();
117  for (unsigned int i=0; i != n_fine_dofs; ++i)
118  {
119  laplacian_fine += d2phi_fine[i][qp](0,0) * Ufine(i);
120  if (dim > 1)
121  laplacian_fine += d2phi_fine[i][qp](1,1) * Ufine(i);
122  if (dim > 2)
123  laplacian_fine += d2phi_fine[i][qp](2,2) * Ufine(i);
124  }
125 
126 
127  // Find the jump in the Laplacian
128  // at this quadrature point
129  const Number jump = laplacian_fine - laplacian_coarse;
130  const Real jump2 = TensorTools::norm_sq(jump);
131 
132  // Accumulate the jump integral
133  error += JxW_face[qp] * jump2;
134  }
135 
136  // Add the h-weighted jump integral to each error term
137  fine_error =
138  error * fine_elem.hmax() * error_norm.weight(var);
139  coarse_error =
140  error * coarse_elem.hmax() * error_norm.weight(var);
141 }
142 
143 } // namespace libMesh
144 
145 #else // defined (LIBMESH_ENABLE_SECOND_DERIVATIVES)
146 
147 #include "libmesh/fourth_error_estimators.h"
148 
149 namespace libMesh
150 {
151 
152 void
154 {
155  libmesh_error_msg("Error: LaplacianErrorEstimator requires second " \
156  << "derivative support; try configuring libmesh with " \
157  << "--enable-second");
158 }
159 
160 
161 void
163 {
164  libmesh_error_msg("Error: LaplacianErrorEstimator requires second " \
165  << "derivative support; try configuring libmesh with " \
166  << "--enable-second");
167 }
168 
169 } // namespace libMesh
170 
171 #endif // defined (LIBMESH_ENABLE_SECOND_DERIVATIVES)
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.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
virtual void init_context(FEMContext &c) libmesh_override
An initialization function, for requesting specific data from the FE objects.
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
virtual void internal_side_integration() libmesh_override
The function which calculates a laplacian jump based error term on an internal side.
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 Real hmax() const
Definition: elem.C:475
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const std::vector< std::vector< OutputTensor > > & get_d2phi() const
Definition: fe_base.h:290
virtual unsigned int dim() const =0
unsigned int var
The variable number currently being evaluated.
UniquePtr< FEMContext > coarse_context
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