libMesh
Public Member Functions | Static Public Member Functions | Private Attributes | List of all members
Biharmonic::JR Class Reference

Biharmonic's friend class definition. More...

#include <biharmonic_jr.h>

Inheritance diagram for Biharmonic::JR:
[legend]

Public Member Functions

 JR (EquationSystems &eqSys, const std::string &name, const unsigned int number)
 Constructor. More...
 
void initialize () override
 
void residual (const NumericVector< Number > &u, NumericVector< Number > &R, NonlinearImplicitSystem &sys) override
 The residual assembly function for the Biharmonic system. More...
 
void residual_and_jacobian (const NumericVector< Number > &u, NumericVector< Number > *R, SparseMatrix< Number > *J, NonlinearImplicitSystem &) override
 The residual and Jacobian assembly function for the Biharmonic system. More...
 
void bounds (NumericVector< Number > &XL, NumericVector< Number > &XU, NonlinearImplicitSystem &) override
 Function defining the bounds of the Biharmonic system. More...
 

Static Public Member Functions

static Number InitialDensityBall (const Point &p, const Parameters &parameters, const std::string &, const std::string &)
 Static functions to be used for initialization. More...
 
static Number InitialDensityRod (const Point &p, const Parameters &parameters, const std::string &, const std::string &)
 
static Number InitialDensityStrip (const Point &p, const Parameters &parameters, const std::string &, const std::string &)
 
static Gradient InitialGradientZero (const Point &, const Parameters &, const std::string &, const std::string &)
 

Private Attributes

Biharmonic_biharmonic
 

Detailed Description

Biharmonic's friend class definition.

Definition at line 29 of file biharmonic_jr.h.

Constructor & Destructor Documentation

◆ JR()

Biharmonic::JR::JR ( EquationSystems eqSys,
const std::string &  name,
const unsigned int  number 
)

Constructor.

Definition at line 19 of file biharmonic_jr.C.

References _biharmonic, Biharmonic::_dim, libMesh::DofMap::add_periodic_boundary(), libMesh::HERMITE, libMesh::PeriodicBoundaryBase::myboundary, libMesh::PeriodicBoundaryBase::pairedboundary, and libMesh::THIRD.

21  :
22  TransientNonlinearImplicitSystem(eqSys, name_in, number_in),
23  _biharmonic(dynamic_cast<Biharmonic &>(eqSys))
24 {
25  // Check that we can actually compute second derivatives
26 #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
27  libmesh_error_msg("Must have second derivatives enabled");
28 #endif
29 
30 #ifdef LIBMESH_ENABLE_PERIODIC
31  // Add periodicity to the mesh
32  DofMap & dof_map = get_dof_map();
33  PeriodicBoundary xbdry(RealVectorValue(1.0, 0.0, 0.0));
34 #if LIBMESH_DIM > 1
35  PeriodicBoundary ybdry(RealVectorValue(0.0, 1.0, 0.0));
36 #endif
37 #if LIBMESH_DIM > 2
38  PeriodicBoundary zbdry(RealVectorValue(0.0, 0.0, 1.0));
39 #endif
40 
41  switch(_biharmonic._dim)
42  {
43  case 1:
44  xbdry.myboundary = 0;
45  xbdry.pairedboundary = 1;
46  dof_map.add_periodic_boundary(xbdry);
47  break;
48 #if LIBMESH_DIM > 1
49  case 2:
50  xbdry.myboundary = 3;
51  xbdry.pairedboundary = 1;
52  dof_map.add_periodic_boundary(xbdry);
53  ybdry.myboundary = 0;
54  ybdry.pairedboundary = 2;
55  dof_map.add_periodic_boundary(ybdry);
56  break;
57 #endif
58 #if LIBMESH_DIM > 2
59  case 3:
60  xbdry.myboundary = 4;
61  xbdry.pairedboundary = 2;
62  dof_map.add_periodic_boundary(xbdry);
63  ybdry.myboundary = 1;
64  ybdry.pairedboundary = 3;
65  dof_map.add_periodic_boundary(ybdry);
66  zbdry.myboundary = 0;
67  zbdry.pairedboundary = 5;
68  dof_map.add_periodic_boundary(zbdry);
69  break;
70 #endif
71  default:
72  libmesh_error_msg("Invalid dimension = " << _biharmonic._dim);
73  }
74 #endif // LIBMESH_ENABLE_PERIODIC
75 
76  // Adds the variable "u" to the system.
77  // u will be approximated using Hermite elements
78  add_variable("u", THIRD, HERMITE);
79 
80  // Give the system an object to compute the initial state.
81  attach_init_object(*this);
82 
83  // Attach the R & J calculation object
84  nonlinear_solver->residual_and_jacobian_object = this;
85 
86  // And the R object for matrix-free evaluations
87  nonlinear_solver->mffd_residual_object = this;
88 
89  // Attach the bounds calculation object
90  nonlinear_solver->bounds_object = this;
91 }
Biharmonic & _biharmonic
Definition: biharmonic_jr.h:91
void add_periodic_boundary(const PeriodicBoundaryBase &periodic_boundary)
Adds a copy of the specified periodic boundary to the system.
VectorValue< Real > RealVectorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
TransientSystem< NonlinearImplicitSystem > TransientNonlinearImplicitSystem
The definition of a periodic boundary.
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
unsigned int _dim
Definition: biharmonic.h:92

Member Function Documentation

◆ bounds()

void Biharmonic::JR::bounds ( NumericVector< Number > &  XL,
NumericVector< Number > &  XU,
NonlinearImplicitSystem  
)
override

Function defining the bounds of the Biharmonic system.

Definition at line 498 of file biharmonic_jr.C.

References std::abs(), libMesh::DofMap::dof_indices(), libMesh::NumericVector< T >::insert(), libMesh::EquationSystems::n_dofs(), libMesh::Real, libMesh::DenseVector< T >::resize(), and libMesh::DofMap::variable_type().

501 {
502  // sys is actually ignored, since it should be the same as *this.
503 
504  // Declare a performance log. Give it a descriptive
505  // string to identify what part of the code we are
506  // logging, since there may be many PerfLogs in an
507  // application.
508  PerfLog perf_log ("Biharmonic bounds", false);
509 
510  // A reference to the DofMap object for this system. The DofMap
511  // object handles the index translation from node and element numbers
512  // to degree of freedom numbers. We will talk more about the DofMap
513  // in future examples.
514  const DofMap & dof_map = get_dof_map();
515 
516  // Get a constant reference to the Finite Element type
517  // for the first (and only) variable in the system.
518  FEType fe_type = dof_map.variable_type(0);
519 
520  // Build a Finite Element object of the specified type. Since the
521  // FEBase::build() member dynamically creates memory we will
522  // store the object as a std::unique_ptr<FEBase>. This can be thought
523  // of as a pointer that will clean up after itself.
524  std::unique_ptr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
525 
526  // Define data structures to contain the bound vectors contributions.
527  DenseVector<Number> XLe, XUe;
528 
529  // These vector will hold the degree of freedom indices for
530  // the element. These define where in the global system
531  // the element degrees of freedom get mapped.
532  std::vector<dof_id_type> dof_indices;
533 
534  for (const auto & elem : _biharmonic._mesh.active_local_element_ptr_range())
535  {
536  // Extract the shape function to be evaluated at the nodes
537  const std::vector<std::vector<Real>> & phi = fe->get_phi();
538 
539  // Get the degree of freedom indices for the current element.
540  // They are in 1-1 correspondence with shape functions phi
541  // and define where in the global vector this element will.
542  dof_map.dof_indices (elem, dof_indices);
543 
544  const unsigned int n_dofs =
545  cast_int<unsigned int>(dof_indices.size());
546 
547  // Resize the local bounds vectors (zeroing them out in the process).
548  XLe.resize(n_dofs);
549  XUe.resize(n_dofs);
550 
551  // Extract the element node coordinates in the reference frame
552  std::vector<Point> nodes;
553  fe->get_refspace_nodes(elem->type(), nodes);
554 
555  // Evaluate the shape functions at the nodes
556  fe->reinit(elem, &nodes);
557 
558  // Construct the bounds based on the value of the i-th phi at the nodes.
559  // Observe that this doesn't really work in general: we rely on the fact
560  // that for Hermite elements each shape function is nonzero at most at a
561  // single node.
562  // More generally the bounds must be constructed by inspecting a "mass-like"
563  // matrix (m_{ij}) of the shape functions (i) evaluated at their corresponding nodes (j).
564  // The constraints imposed on the dofs (d_i) are then are -1 \leq \sum_i d_i m_{ij} \leq 1,
565  // since \sum_i d_i m_{ij} is the value of the solution at the j-th node.
566  // Auxiliary variables will need to be introduced to reduce this to a "box" constraint.
567  // Additional complications will arise since m might be singular (as is the case for Hermite,
568  // which, however, is easily handled by inspection).
569  for (unsigned int i=0; i<n_dofs; ++i)
570  {
571  // FIXME: should be able to define INF and pass it to the solve
572  Real infinity = 1.0e20;
573  Real bound = infinity;
574  for (unsigned int j = 0; j < nodes.size(); ++j)
575  {
576  if (phi[i][j])
577  {
578  bound = 1.0/std::abs(phi[i][j]);
579  break;
580  }
581  }
582 
583  // The value of the solution at this node must be between 1.0 and -1.0.
584  // Based on the value of phi(i)(i) the nodal coordinate must be between 1.0/phi(i)(i) and its negative.
585  XLe(i) = -bound;
586  XUe(i) = bound;
587  }
588  // The element bound vectors are now built for this element.
589  // Insert them into the global vectors, potentially overwriting
590  // the same dof contributions from other elements: no matter --
591  // the bounds are always -1.0 and 1.0.
592  XL.insert(XLe, dof_indices);
593  XU.insert(XUe, dof_indices);
594  }
595 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
Inserts the entries of v in *this at the locations specified by v.
Biharmonic & _biharmonic
Definition: biharmonic_jr.h:91
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1992
std::size_t n_dofs() const
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2144
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:57
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:145
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
unsigned int _dim
Definition: biharmonic.h:92
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
ReplicatedMesh & _mesh
Definition: biharmonic.h:111

◆ InitialDensityBall()

Number Biharmonic::JR::InitialDensityBall ( const Point p,
const Parameters parameters,
const std::string &  ,
const std::string &   
)
static

Static functions to be used for initialization.

Definition at line 127 of file biharmonic_jr.C.

References libMesh::Parameters::get(), libMesh::EquationSystems::parameters, and libMesh::Real.

Referenced by initialize().

131 {
132  // Initialize with a ball in the middle, which is a segment in 1D, a disk in 2D and a ball in 3D.
133  Point center = parameters.get<Point>("center");
134  Real width = parameters.get<Real>("width");
135  Point pc = p-center;
136  Real r = pc.norm();
137  return (r < width) ? 1.0 : -0.5;
138 }
const T & get(std::string_view) const
Definition: parameters.h:426
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Parameters parameters
Data structure holding arbitrary parameters.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ InitialDensityRod()

Number Biharmonic::JR::InitialDensityRod ( const Point p,
const Parameters parameters,
const std::string &  ,
const std::string &   
)
static

Definition at line 143 of file biharmonic_jr.C.

References libMesh::Parameters::get(), libMesh::EquationSystems::parameters, and libMesh::Real.

Referenced by initialize().

147 {
148  // Initialize with a rod in the middle so that we have a z-homogeneous system to model the 2D disk.
149  Point center = parameters.get<Point>("center");
150  Real width = parameters.get<Real>("width");
151  Point pc = p-center;
152 #if LIBMESH_DIM > 2
153  pc(2) = 0;
154 #endif
155  Real r = pc.norm();
156  return (r < width) ? 1.0 : -0.5;
157 }
const T & get(std::string_view) const
Definition: parameters.h:426
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Parameters parameters
Data structure holding arbitrary parameters.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ InitialDensityStrip()

Number Biharmonic::JR::InitialDensityStrip ( const Point p,
const Parameters parameters,
const std::string &  ,
const std::string &   
)
static

Definition at line 163 of file biharmonic_jr.C.

References libMesh::Parameters::get(), libMesh::EquationSystems::parameters, libMesh::Real, and std::sqrt().

Referenced by initialize().

167 {
168  // Initialize with a wide strip in the middle so that we have a yz-homogeneous system to model the 1D.
169  Point center = parameters.get<Point>("center");
170  Real width = parameters.get<Real>("width");
171  Real r = sqrt((p(0)-center(0))*(p(0)-center(0)));
172  return (r < width) ? 1.0 : -0.5;
173 }
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
Definition: type_vector.h:53
const T & get(std::string_view) const
Definition: parameters.h:426
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Parameters parameters
Data structure holding arbitrary parameters.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ InitialGradientZero()

Gradient Biharmonic::JR::InitialGradientZero ( const Point ,
const Parameters ,
const std::string &  ,
const std::string &   
)
static

Definition at line 178 of file biharmonic_jr.C.

Referenced by initialize().

182 {
183  return Gradient(0.0, 0.0, 0.0);
184 }
NumberVectorValue Gradient

◆ initialize()

void Biharmonic::JR::initialize ( )
override

Definition at line 97 of file biharmonic_jr.C.

References Biharmonic::BALL, InitialDensityBall(), InitialDensityRod(), InitialDensityStrip(), InitialGradientZero(), libMesh::out, libMesh::EquationSystems::parameters, libMesh::Real, Biharmonic::ROD, libMesh::Parameters::set(), and Biharmonic::STRIP.

98 {
100  libMesh::out << ">>> Initializing Biharmonic::JR\n";
101 
105 
108 
111 
114 
115  // both states are equal
116  *(old_local_solution) = *(current_local_solution);
117 
118  if (_biharmonic._verbose)
119  libMesh::out << "<<< Initializing Biharmonic::JR\n";
120 }
Biharmonic & _biharmonic
Definition: biharmonic_jr.h:91
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:67
bool _verbose
Definition: biharmonic.h:98
Point _initialCenter
Definition: biharmonic.h:100
static Gradient InitialGradientZero(const Point &, const Parameters &, const std::string &, const std::string &)
static Number InitialDensityRod(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
static Number InitialDensityStrip(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
Real _initialWidth
Definition: biharmonic.h:101
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static Number InitialDensityBall(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
Static functions to be used for initialization.
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
Parameters parameters
Data structure holding arbitrary parameters.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
InitialStateEnum _initialState
Definition: biharmonic.h:99

◆ residual()

void Biharmonic::JR::residual ( const NumericVector< Number > &  u,
NumericVector< Number > &  R,
NonlinearImplicitSystem sys 
)
override

The residual assembly function for the Biharmonic system.

Definition at line 186 of file biharmonic_jr.C.

189 {
190  residual_and_jacobian(u, &R, nullptr, sys);
191 }
void residual_and_jacobian(const NumericVector< Number > &u, NumericVector< Number > *R, SparseMatrix< Number > *J, NonlinearImplicitSystem &) override
The residual and Jacobian assembly function for the Biharmonic system.

◆ residual_and_jacobian()

void Biharmonic::JR::residual_and_jacobian ( const NumericVector< Number > &  u,
NumericVector< Number > *  R,
SparseMatrix< Number > *  J,
NonlinearImplicitSystem  
)
override

The residual and Jacobian assembly function for the Biharmonic system.

Definition at line 193 of file biharmonic_jr.C.

References libMesh::SparseMatrix< T >::add_matrix(), libMesh::NumericVector< T >::add_vector(), libMesh::DofMap::constrain_element_matrix(), libMesh::DofMap::constrain_element_vector(), libMesh::FEType::default_quadrature_rule(), libMesh::DofMap::dof_indices(), Biharmonic::DOUBLE_OBSTACLE, Biharmonic::DOUBLE_WELL, libMesh::libmesh_ignore(), Biharmonic::LOG_DOUBLE_OBSTACLE, Biharmonic::LOG_DOUBLE_WELL, libMesh::EquationSystems::n_dofs(), libMesh::DenseVector< T >::resize(), libMesh::DenseMatrix< T >::resize(), and libMesh::DofMap::variable_type().

197 {
198  libmesh_ignore(u, R, J); // if we don't use --enable-second
199 
200 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
201  // Declare a performance log. Give it a descriptive
202  // string to identify what part of the code we are
203  // logging, since there may be many PerfLogs in an
204  // application.
205  PerfLog perf_log ("Biharmonic Residual and Jacobian", false);
206 
207  // A reference to the DofMap object for this system. The DofMap
208  // object handles the index translation from node and element numbers
209  // to degree of freedom numbers. We will talk more about the DofMap
210  // in future examples.
211  const DofMap & dof_map = get_dof_map();
212 
213  // Get a constant reference to the Finite Element type
214  // for the first (and only) variable in the system.
215  FEType fe_type = dof_map.variable_type(0);
216 
217  // Build a Finite Element object of the specified type. Since the
218  // FEBase::build() member dynamically creates memory we will
219  // store the object as a std::unique_ptr<FEBase>. This can be thought
220  // of as a pointer that will clean up after itself.
221  std::unique_ptr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
222 
223  // Quadrature rule for numerical integration.
224  // With 2D triangles, the Clough quadrature rule puts a Gaussian
225  // quadrature rule on each of the 3 subelements
226  std::unique_ptr<QBase> qrule(fe_type.default_quadrature_rule(_biharmonic._dim));
227 
228  // Tell the finite element object to use our quadrature rule.
229  fe->attach_quadrature_rule (qrule.get());
230 
231  // Here we define some references to element-specific data that
232  // will be used to assemble the linear system.
233  // We begin with the element Jacobian * quadrature weight at each
234  // integration point.
235  const std::vector<Real> & JxW = fe->get_JxW();
236 
237  // The element shape functions evaluated at the quadrature points.
238  const std::vector<std::vector<Real>> & phi = fe->get_phi();
239 
240  // The element shape functions' derivatives evaluated at the quadrature points.
241  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
242 
243  // The element shape functions' second derivatives evaluated at the quadrature points.
244  const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
245 
246  // For efficiency we will compute shape function laplacians n times,
247  // not n^2
248  std::vector<Real> Laplacian_phi_qp;
249 
250  // Define data structures to contain the element matrix
251  // and right-hand-side vector contribution. Following
252  // basic finite element terminology we will denote these
253  // "Je" and "Re". More detail is in example 3.
256 
257  // This vector will hold the degree of freedom indices for
258  // the element. These define where in the global system
259  // the element degrees of freedom get mapped.
260  std::vector<dof_id_type> dof_indices;
261 
262  // Old solution
263  const NumericVector<Number> & u_old = *old_local_solution;
264 
265  // Now we will loop over all the elements in the mesh. We will
266  // compute the element matrix and right-hand-side contribution. See
267  // example 3 for a discussion of the element iterators.
268  for (const auto & elem : _biharmonic._mesh.active_local_element_ptr_range())
269  {
270  // Get the degree of freedom indices for the
271  // current element. These define where in the global
272  // matrix and right-hand-side this element will
273  // contribute to.
274  dof_map.dof_indices (elem, dof_indices);
275 
276  const unsigned int n_dofs =
277  cast_int<unsigned int>(dof_indices.size());
278 
279  // Compute the element-specific data for the current
280  // element. This involves computing the location of the
281  // quadrature points (q_point) and the shape function
282  // values/derivatives (phi, dphi, d2phi) for the current element.
283  fe->reinit (elem);
284 
285  // Zero the element matrix, the right-hand side and the Laplacian matrix
286  // before summing them.
287  if (J)
288  Je.resize(n_dofs, n_dofs);
289 
290  if (R)
291  Re.resize(n_dofs);
292 
293  Laplacian_phi_qp.resize(n_dofs);
294 
295  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
296  {
297  // AUXILIARY QUANTITIES:
298  // Residual and Jacobian share a few calculations:
299  // at the very least, in the case of interfacial energy only with a constant mobility,
300  // both calculations use Laplacian_phi_qp; more is shared the case of a concentration-dependent
301  // mobility and bulk potentials.
302  Number
303  u_qp = 0.0,
304  u_old_qp = 0.0,
305  Laplacian_u_qp = 0.0,
306  Laplacian_u_old_qp = 0.0;
307 
308  Gradient
309  grad_u_qp(0.0, 0.0, 0.0),
310  grad_u_old_qp(0.0, 0.0, 0.0);
311 
312  Number
313  M_qp = 1.0,
314  M_old_qp = 1.0,
315  M_prime_qp = 0.0,
316  M_prime_old_qp = 0.0;
317 
318  for (unsigned int i=0; i<n_dofs; i++)
319  {
320  Laplacian_phi_qp[i] = d2phi[i][qp](0, 0);
321  grad_u_qp(0) += u(dof_indices[i])*dphi[i][qp](0);
322  grad_u_old_qp(0) += u_old(dof_indices[i])*dphi[i][qp](0);
323 
324  if (_biharmonic._dim > 1)
325  {
326  Laplacian_phi_qp[i] += d2phi[i][qp](1, 1);
327  grad_u_qp(1) += u(dof_indices[i])*dphi[i][qp](1);
328  grad_u_old_qp(1) += u_old(dof_indices[i])*dphi[i][qp](1);
329  }
330  if (_biharmonic._dim > 2)
331  {
332  Laplacian_phi_qp[i] += d2phi[i][qp](2, 2);
333  grad_u_qp(2) += u(dof_indices[i])*dphi[i][qp](2);
334  grad_u_old_qp(2) += u_old(dof_indices[i])*dphi[i][qp](2);
335  }
336  u_qp += phi[i][qp]*u(dof_indices[i]);
337  u_old_qp += phi[i][qp]*u_old(dof_indices[i]);
338  Laplacian_u_qp += Laplacian_phi_qp[i]*u(dof_indices[i]);
339  Laplacian_u_old_qp += Laplacian_phi_qp[i]*u_old(dof_indices[i]);
340  } // for i
341 
343  {
344  M_qp = 1.0 - u_qp*u_qp;
345  M_old_qp = 1.0 - u_old_qp*u_old_qp;
346  M_prime_qp = -2.0*u_qp;
347  M_prime_old_qp = -2.0*u_old_qp;
348  }
349 
350  // ELEMENT RESIDUAL AND JACOBIAN
351  for (unsigned int i=0; i<n_dofs; i++)
352  {
353  // RESIDUAL
354  if (R)
355  {
356  Number ri = 0.0, ri_old = 0.0;
357  ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_u_qp;
358  ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._kappa*Laplacian_u_old_qp;
359 
361  {
362  ri -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp);
363  ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*(_biharmonic._kappa*Laplacian_u_old_qp);
364  }
365 
367  {
369  {
370  ri += Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
371  ri_old += Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
373  {
374  ri += (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
375  ri_old += (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
376  }
377  }// if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
378 
380  {
381  ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*u_qp;
382  ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*u_old_qp;
384  {
385  ri -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*u_qp;
386  ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*u_old_qp;
387  }
388  } // if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
389 
391  {
393  {
394  case 2:
395  break;
396  case 3:
397  break;
398  default:
399  break;
400  }// switch(_biharmonic._log_truncation)
401  }// if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
402  }// if (_biharmonic._cahn_hillard)
403  Re(i) += JxW[qp]*((u_qp-u_old_qp)*phi[i][qp]-_biharmonic._dt*0.5*((2.0-_biharmonic._cnWeight)*ri + _biharmonic._cnWeight*ri_old));
404  } // if (R)
405 
406  // JACOBIAN
407  if (J)
408  {
409  Number M_prime_prime_qp = 0.0;
410  if (_biharmonic._degenerate) M_prime_prime_qp = -2.0;
411  for (unsigned int j=0; j<n_dofs; j++)
412  {
413  Number ri_j = 0.0;
414  ri_j -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_phi_qp[j];
416  {
417  ri_j -=
418  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._kappa*Laplacian_u_qp +
419  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp) +
420  (dphi[i][qp]*grad_u_qp)*(M_prime_prime_qp*phi[j][qp])*(_biharmonic._kappa*Laplacian_u_qp) +
421  (dphi[i][qp]*grad_u_qp)*(M_prime_qp)*(_biharmonic._kappa*Laplacian_phi_qp[j]);
422  }
423 
425  {
427  {
428  ri_j +=
429  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
430  Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp] +
431  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
432  (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
433  (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp];
434  }// if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
435 
437  {
438  ri_j -=
439  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*u_qp +
440  Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*phi[j][qp] +
441  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*u_qp +
442  (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*u_qp +
443  (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*phi[j][qp];
444  } // if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
445 
447  {
449  {
450  case 2:
451  break;
452  case 3:
453  break;
454  default:
455  break;
456  }// switch(_biharmonic._log_truncation)
457  }// if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
458  }// if (_biharmonic._cahn_hillard)
459  Je(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] - 0.5*_biharmonic._dt*(2.0-_biharmonic._cnWeight)*ri_j);
460  } // for j
461  } // if (J)
462  } // for i
463  } // for qp
464 
465  // The element matrix and right-hand-side are now built
466  // for this element. Add them to the global matrix and
467  // right-hand-side vector. The SparseMatrix::add_matrix()
468  // and NumericVector::add_vector() members do this for us.
469  // Start logging the insertion of the local (element)
470  // matrix and vector into the global matrix and vector
471  if (R)
472  {
473  // If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
474  dof_map.constrain_element_vector(Re, dof_indices);
475  R->add_vector(Re, dof_indices);
476  }
477 
478  if (J)
479  {
480  if (R)
481  // The call to constrain_element_vector can modify the dof_indices based on constraints,
482  // but the input to constrain_element_matrix relies on using the unmodified indices, so
483  // we reset them here
484  dof_map.dof_indices(elem, dof_indices);
485 
486  // If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
487  dof_map.constrain_element_matrix(Je, dof_indices);
488  J->add_matrix(Je, dof_indices);
489  }
490  } // for el
491 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
492 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
Biharmonic & _biharmonic
Definition: biharmonic_jr.h:91
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1992
std::size_t n_dofs() const
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:374
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2144
Real _cnWeight
Definition: biharmonic.h:103
Real _kappa
Definition: biharmonic.h:93
Real _theta_c
Definition: biharmonic.h:93
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:145
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:169
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:34
bool _cahn_hillard
Definition: biharmonic.h:95
void libmesh_ignore(const Args &...)
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
FreeEnergyEnum _energy
Definition: biharmonic.h:96
int _log_truncation
Definition: biharmonic.h:97
unsigned int _dim
Definition: biharmonic.h:92
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
Constrains the element vector.
Definition: dof_map.h:2250
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
ReplicatedMesh & _mesh
Definition: biharmonic.h:111
bool _degenerate
Definition: biharmonic.h:95
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix.
Definition: dof_map.h:2241

Member Data Documentation

◆ _biharmonic

Biharmonic& Biharmonic::JR::_biharmonic
private

Definition at line 91 of file biharmonic_jr.h.

Referenced by JR().


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