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...
 
 HPCoarsenTest (const HPCoarsenTest &)=delete
 This class cannot be (default) copy constructed/assigned because it has unique_ptr members. More...
 
HPCoarsenTestoperator= (const HPCoarsenTest &)=delete
 
 HPCoarsenTest (HPCoarsenTest &&)=default
 Defaulted move ctor, move assignment operator, and destructor. More...
 
HPCoarsenTestoperator= (HPCoarsenTest &&)=default
 
virtual ~HPCoarsenTest ()=default
 
virtual void select_refinement (System &system) 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...
 
void extra_quadrature_order (const int extraorder)
 Increases or decreases the order of the quadrature rule used for numerical integration. 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

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...
 
std::unique_ptr< FEBasefe
 The finite element objects for fine and coarse elements. More...
 
std::unique_ptr< 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
 
std::unique_ptr< 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
 
int _extra_order
 Extra order to use for quadrature rule. More...
 

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 67 of file hp_coarsentest.h.

Constructor & Destructor Documentation

◆ HPCoarsenTest() [1/3]

libMesh::HPCoarsenTest::HPCoarsenTest ( )
inline

Constructor.

Definition at line 74 of file hp_coarsentest.h.

74  : p_weight(1.0), _extra_order(1)
75  {
76  libmesh_experimental();
77  }
int _extra_order
Extra order to use for quadrature rule.
Real p_weight
Because the coarsening test seems to always choose p refinement, we&#39;re providing an option to make h ...

◆ HPCoarsenTest() [2/3]

libMesh::HPCoarsenTest::HPCoarsenTest ( const HPCoarsenTest )
delete

This class cannot be (default) copy constructed/assigned because it has unique_ptr members.

Explicitly deleting these functions is the best way to document this fact.

◆ HPCoarsenTest() [3/3]

libMesh::HPCoarsenTest::HPCoarsenTest ( HPCoarsenTest &&  )
default

Defaulted move ctor, move assignment operator, and destructor.

◆ ~HPCoarsenTest()

virtual libMesh::HPCoarsenTest::~HPCoarsenTest ( )
virtualdefault

Member Function Documentation

◆ add_projection()

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::TypeTensor< T >::add_scaled(), libMesh::TypeVector< 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::index_range(), libMesh::FEMap::inverse_map(), Ke, libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::mesh_dimension(), phi_coarse, qrule, libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseVector< T >::size(), libMesh::Elem::subactive(), Uc, 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  const FEContinuity cont = fe->get_continuity();
64 
65  fe->reinit(elem);
66 
67  dof_map.dof_indices(elem, dof_indices, var);
68 
69  const unsigned int n_dofs =
70  cast_int<unsigned int>(dof_indices.size());
71 
72  FEMap::inverse_map (system.get_mesh().mesh_dimension(), coarse,
74 
75  fe_coarse->reinit(coarse, &coarse_qpoints);
76 
77  const unsigned int n_coarse_dofs =
78  cast_int<unsigned int>(phi_coarse->size());
79 
80  if (Uc.size() == 0)
81  {
82  Ke.resize(n_coarse_dofs, n_coarse_dofs);
83  Ke.zero();
84  Fe.resize(n_coarse_dofs);
85  Fe.zero();
86  Uc.resize(n_coarse_dofs);
87  Uc.zero();
88  }
89  libmesh_assert_equal_to (Uc.size(), phi_coarse->size());
90 
91  // Loop over the quadrature points
92  for (auto qp : make_range(qrule->n_points()))
93  {
94  // The solution value at the quadrature point
95  Number val = libMesh::zero;
96  Gradient grad;
97  Tensor hess;
98 
99  for (unsigned int i=0; i != n_dofs; i++)
100  {
101  dof_id_type dof_num = dof_indices[i];
102  val += (*phi)[i][qp] *
103  system.current_solution(dof_num);
104  if (cont == C_ZERO || cont == C_ONE)
105  grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
106  if (cont == C_ONE)
107  hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
108  }
109 
110  // The projection matrix and vector
111  for (auto i : index_range(Fe))
112  {
113  Fe(i) += (*JxW)[qp] *
114  (*phi_coarse)[i][qp]*val;
115  if (cont == C_ZERO || cont == C_ONE)
116  Fe(i) += (*JxW)[qp] *
117  (grad*(*dphi_coarse)[i][qp]);
118  if (cont == C_ONE)
119  Fe(i) += (*JxW)[qp] *
120  hess.contract((*d2phi_coarse)[i][qp]);
121 
122  for (auto j : index_range(Fe))
123  {
124  Ke(i,j) += (*JxW)[qp] *
125  (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
126  if (cont == C_ZERO || cont == C_ONE)
127  Ke(i,j) += (*JxW)[qp] *
128  (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
129  if (cont == C_ONE)
130  Ke(i,j) += (*JxW)[qp] *
131  ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
132  }
133  }
134  }
135 }
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() override final
Set every element in the vector to 0.
Definition: dense_vector.h:398
DenseVector< Number > Fe
virtual void zero() override final
Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously...
Definition: dense_matrix.h:911
static Point inverse_map(const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true, const bool extra_checks=true)
Definition: fe_map.C:1626
std::unique_ptr< FEBase > fe_coarse
const std::vector< std::vector< Real > > * phi_coarse
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
const Number zero
.
Definition: libmesh.h:280
const std::vector< Point > * xyz_values
Quadrature locations.
DenseVector< Number > Uc
Coefficients for projected coarse and projected p-derefined solutions.
NumberVectorValue Gradient
libmesh_assert(ctx)
std::unique_ptr< FEBase > fe
The finite element objects for fine and coarse elements.
NumberTensorValue Tensor
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:134
virtual unsigned int size() const override final
Definition: dense_vector.h:104
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
Elem * coarse
The coarse element on which a solution projection is cached.
std::unique_ptr< QBase > qrule
The quadrature rule for the fine element.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:111
uint8_t dof_id_type
Definition: id_types.h:67
DenseMatrix< Number > Ke
Linear system for projections.

◆ extra_quadrature_order()

void libMesh::HPCoarsenTest::extra_quadrature_order ( const int  extraorder)
inline

Increases or decreases the order of the quadrature rule used for numerical integration.

The default extraorder is 1, because properly integrating L2 error requires integrating the squares of terms with order p+1, and 2p+2 is 1 higher than what we default to using for reasonable mass matrix integration.

Definition at line 115 of file hp_coarsentest.h.

References _extra_order.

116  { _extra_order = extraorder; }
int _extra_order
Extra order to use for quadrature rule.

◆ operator=() [1/2]

HPCoarsenTest& libMesh::HPCoarsenTest::operator= ( const HPCoarsenTest )
delete

◆ operator=() [2/2]

HPCoarsenTest& libMesh::HPCoarsenTest::operator= ( HPCoarsenTest &&  )
default

◆ select_refinement()

void libMesh::HPCoarsenTest::select_refinement ( System system)
overridevirtual

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 137 of file hp_coarsentest.C.

References _extra_order, add_projection(), libMesh::TypeTensor< T >::add_scaled(), libMesh::TypeVector< 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::Elem::hack_p_level(), libMesh::index_range(), libMesh::FEMap::inverse_map(), JxW, Ke, libMesh::libmesh_assert(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::make_range(), mesh, 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::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(), std::sqrt(), libMesh::TypeTensor< T >::subtract_scaled(), libMesh::TypeVector< 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().

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

Member Data Documentation

◆ _extra_order

int libMesh::HPCoarsenTest::_extra_order
protected

Extra order to use for quadrature rule.

Definition at line 178 of file hp_coarsentest.h.

Referenced by extra_quadrature_order(), and select_refinement().

◆ coarse

Elem* libMesh::HPCoarsenTest::coarse
protected

The coarse element on which a solution projection is cached.

Definition at line 128 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ coarse_qpoints

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

Definition at line 156 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ component_scale

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 82 of file hp_selector.h.

Referenced by select_refinement().

◆ d2phi

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

Definition at line 145 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ d2phi_coarse

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

Definition at line 145 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ dof_indices

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

Global DOF indices for fine elements.

Definition at line 133 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ dphi

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

Definition at line 144 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ dphi_coarse

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

Definition at line 144 of file hp_coarsentest.h.

Referenced by select_refinement().

◆ fe

std::unique_ptr<FEBase> libMesh::HPCoarsenTest::fe
protected

The finite element objects for fine and coarse elements.

Definition at line 138 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ Fe

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

Definition at line 167 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ fe_coarse

std::unique_ptr<FEBase> libMesh::HPCoarsenTest::fe_coarse
protected

Definition at line 138 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ JxW

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

Mapping jacobians.

Definition at line 150 of file hp_coarsentest.h.

Referenced by select_refinement().

◆ Ke

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

Linear system for projections.

Definition at line 166 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ p_weight

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 106 of file hp_coarsentest.h.

Referenced by select_refinement().

◆ phi

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

The shape functions and their derivatives.

Definition at line 143 of file hp_coarsentest.h.

Referenced by select_refinement().

◆ phi_coarse

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

Definition at line 143 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ qrule

std::unique_ptr<QBase> libMesh::HPCoarsenTest::qrule
protected

The quadrature rule for the fine element.

Definition at line 161 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ Uc

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

Coefficients for projected coarse and projected p-derefined solutions.

Definition at line 172 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().

◆ Up

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

Definition at line 173 of file hp_coarsentest.h.

Referenced by select_refinement().

◆ xyz_values

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

Quadrature locations.

Definition at line 155 of file hp_coarsentest.h.

Referenced by add_projection(), and select_refinement().


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