libMesh
hp_coarsentest.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // C++ includes
19 #include <limits> // for std::numeric_limits::max
20 #include <math.h> // for sqrt
21 
22 
23 // Local Includes
24 #include "libmesh/hp_coarsentest.h"
25 #include "libmesh/dense_matrix.h"
26 #include "libmesh/dense_vector.h"
27 #include "libmesh/dof_map.h"
28 #include "libmesh/fe_base.h"
29 #include "libmesh/fe_interface.h"
30 #include "libmesh/libmesh_logging.h"
31 #include "libmesh/elem.h"
32 #include "libmesh/error_vector.h"
33 #include "libmesh/mesh_base.h"
34 #include "libmesh/mesh_refinement.h"
35 #include "libmesh/quadrature.h"
36 #include "libmesh/system.h"
37 #include "libmesh/tensor_value.h"
38 #include "libmesh/smoothness_estimator.h"
39 
40 #ifdef LIBMESH_ENABLE_AMR
41 
42 namespace libMesh
43 {
44 
45 //-----------------------------------------------------------------
46 // HPCoarsenTest implementations
47 
49  const Elem * elem,
50  unsigned int var)
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 
75 
76  fe_coarse->reinit(coarse, &coarse_qpoints);
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
96  Number val = libMesh::zero;
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 }
137 
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
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 
438  fe_coarse->reinit(coarse, &coarse_qpoints);
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.
557 }
558 
559 } // namespace libMesh
560 
561 #endif // #ifdef LIBMESH_ENABLE_AMR
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
const std::vector< std::vector< Real > > * phi
The shape functions and their derivatives.
const Elem * parent() const
Definition: elem.h:3030
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
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...
A Node is like a Point, but with more information.
Definition: node.h:52
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:355
virtual void zero() override final
Set every element in the vector to 0.
Definition: dense_vector.h:420
void subtract_scaled(const TypeVector< T2 > &, const T &)
Subtract a scaled value from this vector without creating a temporary.
Definition: type_vector.h:703
DenseVector< Number > Fe
virtual void select_refinement(System &system) override
This pure virtual function must be redefined in derived classes to take a mesh flagged for h refineme...
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2229
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:1628
The ErrorVector is a specialization of the StatisticsVector for error data computed on a finite eleme...
Definition: error_vector.h:50
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:396
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2230
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
auto norm_sq() const -> decltype(std::norm(T()))
Definition: type_tensor.h:1348
MeshBase & mesh
const std::vector< std::vector< RealGradient > > * dphi_coarse
unsigned int p_level() const
Definition: elem.h:3108
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:304
const std::vector< Point > * xyz_values
Quadrature locations.
const MeshBase & get_mesh() const
Definition: system.h:2358
DIE A HORRIBLE DEATH HERE typedef float ErrorVectorReal
DenseVector< Number > Uc
Coefficients for projected coarse and projected p-derefined solutions.
This is the MeshBase class.
Definition: mesh_base.h:75
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition: elem.h:2346
Number current_solution(const dof_id_type global_dof_number) const
Definition: system.C:165
const std::vector< Real > * JxW
Mapping jacobians.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:34
const dof_id_type n_nodes
Definition: tecplot_io.C:67
unsigned int number() const
Definition: system.h:2350
auto norm_sq() const -> decltype(std::norm(T()))
Definition: type_vector.h:926
void subtract_scaled(const TypeTensor< T2 > &, const T &)
Subtract a scaled value from this tensor without creating a temporary.
Definition: type_tensor.h:921
Implements (adaptive) mesh refinement algorithms for a MeshBase.
unsigned int n_vars
This class implements the Smoothness estimate.
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libmesh_assert(ctx)
bool make_flags_parallel_consistent()
Copy refinement flags on ghost elements from their local processors.
virtual void estimate_smoothness(const System &system, ErrorVector &smoothness_per_cell, const NumericVector< Number > *solution_vector=nullptr)
This function uses the Legendre expansion of solution to estimate coefficient decay to quantify the s...
void add_scaled(const TypeTensor< T2 > &, const T &)
Add a scaled tensor to this tensor without creating a temporary.
Definition: type_tensor.h:851
DenseVector< Number > Up
CompareTypes< T, T2 >::supertype contract(const TypeTensor< T2 > &) const
Multiply 2 tensors together to return a scalar, i.e.
Definition: type_tensor.h:1268
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
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
bool subactive() const
Definition: elem.h:2959
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:140
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
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
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
const std::vector< std::vector< RealGradient > > * dphi
unsigned int n_vars() const
Definition: system.h:2430
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
bool active() const
Definition: elem.h:2941
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:3264
const DofMap & get_dof_map() const
Definition: system.h:2374
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:117
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.
uint8_t dof_id_type
Definition: id_types.h:67
virtual Real mean() const override
Definition: error_vector.C:69
DenseMatrix< Number > Ke
Linear system for projections.