libMesh::HPCoarsenTest Class Reference

#include <hp_coarsentest.h>

Inheritance diagram for libMesh::HPCoarsenTest:

Public Member Functions

 HPCoarsenTest ()
 
virtual ~HPCoarsenTest ()
 
virtual void select_refinement (System &system) libmesh_override
 

Public Attributes

Real p_weight
 
std::vector< float > component_scale
 

Protected Member Functions

void add_projection (const System &, const Elem *, unsigned int var)
 

Protected Attributes

const Elemcoarse
 
std::vector< dof_id_typedof_indices
 
UniquePtr< FEBasefe
 
UniquePtr< FEBasefe_coarse
 
const std::vector< std::vector< Real > > * phi
 
const std::vector< std::vector< Real > > * phi_coarse
 
const std::vector< std::vector< RealGradient > > * dphi
 
const std::vector< std::vector< RealGradient > > * dphi_coarse
 
const std::vector< std::vector< RealTensor > > * d2phi
 
const std::vector< std::vector< RealTensor > > * d2phi_coarse
 
const std::vector< Real > * JxW
 
const std::vector< Point > * xyz_values
 
std::vector< Pointcoarse_qpoints
 
UniquePtr< QBaseqrule
 
DenseMatrix< NumberKe
 
DenseVector< NumberFe
 
DenseVector< NumberUc
 
DenseVector< NumberUp
 

Detailed Description

This class uses the error estimate given by different types of derefinement (h coarsening or p reduction) to choose between h refining and p elevation. Currently we assume that a set of elements has already been flagged for h refinement, and we may want to change some of those elements to be flagged for p refinement.

This code is currently experimental and will not produce optimal hp meshes without significant improvement.

Author
Roy H. Stogner
Date
2006

Definition at line 68 of file hp_coarsentest.h.

Constructor & Destructor Documentation

libMesh::HPCoarsenTest::HPCoarsenTest ( )
inline

Constructor.

Definition at line 75 of file hp_coarsentest.h.

75  : p_weight(1.0)
76  {
77  libmesh_experimental();
78  }
virtual libMesh::HPCoarsenTest::~HPCoarsenTest ( )
inlinevirtual

Destructor.

Definition at line 83 of file hp_coarsentest.h.

References select_refinement().

83 {}

Member Function Documentation

void libMesh::HPCoarsenTest::add_projection ( const System system,
const Elem elem,
unsigned int  var 
)
protected

The helper function which adds individual fine element data to the coarse element projection

Definition at line 46 of file hp_coarsentest.C.

References libMesh::Elem::active(), libMesh::TypeVector< T >::add_scaled(), libMesh::TypeTensor< T >::add_scaled(), libMesh::C_ONE, libMesh::C_ZERO, libMesh::Elem::child_ptr(), coarse, coarse_qpoints, libMesh::TypeTensor< T >::contract(), libMesh::System::current_solution(), d2phi, d2phi_coarse, dof_indices, libMesh::DofMap::dof_indices(), dphi, fe, Fe, fe_coarse, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::FEInterface::inverse_map(), Ke, libMesh::libmesh_assert(), libMesh::MeshBase::mesh_dimension(), libMesh::Elem::n_children(), phi_coarse, qrule, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseVector< T >::size(), libMesh::Elem::subactive(), Uc, libMesh::DofMap::variable_type(), xyz_values, libMesh::DenseMatrix< T >::zero(), libMesh::DenseVector< T >::zero(), and libMesh::zero.

Referenced by select_refinement().

49 {
50  // If we have children, we need to add their projections instead
51  if (!elem->active())
52  {
53  libmesh_assert(!elem->subactive());
54  for (unsigned int c = 0; c != elem->n_children(); ++c)
55  this->add_projection(system, elem->child_ptr(c), var);
56  return;
57  }
58 
59  // The DofMap for this system
60  const DofMap & dof_map = system.get_dof_map();
61 
62  // The type of finite element to use for this variable
63  const FEType & fe_type = dof_map.variable_type (var);
64 
65  const FEContinuity cont = fe->get_continuity();
66 
67  fe->reinit(elem);
68 
69  dof_map.dof_indices(elem, dof_indices, var);
70 
71  const unsigned int n_dofs =
72  cast_int<unsigned int>(dof_indices.size());
73 
74  FEInterface::inverse_map (system.get_mesh().mesh_dimension(),
75  fe_type, coarse, *xyz_values, coarse_qpoints);
76 
77  fe_coarse->reinit(coarse, &coarse_qpoints);
78 
79  const unsigned int n_coarse_dofs =
80  cast_int<unsigned int>(phi_coarse->size());
81 
82  if (Uc.size() == 0)
83  {
84  Ke.resize(n_coarse_dofs, n_coarse_dofs);
85  Ke.zero();
86  Fe.resize(n_coarse_dofs);
87  Fe.zero();
88  Uc.resize(n_coarse_dofs);
89  Uc.zero();
90  }
91  libmesh_assert_equal_to (Uc.size(), phi_coarse->size());
92 
93  // Loop over the quadrature points
94  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
95  {
96  // The solution value at the quadrature point
97  Number val = libMesh::zero;
98  Gradient grad;
99  Tensor hess;
100 
101  for (unsigned int i=0; i != n_dofs; i++)
102  {
103  dof_id_type dof_num = dof_indices[i];
104  val += (*phi)[i][qp] *
105  system.current_solution(dof_num);
106  if (cont == C_ZERO || cont == C_ONE)
107  grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
108  // grad += (*dphi)[i][qp] *
109  // system.current_solution(dof_num);
110  if (cont == C_ONE)
111  hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
112  // hess += (*d2phi)[i][qp] *
113  // system.current_solution(dof_num);
114  }
115 
116  // The projection matrix and vector
117  for (unsigned int i=0; i != Fe.size(); ++i)
118  {
119  Fe(i) += (*JxW)[qp] *
120  (*phi_coarse)[i][qp]*val;
121  if (cont == C_ZERO || cont == C_ONE)
122  Fe(i) += (*JxW)[qp] *
123  (grad*(*dphi_coarse)[i][qp]);
124  if (cont == C_ONE)
125  Fe(i) += (*JxW)[qp] *
126  hess.contract((*d2phi_coarse)[i][qp]);
127  // Fe(i) += (*JxW)[qp] *
128  // (*d2phi_coarse)[i][qp].contract(hess);
129 
130  for (unsigned int j=0; j != Fe.size(); ++j)
131  {
132  Ke(i,j) += (*JxW)[qp] *
133  (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
134  if (cont == C_ZERO || cont == C_ONE)
135  Ke(i,j) += (*JxW)[qp] *
136  (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
137  if (cont == C_ONE)
138  Ke(i,j) += (*JxW)[qp] *
139  ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
140  }
141  }
142  }
143 }
void add_projection(const System &, const Elem *, unsigned int var)
virtual void zero() libmesh_override
Definition: dense_matrix.h:792
DenseVector< Number > Fe
const std::vector< std::vector< Real > > * phi_coarse
void resize(const unsigned int n)
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Definition: dense_vector.h:374
UniquePtr< QBase > qrule
UniquePtr< FEBase > fe
const Number zero
Definition: libmesh.h:178
const std::vector< Point > * xyz_values
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
DenseVector< Number > Uc
libmesh_assert(j)
NumberVectorValue Gradient
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:556
UniquePtr< FEBase > fe_coarse
NumberTensorValue Tensor
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
const std::vector< std::vector< RealTensor > > * d2phi_coarse
const std::vector< std::vector< RealTensor > > * d2phi
const std::vector< std::vector< RealGradient > > * dphi
std::vector< dof_id_type > dof_indices
std::vector< Point > coarse_qpoints
uint8_t dof_id_type
Definition: id_types.h:64
DenseMatrix< Number > Ke
void libMesh::HPCoarsenTest::select_refinement ( System system)
virtual

This pure virtual function must be redefined in derived classes to take a mesh flagged for h refinement and potentially change the desired refinement type.

Implements libMesh::HPSelector.

Definition at line 145 of file hp_coarsentest.C.

References libMesh::MeshBase::active_local_elements_begin(), libMesh::MeshBase::active_local_elements_end(), add_projection(), libMesh::TypeVector< T >::add_scaled(), libMesh::TypeTensor< T >::add_scaled(), libMesh::FEGenericBase< OutputType >::build(), libMesh::C_ONE, libMesh::C_ZERO, libMesh::DenseMatrix< T >::cholesky_solve(), coarse, coarse_qpoints, libMesh::HPSelector::component_scale, libMesh::TypeTensor< T >::contract(), libMesh::System::current_solution(), d2phi, d2phi_coarse, libMesh::FEType::default_quadrature_rule(), libMesh::DISCONTINUOUS, libMesh::Elem::DO_NOTHING, dof_indices, libMesh::DofMap::dof_indices(), libMesh::DofObject::dof_number(), dphi, dphi_coarse, libMesh::ErrorVectorReal, fe, Fe, fe_coarse, libMesh::System::get_dof_map(), libMesh::System::get_mesh(), libMesh::DofObject::id(), libMesh::FEInterface::inverse_map(), libMesh::Elem::is_vertex(), JxW, Ke, libMesh::libmesh_assert(), libmesh_nullptr, std::max(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::Elem::n_children(), libMesh::FEInterface::n_dofs(), n_nodes, libMesh::Elem::n_nodes(), n_vars, libMesh::System::n_vars(), libMesh::Elem::node_ref(), libMesh::TensorTools::norm_sq(), libMesh::TypeVector< T >::norm_sq(), libMesh::TypeTensor< T >::norm_sq(), libMesh::System::number(), libMesh::FEType::order, libMesh::Elem::p_level(), p_weight, libMesh::Elem::parent(), phi, phi_coarse, qrule, libMesh::Real, libMesh::Elem::REFINE, libMesh::Elem::refinement_flag(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::Elem::set_p_refinement_flag(), libMesh::Elem::set_refinement_flag(), libMesh::DenseVector< T >::size(), libMesh::TypeVector< T >::subtract_scaled(), libMesh::TypeTensor< T >::subtract_scaled(), libMesh::Elem::type(), Uc, Up, libMesh::DofMap::variable_type(), xyz_values, libMesh::DenseMatrix< T >::zero(), libMesh::DenseVector< T >::zero(), and libMesh::zero.

Referenced by ~HPCoarsenTest().

146 {
147  LOG_SCOPE("select_refinement()", "HPCoarsenTest");
148 
149  // The current mesh
150  MeshBase & mesh = system.get_mesh();
151 
152  // The dimensionality of the mesh
153  const unsigned int dim = mesh.mesh_dimension();
154 
155  // The number of variables in the system
156  const unsigned int n_vars = system.n_vars();
157 
158  // The DofMap for this system
159  const DofMap & dof_map = system.get_dof_map();
160 
161  // The system number (for doing bad hackery)
162  const unsigned int sys_num = system.number();
163 
164  // Check for a valid component_scale
165  if (!component_scale.empty())
166  {
167  if (component_scale.size() != n_vars)
168  libmesh_error_msg("ERROR: component_scale is the wrong size:\n" \
169  << " component_scale.size()=" \
170  << component_scale.size() \
171  << "\n n_vars=" \
172  << n_vars);
173  }
174  else
175  {
176  // No specified scaling. Scale all variables by one.
177  component_scale.resize (n_vars, 1.0);
178  }
179 
180  // Resize the error_per_cell vectors to handle
181  // the number of elements, initialize them to 0.
182  std::vector<ErrorVectorReal> h_error_per_cell(mesh.max_elem_id(), 0.);
183  std::vector<ErrorVectorReal> p_error_per_cell(mesh.max_elem_id(), 0.);
184 
185  // Loop over all the variables in the system
186  for (unsigned int var=0; var<n_vars; var++)
187  {
188  // Possibly skip this variable
189  if (!component_scale.empty())
190  if (component_scale[var] == 0.0) continue;
191 
192  // The type of finite element to use for this variable
193  const FEType & fe_type = dof_map.variable_type (var);
194 
195  // Finite element objects for a fine (and probably a coarse)
196  // element will be needed
197  fe = FEBase::build (dim, fe_type);
198  fe_coarse = FEBase::build (dim, fe_type);
199 
200  // Any cached coarse element results have expired
202  unsigned int cached_coarse_p_level = 0;
203 
204  const FEContinuity cont = fe->get_continuity();
205  libmesh_assert (cont == DISCONTINUOUS || cont == C_ZERO ||
206  cont == C_ONE);
207 
208  // Build an appropriate quadrature rule
209  qrule = fe_type.default_quadrature_rule(dim);
210 
211  // Tell the refined finite element about the quadrature
212  // rule. The coarse finite element need not know about it
213  fe->attach_quadrature_rule (qrule.get());
214 
215  // We will always do the integration
216  // on the fine elements. Get their Jacobian values, etc..
217  JxW = &(fe->get_JxW());
218  xyz_values = &(fe->get_xyz());
219 
220  // The shape functions
221  phi = &(fe->get_phi());
222  phi_coarse = &(fe_coarse->get_phi());
223 
224  // The shape function derivatives
225  if (cont == C_ZERO || cont == C_ONE)
226  {
227  dphi = &(fe->get_dphi());
228  dphi_coarse = &(fe_coarse->get_dphi());
229  }
230 
231 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
232  // The shape function second derivatives
233  if (cont == C_ONE)
234  {
235  d2phi = &(fe->get_d2phi());
236  d2phi_coarse = &(fe_coarse->get_d2phi());
237  }
238 #endif // defined (LIBMESH_ENABLE_SECOND_DERIVATIVES)
239 
240  // Iterate over all the active elements in the mesh
241  // that live on this processor.
242 
243  MeshBase::const_element_iterator elem_it =
244  mesh.active_local_elements_begin();
245  const MeshBase::const_element_iterator elem_end =
246  mesh.active_local_elements_end();
247 
248  for (; elem_it != elem_end; ++elem_it)
249  {
250  const Elem * elem = *elem_it;
251 
252  // We're only checking elements that are already flagged for h
253  // refinement
254  if (elem->refinement_flag() != Elem::REFINE)
255  continue;
256 
257  const dof_id_type e_id = elem->id();
258 
259  // Find the projection onto the parent element,
260  // if necessary
261  if (elem->parent() &&
262  (coarse != elem->parent() ||
263  cached_coarse_p_level != elem->p_level()))
264  {
265  Uc.resize(0);
266 
267  coarse = elem->parent();
268  cached_coarse_p_level = elem->p_level();
269 
270  unsigned int old_parent_level = coarse->p_level();
271  (const_cast<Elem *>(coarse))->hack_p_level(elem->p_level());
272 
273  this->add_projection(system, coarse, var);
274 
275  (const_cast<Elem *>(coarse))->hack_p_level(old_parent_level);
276 
277  // Solve the h-coarsening projection problem
279  }
280 
281  fe->reinit(elem);
282 
283  // Get the DOF indices for the fine element
284  dof_map.dof_indices (elem, dof_indices, var);
285 
286  // The number of quadrature points
287  const unsigned int n_qp = qrule->n_points();
288 
289  // The number of DOFS on the fine element
290  const unsigned int n_dofs =
291  cast_int<unsigned int>(dof_indices.size());
292 
293  // The number of nodes on the fine element
294  const unsigned int n_nodes = elem->n_nodes();
295 
296  // The average element value (used as an ugly hack
297  // when we have nothing p-coarsened to compare to)
298  // Real average_val = 0.;
299  Number average_val = 0.;
300 
301  // Calculate this variable's contribution to the p
302  // refinement error
303 
304  if (elem->p_level() == 0)
305  {
306  unsigned int n_vertices = 0;
307  for (unsigned int n = 0; n != n_nodes; ++n)
308  if (elem->is_vertex(n))
309  {
310  n_vertices++;
311  const Node & node = elem->node_ref(n);
312  average_val += system.current_solution
313  (node.dof_number(sys_num,var,0));
314  }
315  average_val /= n_vertices;
316  }
317  else
318  {
319  unsigned int old_elem_level = elem->p_level();
320  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level - 1);
321 
322  fe_coarse->reinit(elem, &(qrule->get_points()));
323 
324  const unsigned int n_coarse_dofs =
325  cast_int<unsigned int>(phi_coarse->size());
326 
327  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
328 
329  Ke.resize(n_coarse_dofs, n_coarse_dofs);
330  Ke.zero();
331  Fe.resize(n_coarse_dofs);
332  Fe.zero();
333 
334  // Loop over the quadrature points
335  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
336  {
337  // The solution value at the quadrature point
338  Number val = libMesh::zero;
339  Gradient grad;
340  Tensor hess;
341 
342  for (unsigned int i=0; i != n_dofs; i++)
343  {
344  dof_id_type dof_num = dof_indices[i];
345  val += (*phi)[i][qp] *
346  system.current_solution(dof_num);
347  if (cont == C_ZERO || cont == C_ONE)
348  grad.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
349  // grad += (*dphi)[i][qp] *
350  // system.current_solution(dof_num);
351  if (cont == C_ONE)
352  hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
353  // hess += (*d2phi)[i][qp] *
354  // system.current_solution(dof_num);
355  }
356 
357  // The projection matrix and vector
358  for (unsigned int i=0; i != Fe.size(); ++i)
359  {
360  Fe(i) += (*JxW)[qp] *
361  (*phi_coarse)[i][qp]*val;
362  if (cont == C_ZERO || cont == C_ONE)
363  Fe(i) += (*JxW)[qp] *
364  grad * (*dphi_coarse)[i][qp];
365  if (cont == C_ONE)
366  Fe(i) += (*JxW)[qp] *
367  hess.contract((*d2phi_coarse)[i][qp]);
368 
369  for (unsigned int j=0; j != Fe.size(); ++j)
370  {
371  Ke(i,j) += (*JxW)[qp] *
372  (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
373  if (cont == C_ZERO || cont == C_ONE)
374  Ke(i,j) += (*JxW)[qp] *
375  (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
376  if (cont == C_ONE)
377  Ke(i,j) += (*JxW)[qp] *
378  ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
379  }
380  }
381  }
382 
383  // Solve the p-coarsening projection problem
385  }
386 
387  // loop over the integration points on the fine element
388  for (unsigned int qp=0; qp<n_qp; qp++)
389  {
390  Number value_error = 0.;
391  Gradient grad_error;
392  Tensor hessian_error;
393  for (unsigned int i=0; i<n_dofs; i++)
394  {
395  const dof_id_type dof_num = dof_indices[i];
396  value_error += (*phi)[i][qp] *
397  system.current_solution(dof_num);
398  if (cont == C_ZERO || cont == C_ONE)
399  grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
400  // grad_error += (*dphi)[i][qp] *
401  // system.current_solution(dof_num);
402  if (cont == C_ONE)
403  hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
404  // hessian_error += (*d2phi)[i][qp] *
405  // system.current_solution(dof_num);
406  }
407  if (elem->p_level() == 0)
408  {
409  value_error -= average_val;
410  }
411  else
412  {
413  for (unsigned int i=0; i<Up.size(); i++)
414  {
415  value_error -= (*phi_coarse)[i][qp] * Up(i);
416  if (cont == C_ZERO || cont == C_ONE)
417  grad_error.subtract_scaled((*dphi_coarse)[i][qp], Up(i));
418  // grad_error -= (*dphi_coarse)[i][qp] * Up(i);
419  if (cont == C_ONE)
420  hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Up(i));
421  // hessian_error -= (*d2phi_coarse)[i][qp] * Up(i);
422  }
423  }
424 
425  p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
426  (component_scale[var] *
427  (*JxW)[qp] * TensorTools::norm_sq(value_error));
428  if (cont == C_ZERO || cont == C_ONE)
429  p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
430  (component_scale[var] *
431  (*JxW)[qp] * grad_error.norm_sq());
432  if (cont == C_ONE)
433  p_error_per_cell[e_id] += static_cast<ErrorVectorReal>
434  (component_scale[var] *
435  (*JxW)[qp] * hessian_error.norm_sq());
436  }
437 
438  // Calculate this variable's contribution to the h
439  // refinement error
440 
441  if (!elem->parent())
442  {
443  // For now, we'll always start with an h refinement
444  h_error_per_cell[e_id] =
446  }
447  else
448  {
449  FEInterface::inverse_map (dim, fe_type, coarse,
451 
452  unsigned int old_parent_level = coarse->p_level();
453  (const_cast<Elem *>(coarse))->hack_p_level(elem->p_level());
454 
455  fe_coarse->reinit(coarse, &coarse_qpoints);
456 
457  (const_cast<Elem *>(coarse))->hack_p_level(old_parent_level);
458 
459  // The number of DOFS on the coarse element
460  unsigned int n_coarse_dofs =
461  cast_int<unsigned int>(phi_coarse->size());
462 
463  // Loop over the quadrature points
464  for (unsigned int qp=0; qp<n_qp; qp++)
465  {
466  // The solution difference at the quadrature point
467  Number value_error = libMesh::zero;
468  Gradient grad_error;
469  Tensor hessian_error;
470 
471  for (unsigned int i=0; i != n_dofs; ++i)
472  {
473  const dof_id_type dof_num = dof_indices[i];
474  value_error += (*phi)[i][qp] *
475  system.current_solution(dof_num);
476  if (cont == C_ZERO || cont == C_ONE)
477  grad_error.add_scaled((*dphi)[i][qp], system.current_solution(dof_num));
478  // grad_error += (*dphi)[i][qp] *
479  // system.current_solution(dof_num);
480  if (cont == C_ONE)
481  hessian_error.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
482  // hessian_error += (*d2phi)[i][qp] *
483  // system.current_solution(dof_num);
484  }
485 
486  for (unsigned int i=0; i != n_coarse_dofs; ++i)
487  {
488  value_error -= (*phi_coarse)[i][qp] * Uc(i);
489  if (cont == C_ZERO || cont == C_ONE)
490  // grad_error -= (*dphi_coarse)[i][qp] * Uc(i);
491  grad_error.subtract_scaled((*dphi_coarse)[i][qp], Uc(i));
492  if (cont == C_ONE)
493  hessian_error.subtract_scaled((*d2phi_coarse)[i][qp], Uc(i));
494  // hessian_error -= (*d2phi_coarse)[i][qp] * Uc(i);
495  }
496 
497  h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
498  (component_scale[var] *
499  (*JxW)[qp] * TensorTools::norm_sq(value_error));
500  if (cont == C_ZERO || cont == C_ONE)
501  h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
502  (component_scale[var] *
503  (*JxW)[qp] * grad_error.norm_sq());
504  if (cont == C_ONE)
505  h_error_per_cell[e_id] += static_cast<ErrorVectorReal>
506  (component_scale[var] *
507  (*JxW)[qp] * hessian_error.norm_sq());
508  }
509 
510  }
511  }
512  }
513 
514  // Now that we've got our approximations for p_error and h_error, let's see
515  // if we want to switch any h refinement flags to p refinement
516 
517  // Iterate over all the active elements in the mesh
518  // that live on this processor.
519 
520  MeshBase::element_iterator elem_it =
521  mesh.active_local_elements_begin();
522  const MeshBase::element_iterator elem_end =
523  mesh.active_local_elements_end();
524 
525  for (; elem_it != elem_end; ++elem_it)
526  {
527  Elem * elem = *elem_it;
528 
529  // We're only checking elements that are already flagged for h
530  // refinement
531  if (elem->refinement_flag() != Elem::REFINE)
532  continue;
533 
534  const dof_id_type e_id = elem->id();
535 
536  unsigned int dofs_per_elem = 0, dofs_per_p_elem = 0;
537 
538  // Loop over all the variables in the system
539  for (unsigned int var=0; var<n_vars; var++)
540  {
541  // The type of finite element to use for this variable
542  const FEType & fe_type = dof_map.variable_type (var);
543 
544  // FIXME: we're overestimating the number of DOFs added by h
545  // refinement
546  FEType elem_fe_type = fe_type;
547  elem_fe_type.order =
548  static_cast<Order>(fe_type.order + elem->p_level());
549  dofs_per_elem +=
550  FEInterface::n_dofs(dim, elem_fe_type, elem->type());
551 
552  elem_fe_type.order =
553  static_cast<Order>(fe_type.order + elem->p_level() + 1);
554  dofs_per_p_elem +=
555  FEInterface::n_dofs(dim, elem_fe_type, elem->type());
556  }
557 
558  const unsigned int new_h_dofs = dofs_per_elem *
559  (elem->n_children() - 1);
560 
561  const unsigned int new_p_dofs = dofs_per_p_elem -
562  dofs_per_elem;
563 
564  /*
565  libMesh::err << "Cell " << e_id << ": h = " << elem->hmax()
566  << ", p = " << elem->p_level() + 1 << "," << std::endl
567  << " h_error = " << h_error_per_cell[e_id]
568  << ", p_error = " << p_error_per_cell[e_id] << std::endl
569  << " new_h_dofs = " << new_h_dofs
570  << ", new_p_dofs = " << new_p_dofs << std::endl;
571  */
572  const Real p_value =
573  std::sqrt(p_error_per_cell[e_id]) * p_weight / new_p_dofs;
574  const Real h_value =
575  std::sqrt(h_error_per_cell[e_id]) /
576  static_cast<Real>(new_h_dofs);
577  if (p_value > h_value)
578  {
579  elem->set_p_refinement_flag(Elem::REFINE);
580  elem->set_refinement_flag(Elem::DO_NOTHING);
581  }
582  }
583 }
const std::vector< std::vector< Real > > * phi
std::vector< float > component_scale
Definition: hp_selector.h:78
void add_projection(const System &, const Elem *, unsigned int var)
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:414
virtual void zero() libmesh_override
Definition: dense_matrix.h:792
DenseVector< Number > Fe
unsigned int p_level() const
Definition: elem.h:2224
const std::vector< std::vector< Real > > * phi_coarse
void resize(const unsigned int n)
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Definition: dense_vector.h:374
const Elem * parent() const
Definition: elem.h:2148
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
const std::vector< std::vector< RealGradient > > * dphi_coarse
UniquePtr< QBase > qrule
const unsigned int n_vars
Definition: tecplot_io.C:68
UniquePtr< FEBase > fe
const Number zero
Definition: libmesh.h:178
const std::vector< Point > * xyz_values
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
long double max(long double a, double b)
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
DenseVector< Number > Uc
libmesh_assert(j)
const std::vector< Real > * JxW
const dof_id_type n_nodes
Definition: tecplot_io.C:67
NumberVectorValue Gradient
static Point inverse_map(const unsigned int dim, const FEType &fe_t, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_interface.C:556
DenseVector< Number > Up
UniquePtr< FEBase > fe_coarse
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
NumberTensorValue Tensor
void resize(const unsigned int new_m, const unsigned int new_n)
Definition: dense_matrix.h:776
const std::vector< std::vector< RealTensor > > * d2phi_coarse
const std::vector< std::vector< RealTensor > > * d2phi
const std::vector< std::vector< RealGradient > > * dphi
std::vector< dof_id_type > dof_indices
std::vector< Point > coarse_qpoints
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
Definition: dense_matrix.C:911
uint8_t dof_id_type
Definition: id_types.h:64
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
DenseMatrix< Number > Ke

Member Data Documentation

const Elem* libMesh::HPCoarsenTest::coarse
protected

The coarse element on which a solution projection is cached

Definition at line 110 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

std::vector<Point> libMesh::HPCoarsenTest::coarse_qpoints
protected

Definition at line 138 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

std::vector<float> libMesh::HPSelector::component_scale
inherited

This vector can be used to "scale" certain variables in a system. If the mask is not empty, the consideration given to each component's h and p error estimates will be scaled by component_scale[c].

Definition at line 78 of file hp_selector.h.

Referenced by select_refinement().

const std::vector<std::vector<RealTensor> >* libMesh::HPCoarsenTest::d2phi
protected

Definition at line 127 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

const std::vector<std::vector<RealTensor> > * libMesh::HPCoarsenTest::d2phi_coarse
protected

Definition at line 127 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

std::vector<dof_id_type> libMesh::HPCoarsenTest::dof_indices
protected

Global DOF indices for fine elements

Definition at line 115 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

const std::vector<std::vector<RealGradient> >* libMesh::HPCoarsenTest::dphi
protected

Definition at line 126 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

const std::vector<std::vector<RealGradient> > * libMesh::HPCoarsenTest::dphi_coarse
protected

Definition at line 126 of file hp_coarsentest.h.

Referenced by select_refinement().

UniquePtr<FEBase> libMesh::HPCoarsenTest::fe
protected

The finite element objects for fine and coarse elements

Definition at line 120 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

DenseVector<Number> libMesh::HPCoarsenTest::Fe
protected

Definition at line 149 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

UniquePtr<FEBase> libMesh::HPCoarsenTest::fe_coarse
protected

Definition at line 120 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

const std::vector<Real>* libMesh::HPCoarsenTest::JxW
protected

Mapping jacobians

Definition at line 132 of file hp_coarsentest.h.

Referenced by select_refinement().

DenseMatrix<Number> libMesh::HPCoarsenTest::Ke
protected

Linear system for projections

Definition at line 148 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

Real libMesh::HPCoarsenTest::p_weight

Because the coarsening test seems to always choose p refinement, we're providing an option to make h refinement more likely

Definition at line 98 of file hp_coarsentest.h.

Referenced by select_refinement().

const std::vector<std::vector<Real> >* libMesh::HPCoarsenTest::phi
protected

The shape functions and their derivatives

Definition at line 125 of file hp_coarsentest.h.

Referenced by select_refinement().

const std::vector<std::vector<Real> > * libMesh::HPCoarsenTest::phi_coarse
protected

Definition at line 125 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

UniquePtr<QBase> libMesh::HPCoarsenTest::qrule
protected

The quadrature rule for the fine element

Definition at line 143 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

DenseVector<Number> libMesh::HPCoarsenTest::Uc
protected

Coefficients for projected coarse and projected p-derefined solutions

Definition at line 154 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

DenseVector<Number> libMesh::HPCoarsenTest::Up
protected

Definition at line 155 of file hp_coarsentest.h.

Referenced by select_refinement().

const std::vector<Point>* libMesh::HPCoarsenTest::xyz_values
protected

Quadrature locations

Definition at line 137 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().


The documentation for this class was generated from the following files: