libMesh
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
libMesh::HPCoarsenTest Class Reference

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. More...

#include <hp_coarsentest.h>

Inheritance diagram for libMesh::HPCoarsenTest:
[legend]

Public Member Functions

 HPCoarsenTest ()
 Constructor. More...
 
virtual ~HPCoarsenTest ()
 Destructor. More...
 
virtual void select_refinement (System &system) libmesh_override
 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. More...
 

Public Attributes

Real p_weight
 Because the coarsening test seems to always choose p refinement, we're providing an option to make h refinement more likely. More...
 
std::vector< float > component_scale
 This vector can be used to "scale" certain variables in a system. More...
 

Protected Member Functions

void add_projection (const System &, const Elem *, unsigned int var)
 The helper function which adds individual fine element data to the coarse element projection. More...
 

Protected Attributes

const Elemcoarse
 The coarse element on which a solution projection is cached. More...
 
std::vector< dof_id_typedof_indices
 Global DOF indices for fine elements. More...
 
UniquePtr< FEBasefe
 The finite element objects for fine and coarse elements. More...
 
UniquePtr< FEBasefe_coarse
 
const std::vector< std::vector< Real > > * phi
 The shape functions and their derivatives. More...
 
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
 Mapping jacobians. More...
 
const std::vector< Point > * xyz_values
 Quadrature locations. More...
 
std::vector< Pointcoarse_qpoints
 
UniquePtr< QBaseqrule
 The quadrature rule for the fine element. More...
 
DenseMatrix< NumberKe
 Linear system for projections. More...
 
DenseVector< NumberFe
 
DenseVector< NumberUc
 Coefficients for projected coarse and projected p-derefined solutions. More...
 
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 ( )

Constructor.

Definition at line 75 of file hp_coarsentest.h.

75  : p_weight(1.0)
76  {
77  libmesh_experimental();
78  }
Real p_weight
Because the coarsening test seems to always choose p refinement, we&#39;re providing an option to make h ...
virtual libMesh::HPCoarsenTest::~HPCoarsenTest ( )
virtual

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 47 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_ref_range(), 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(), 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().

50 {
51  // If we have children, we need to add their projections instead
52  if (!elem->active())
53  {
54  libmesh_assert(!elem->subactive());
55  for (auto & child : elem->child_ref_range())
56  this->add_projection(system, &child, var);
57  return;
58  }
59 
60  // The DofMap for this system
61  const DofMap & dof_map = system.get_dof_map();
62 
63  // The type of finite element to use for this variable
64  const FEType & fe_type = dof_map.variable_type (var);
65 
66  const FEContinuity cont = fe->get_continuity();
67 
68  fe->reinit(elem);
69 
70  dof_map.dof_indices(elem, dof_indices, var);
71 
72  const unsigned int n_dofs =
73  cast_int<unsigned int>(dof_indices.size());
74 
75  FEInterface::inverse_map (system.get_mesh().mesh_dimension(),
76  fe_type, coarse, *xyz_values, coarse_qpoints);
77 
78  fe_coarse->reinit(coarse, &coarse_qpoints);
79 
80  const unsigned int n_coarse_dofs =
81  cast_int<unsigned int>(phi_coarse->size());
82 
83  if (Uc.size() == 0)
84  {
85  Ke.resize(n_coarse_dofs, n_coarse_dofs);
86  Ke.zero();
87  Fe.resize(n_coarse_dofs);
88  Fe.zero();
89  Uc.resize(n_coarse_dofs);
90  Uc.zero();
91  }
92  libmesh_assert_equal_to (Uc.size(), phi_coarse->size());
93 
94  // Loop over the quadrature points
95  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
96  {
97  // The solution value at the quadrature point
98  Number val = libMesh::zero;
99  Gradient grad;
100  Tensor hess;
101 
102  for (unsigned int i=0; i != n_dofs; i++)
103  {
104  dof_id_type dof_num = dof_indices[i];
105  val += (*phi)[i][qp] *
106  system.current_solution(dof_num);
107  if (cont == C_ZERO || cont == C_ONE)
108  grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
109  // grad += (*dphi)[i][qp] *
110  // system.current_solution(dof_num);
111  if (cont == C_ONE)
112  hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
113  // hess += (*d2phi)[i][qp] *
114  // system.current_solution(dof_num);
115  }
116 
117  // The projection matrix and vector
118  for (unsigned int i=0; i != Fe.size(); ++i)
119  {
120  Fe(i) += (*JxW)[qp] *
121  (*phi_coarse)[i][qp]*val;
122  if (cont == C_ZERO || cont == C_ONE)
123  Fe(i) += (*JxW)[qp] *
124  (grad*(*dphi_coarse)[i][qp]);
125  if (cont == C_ONE)
126  Fe(i) += (*JxW)[qp] *
127  hess.contract((*d2phi_coarse)[i][qp]);
128  // Fe(i) += (*JxW)[qp] *
129  // (*d2phi_coarse)[i][qp].contract(hess);
130 
131  for (unsigned int j=0; j != Fe.size(); ++j)
132  {
133  Ke(i,j) += (*JxW)[qp] *
134  (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
135  if (cont == C_ZERO || cont == C_ONE)
136  Ke(i,j) += (*JxW)[qp] *
137  (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
138  if (cont == C_ONE)
139  Ke(i,j) += (*JxW)[qp] *
140  ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
141  }
142  }
143  }
144 }
void add_projection(const System &, const Elem *, unsigned int var)
The helper function which adds individual fine element data to the coarse element projection...
virtual void zero() libmesh_override
Set every element in the matrix to 0.
Definition: dense_matrix.h:792
DenseVector< Number > Fe
const std::vector< std::vector< Real > > * phi_coarse
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
UniquePtr< QBase > qrule
The quadrature rule for the fine element.
UniquePtr< FEBase > fe
The finite element objects for fine and coarse elements.
const Number zero
.
Definition: libmesh.h:178
const std::vector< Point > * xyz_values
Quadrature locations.
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
DenseVector< Number > Uc
Coefficients for projected coarse and projected p-derefined solutions.
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:569
UniquePtr< FEBase > fe_coarse
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
NumberTensorValue Tensor
const Elem * coarse
The coarse element on which a solution projection is cached.
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
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
Global DOF indices for fine elements.
std::vector< Point > coarse_qpoints
uint8_t dof_id_type
Definition: id_types.h:64
DenseMatrix< Number > Ke
Linear system for projections.
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 146 of file hp_coarsentest.C.

References libMesh::MeshBase::active_local_element_ptr_range(), 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(), dim, 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::FEInterface::inverse_map(), JxW, Ke, libMesh::libmesh_assert(), libmesh_nullptr, libMesh::MeshRefinement::make_flags_parallel_consistent(), std::max(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::FEInterface::n_dofs(), n_nodes, n_vars, libMesh::System::n_vars(), 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::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseVector< T >::size(), libMesh::TypeVector< T >::subtract_scaled(), libMesh::TypeTensor< T >::subtract_scaled(), Uc, Up, libMesh::DofMap::variable_type(), xyz_values, libMesh::DenseMatrix< T >::zero(), libMesh::DenseVector< T >::zero(), and libMesh::zero.

Referenced by main(), and ~HPCoarsenTest().

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

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: