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

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.
 
std::vector< float > component_scale
 This vector can be used to "scale" certain variables in a system.
 

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.
 

Protected Attributes

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

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

51{
52 // If we have children, we need to add their projections instead
53 if (!elem->active())
54 {
55 libmesh_assert(!elem->subactive());
56 for (auto & child : elem->child_ref_range())
57 this->add_projection(system, &child, var);
58 return;
59 }
60
61 // The DofMap for this system
62 const DofMap & dof_map = system.get_dof_map();
63
64 const FEContinuity cont = fe->get_continuity();
65
66 fe->reinit(elem);
67
68 dof_map.dof_indices(elem, dof_indices, var);
69
70 const unsigned int n_dofs =
71 cast_int<unsigned int>(dof_indices.size());
72
73 FEMap::inverse_map (system.get_mesh().mesh_dimension(), coarse,
75
77
78 const unsigned int n_coarse_dofs =
79 cast_int<unsigned int>(phi_coarse->size());
80
81 if (Uc.size() == 0)
82 {
83 Ke.resize(n_coarse_dofs, n_coarse_dofs);
84 Ke.zero();
85 Fe.resize(n_coarse_dofs);
86 Fe.zero();
87 Uc.resize(n_coarse_dofs);
88 Uc.zero();
89 }
90 libmesh_assert_equal_to (Uc.size(), phi_coarse->size());
91
92 // Loop over the quadrature points
93 for (auto qp : make_range(qrule->n_points()))
94 {
95 // The solution value at the quadrature point
97 Gradient grad;
98 Tensor hess;
99
100 for (unsigned int i=0; i != n_dofs; i++)
101 {
102 dof_id_type dof_num = dof_indices[i];
103 val += (*phi)[i][qp] *
104 system.current_solution(dof_num);
105 if (cont == C_ZERO || cont == C_ONE)
106 grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
107 if (cont == C_ONE)
108 hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
109 }
110
111 // The projection matrix and vector
112 for (auto i : index_range(Fe))
113 {
114 Fe(i) += (*JxW)[qp] *
115 (*phi_coarse)[i][qp]*val;
116 if (cont == C_ZERO || cont == C_ONE)
117 Fe(i) += (*JxW)[qp] *
118 (grad*(*dphi_coarse)[i][qp]);
119 if (cont == C_ONE)
120 Fe(i) += (*JxW)[qp] *
121 hess.contract((*d2phi_coarse)[i][qp]);
122
123 for (auto j : index_range(Fe))
124 {
125 Ke(i,j) += (*JxW)[qp] *
126 (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
127 if (cont == C_ZERO || cont == C_ONE)
128 Ke(i,j) += (*JxW)[qp] *
129 (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
130 if (cont == C_ONE)
131 Ke(i,j) += (*JxW)[qp] *
132 ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
133 }
134 }
135 }
136}
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
virtual void zero() override final
Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously...
void resize(const unsigned int n)
Resize the vector.
virtual void zero() override final
Set every element in the vector to 0.
virtual unsigned int size() const override final
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:1512
DenseVector< Number > Fe
std::vector< dof_id_type > dof_indices
Global DOF indices for fine elements.
DenseMatrix< Number > Ke
Linear system for projections.
DenseVector< Number > Uc
Coefficients for projected coarse and projected p-derefined solutions.
const std::vector< Point > * xyz_values
Quadrature locations.
std::unique_ptr< FEBase > fe
The finite element objects for fine and coarse elements.
std::unique_ptr< FEBase > fe_coarse
const std::vector< std::vector< RealTensor > > * d2phi_coarse
const std::vector< std::vector< Real > > * phi_coarse
void add_projection(const System &, const Elem *, unsigned int var)
The helper function which adds individual fine element data to the coarse element projection.
Elem * coarse
The coarse element on which a solution projection is cached.
std::unique_ptr< QBase > qrule
The quadrature rule for the fine element.
const std::vector< std::vector< RealTensor > > * d2phi
std::vector< Point > coarse_qpoints
const std::vector< std::vector< RealGradient > > * dphi
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
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:153
NumberVectorValue Gradient
libmesh_assert(ctx)
NumberTensorValue Tensor
const Number zero
.
Definition libmesh.h:297
uint8_t dof_id_type
Definition id_types.h:67
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:176

References libMesh::Elem::active(), add_projection(), 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, libMesh::DofMap::dof_indices(), 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::DofMap::reinit(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), libMesh::DenseVector< T >::size(), libMesh::Elem::subactive(), Uc, xyz_values, libMesh::zero, libMesh::DenseMatrix< T >::zero(), and libMesh::DenseVector< T >::zero().

Referenced by add_projection(), and select_refinement().

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

116 { _extra_order = extraorder; }

References _extra_order.

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

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

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, libMesh::DofMap::dof_indices(), dof_indices, libMesh::DofObject::dof_number(), dphi, dphi_coarse, libMesh::ErrorVectorReal, libMesh::SmoothnessEstimator::estimate_smoothness(), 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(), libMesh::ErrorVector::mean(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::FEInterface::n_dofs(), n_nodes, n_vars, libMesh::System::n_vars(), libMesh::TypeTensor< T >::norm_sq(), libMesh::TypeVector< T >::norm_sq(), libMesh::TensorTools::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(), libMesh::TypeTensor< T >::subtract_scaled(), libMesh::TypeVector< T >::subtract_scaled(), Uc, Up, libMesh::DofMap::variable_type(), xyz_values, libMesh::zero, libMesh::DenseMatrix< T >::zero(), and libMesh::DenseVector< T >::zero().

Referenced by main().

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: