jump_error_estimator.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/error_vector.h"
30 #include "libmesh/fe_base.h"
31 #include "libmesh/fe_interface.h"
34 #include "libmesh/elem.h"
35 #include "libmesh/mesh_base.h"
36 #include "libmesh/system.h"
37 
38 #include "libmesh/dense_vector.h"
39 #include "libmesh/numeric_vector.h"
40 
41 namespace libMesh
42 {
43 
44 //-----------------------------------------------------------------
45 // JumpErrorEstimator implementations
47  ErrorVector&,
48  bool)
49 {
50 }
51 
52 
53 
55  ErrorVector& error_per_cell,
56  const NumericVector<Number>* solution_vector,
57  bool estimate_parent_error)
58 {
59  START_LOG("estimate_error()", "JumpErrorEstimator");
60  /*
61 
62  Conventions for assigning the direction of the normal:
63 
64  - e & f are global element ids
65 
66  Case (1.) Elements are at the same level, e<f
67  Compute the flux jump on the face and
68  add it as a contribution to error_per_cell[e]
69  and error_per_cell[f]
70 
71  ----------------------
72  | | |
73  | | f |
74  | | |
75  | e |---> n |
76  | | |
77  | | |
78  ----------------------
79 
80 
81  Case (2.) The neighbor is at a higher level.
82  Compute the flux jump on e's face and
83  add it as a contribution to error_per_cell[e]
84  and error_per_cell[f]
85 
86  ----------------------
87  | | | |
88  | | e |---> n |
89  | | | |
90  |-----------| f |
91  | | | |
92  | | | |
93  | | | |
94  ----------------------
95  */
96 
97  // The current mesh
98  const MeshBase& mesh = system.get_mesh();
99 
100  // The dimensionality of the mesh
101  const unsigned int dim = mesh.mesh_dimension();
102 
103  // The number of variables in the system
104  const unsigned int n_vars = system.n_vars();
105 
106  // The DofMap for this system
107  const DofMap& dof_map = system.get_dof_map();
108 
109  // Resize the error_per_cell vector to be
110  // the number of elements, initialize it to 0.
111  error_per_cell.resize (mesh.max_elem_id());
112  std::fill (error_per_cell.begin(), error_per_cell.end(), 0.);
113 
114  // Declare a vector of floats which is as long as
115  // error_per_cell above, and fill with zeros. This vector will be
116  // used to keep track of the number of edges (faces) on each active
117  // element which are either:
118  // 1) an internal edge
119  // 2) an edge on a Neumann boundary for which a boundary condition
120  // function has been specified.
121  // The error estimator can be scaled by the number of flux edges (faces)
122  // which the element actually has to obtain a more uniform measure
123  // of the error. Use floats instead of ints since in case 2 (above)
124  // f gets 1/2 of a flux face contribution from each of his
125  // neighbors
126  std::vector<float> n_flux_faces (error_per_cell.size());
127 
128  // Prepare current_local_solution to localize a non-standard
129  // solution vector if necessary
130  if (solution_vector && solution_vector != system.solution.get())
131  {
132  NumericVector<Number>* newsol =
133  const_cast<NumericVector<Number>*>(solution_vector);
134  System &sys = const_cast<System&>(system);
135  newsol->swap(*sys.solution);
136  sys.update();
137  }
138 
139  // Loop over all the variables in the system
140  for (var=0; var<n_vars; var++)
141  {
142  // Possibly skip this variable
143  if (error_norm.weight(var) == 0.0) continue;
144 
145  // The type of finite element to use for this variable
146  const FEType& fe_type = dof_map.variable_type (var);
147 
148  // Finite element objects for the same face from
149  // different sides
150  fe_fine = FEBase::build (dim, fe_type);
151  fe_coarse = FEBase::build (dim, fe_type);
152 
153  // Build an appropriate Gaussian quadrature rule
154  QGauss qrule (dim-1, fe_type.default_quadrature_order());
155 
156  // Tell the finite element for the fine element about the quadrature
157  // rule. The finite element for the coarse element need not know about it
158  fe_fine->attach_quadrature_rule (&qrule);
159 
160  // By convention we will always do the integration
161  // on the face of element e. We'll need its Jacobian values and
162  // physical point locations, at least
163  fe_fine->get_JxW();
164  fe_fine->get_xyz();
165 
166  // Our derived classes may want to do some initialization here
167  this->initialize(system, error_per_cell, estimate_parent_error);
168 
169  // The global DOF indices for elements e & f
170  std::vector<dof_id_type> dof_indices_fine;
171  std::vector<dof_id_type> dof_indices_coarse;
172 
173 
174 
175  // Iterate over all the active elements in the mesh
176  // that live on this processor.
179 
180  for (; elem_it != elem_end; ++elem_it)
181  {
182  // e is necessarily an active element on the local processor
183  const Elem* e = *elem_it;
184  const dof_id_type e_id = e->id();
185 
186 #ifdef LIBMESH_ENABLE_AMR
187  // See if the parent of element e has been examined yet;
188  // if not, we may want to compute the estimator on it
189  const Elem* parent = e->parent();
190 
191  // We only can compute and only need to compute on
192  // parents with all active children
193  bool compute_on_parent = true;
194  if (!parent || !estimate_parent_error)
195  compute_on_parent = false;
196  else
197  for (unsigned int c=0; c != parent->n_children(); ++c)
198  if (!parent->child(c)->active())
199  compute_on_parent = false;
200 
201  if (compute_on_parent &&
202  !error_per_cell[parent->id()])
203  {
204  // Compute a projection onto the parent
205  DenseVector<Number> Uparent;
207  dof_map, parent, Uparent,
208  var, false);
209 
210  // Loop over the neighbors of the parent
211  for (unsigned int n_p=0; n_p<parent->n_neighbors(); n_p++)
212  {
213  if (parent->neighbor(n_p) != NULL) // parent has a neighbor here
214  {
215  // Find the active neighbors in this direction
216  std::vector<const Elem*> active_neighbors;
217  parent->neighbor(n_p)->
218  active_family_tree_by_neighbor(active_neighbors,
219  parent);
220  // Compute the flux to each active neighbor
221  for (unsigned int a=0;
222  a != active_neighbors.size(); ++a)
223  {
224  const Elem *f = active_neighbors[a];
225  // FIXME - what about when f->level <
226  // parent->level()??
227  if (f->level() >= parent->level())
228  {
229  fine_elem = f;
230  coarse_elem = parent;
231  Ucoarse = Uparent;
232 
233  dof_map.dof_indices (fine_elem, dof_indices_fine, var);
234  const unsigned int n_dofs_fine =
235  libmesh_cast_int<unsigned int>(dof_indices_fine.size());
236  Ufine.resize(n_dofs_fine);
237 
238  for (unsigned int i=0; i<n_dofs_fine; i++)
239  Ufine(i) = system.current_solution(dof_indices_fine[i]);
240  this->reinit_sides();
242 
243  error_per_cell[fine_elem->id()] +=
244  static_cast<ErrorVectorReal>(fine_error);
245  error_per_cell[coarse_elem->id()] +=
246  static_cast<ErrorVectorReal>(coarse_error);
247 
248  // Keep track of the number of internal flux
249  // sides found on each element
250  n_flux_faces[fine_elem->id()]++;
251  n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
252  }
253  }
254  }
255  else if (integrate_boundary_sides)
256  {
257  fine_elem = parent;
258  Ufine = Uparent;
259 
260  // Reinitialize shape functions on the fine element side
261  fe_fine->reinit (fine_elem, fine_side);
262 
263  if (this->boundary_side_integration())
264  {
265  error_per_cell[fine_elem->id()] +=
266  static_cast<ErrorVectorReal>(fine_error);
267  n_flux_faces[fine_elem->id()]++;
268  }
269  }
270  }
271  }
272 #endif // #ifdef LIBMESH_ENABLE_AMR
273 
274  // If we do any more flux integration, e will be the fine element
275  fine_elem = e;
276 
277  // Loop over the neighbors of element e
278  for (unsigned int n_e=0; n_e<e->n_neighbors(); n_e++)
279  {
280  fine_side = n_e;
281 
282  if (e->neighbor(n_e) != NULL) // e is not on the boundary
283  {
284  const Elem* f = e->neighbor(n_e);
285  const dof_id_type f_id = f->id();
286 
287  // Compute flux jumps if we are in case 1 or case 2.
288  if ((f->active() && (f->level() == e->level()) && (e_id < f_id))
289  || (f->level() < e->level()))
290  {
291  // f is now the coarse element
292  coarse_elem = f;
293 
294  // Get the DOF indices for the two elements
295  dof_map.dof_indices (fine_elem, dof_indices_fine, var);
296  dof_map.dof_indices (coarse_elem, dof_indices_coarse, var);
297 
298  // The number of DOFS on each element
299  const unsigned int n_dofs_fine =
300  libmesh_cast_int<unsigned int>(dof_indices_fine.size());
301  const unsigned int n_dofs_coarse =
302  libmesh_cast_int<unsigned int>(dof_indices_coarse.size());
303  Ufine.resize(n_dofs_fine);
304  Ucoarse.resize(n_dofs_coarse);
305 
306  // The local solutions on each element
307  for (unsigned int i=0; i<n_dofs_fine; i++)
308  Ufine(i) = system.current_solution(dof_indices_fine[i]);
309  for (unsigned int i=0; i<n_dofs_coarse; i++)
310  Ucoarse(i) = system.current_solution(dof_indices_coarse[i]);
311 
312  this->reinit_sides();
314 
315  error_per_cell[fine_elem->id()] +=
316  static_cast<ErrorVectorReal>(fine_error);
317  error_per_cell[coarse_elem->id()] +=
318  static_cast<ErrorVectorReal>(coarse_error);
319 
320  // Keep track of the number of internal flux
321  // sides found on each element
322  n_flux_faces[fine_elem->id()]++;
323  n_flux_faces[coarse_elem->id()] += this->coarse_n_flux_faces_increment();
324  } // end if (case1 || case2)
325  } // if (e->neigbor(n_e) != NULL)
326 
327  // Otherwise, e is on the boundary. If it happens to
328  // be on a Dirichlet boundary, we need not do anything.
329  // On the other hand, if e is on a Neumann (flux) boundary
330  // with grad(u).n = g, we need to compute the additional residual
331  // (h * \int |g - grad(u_h).n|^2 dS)^(1/2).
332  // We can only do this with some knowledge of the boundary
333  // conditions, i.e. the user must have attached an appropriate
334  // BC function.
335  else
336  {
338  {
339  // Reinitialize shape functions on the fine element side
340  fe_fine->reinit (fine_elem, fine_side);
341 
342  // Get the DOF indices
343  dof_map.dof_indices (fine_elem, dof_indices_fine, var);
344 
345  // The number of DOFS on each element
346  const unsigned int n_dofs_fine =
347  libmesh_cast_int<unsigned int>(dof_indices_fine.size());
348  Ufine.resize(n_dofs_fine);
349 
350  for (unsigned int i=0; i<n_dofs_fine; i++)
351  Ufine(i) = system.current_solution(dof_indices_fine[i]);
352 
353  if (this->boundary_side_integration())
354  {
355  error_per_cell[fine_elem->id()] +=
356  static_cast<ErrorVectorReal>(fine_error);
357  n_flux_faces[fine_elem->id()]++;
358  }
359  } // end if _bc_function != NULL
360  } // end if (e->neighbor(n_e) == NULL)
361  } // end loop over neighbors
362  } // End loop over active local elements
363  } // End loop over variables
364 
365 
366 
367  // Each processor has now computed the error contribuions
368  // for its local elements. We need to sum the vector
369  // and then take the square-root of each component. Note
370  // that we only need to sum if we are running on multiple
371  // processors, and we only need to take the square-root
372  // if the value is nonzero. There will in general be many
373  // zeros for the inactive elements.
374 
375  // First sum the vector of estimated error values
376  this->reduce_error(error_per_cell, system.comm());
377 
378  // Compute the square-root of each component.
379  for (std::size_t i=0; i<error_per_cell.size(); i++)
380  if (error_per_cell[i] != 0.)
381  error_per_cell[i] = std::sqrt(error_per_cell[i]);
382 
383 
384  if (this->scale_by_n_flux_faces)
385  {
386  // Sum the vector of flux face counts
387  this->reduce_error(n_flux_faces, system.comm());
388 
389  // Sanity check: Make sure the number of flux faces is
390  // always an integer value
391 #ifdef DEBUG
392  for (unsigned int i=0; i<n_flux_faces.size(); ++i)
393  libmesh_assert_equal_to (n_flux_faces[i], static_cast<float>(static_cast<unsigned int>(n_flux_faces[i])) );
394 #endif
395 
396  // Scale the error by the number of flux faces for each element
397  for (unsigned int i=0; i<n_flux_faces.size(); ++i)
398  {
399  if (n_flux_faces[i] == 0.0) // inactive or non-local element
400  continue;
401 
402  //libMesh::out << "Element " << i << " has " << n_flux_faces[i] << " flux faces." << std::endl;
403  error_per_cell[i] /= static_cast<ErrorVectorReal>(n_flux_faces[i]);
404  }
405  }
406 
407  // If we used a non-standard solution before, now is the time to fix
408  // the current_local_solution
409  if (solution_vector && solution_vector != system.solution.get())
410  {
411  NumericVector<Number>* newsol =
412  const_cast<NumericVector<Number>*>(solution_vector);
413  System &sys = const_cast<System&>(system);
414  newsol->swap(*sys.solution);
415  sys.update();
416  }
417 
418  STOP_LOG("estimate_error()", "JumpErrorEstimator");
419 }
420 
421 
422 
423 void
425 {
426  // The master quadrature point locations on the coarse element
427  std::vector<Point> qp_coarse;
428 
429  // Reinitialize shape functions on the fine element side
430  fe_fine->reinit (fine_elem, fine_side);
431 
432  // Get the physical locations of the fine element quadrature points
433  std::vector<Point> qface_point = fe_fine->get_xyz();
434 
435  // Find their locations on the coarse element
437  coarse_elem, qface_point, qp_coarse);
438 
439  // Calculate the coarse element shape functions at those locations
440  fe_coarse->reinit (coarse_elem, &qp_coarse);
441 }
442 
443 
444 
446 {
447  // Keep track of the number of internal flux sides found on each
448  // element
449  unsigned int dim = coarse_elem->dim();
450 
451  const unsigned int divisor =
452  1 << (dim-1)*(fine_elem->level() - coarse_elem->level());
453 
454  // With a difference of n levels between fine and coarse elements,
455  // we compute a fractional flux face for the coarse element by adding:
456  // 1/2^n in 2D
457  // 1/4^n in 3D
458  // each time. This code will get hit 2^n times in 2D and 4^n
459  // times in 3D so that the final flux face count for the coarse
460  // element will be an integer value.
461 
462  return 1.0f / static_cast<float>(divisor);
463 }
464 
465 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo