libMesh
biharmonic_jr.C
Go to the documentation of this file.
1 // Libmesh includes
2 #include "libmesh/mesh.h"
3 #include "libmesh/quadrature.h"
4 #include "libmesh/dense_matrix.h"
5 #include "libmesh/dense_vector.h"
6 #include "libmesh/sparse_matrix.h"
7 #include "libmesh/fourth_error_estimators.h"
8 #include "libmesh/dof_map.h"
9 #include "libmesh/numeric_vector.h"
10 #include "libmesh/periodic_boundaries.h"
11 #include "libmesh/periodic_boundary.h"
12 #include "libmesh/elem.h"
13 
14 // Example includes
15 #include "biharmonic_jr.h"
16 
17 using namespace libMesh;
18 
20  const std::string & name_in,
21  const unsigned int number_in) :
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  // Attache the R & J calculation object
84  nonlinear_solver->residual_and_jacobian_object = this;
85 
86  // Attach the bounds calculation object
87  nonlinear_solver->bounds_object = this;
88 }
89 
90 
91 
92 
93 
95 {
97  libMesh::out << ">>> Initializing Biharmonic::JR\n";
98 
100  parameters.set<Point>("center") = _biharmonic._initialCenter;
101  parameters.set<Real>("width") = _biharmonic._initialWidth;
102 
105 
108 
111 
112  // both states are equal
113  *(old_local_solution) = *(current_local_solution);
114 
115  if (_biharmonic._verbose)
116  libMesh::out << "<<< Initializing Biharmonic::JR\n";
117 }
118 
119 
120 
121 
122 
123 
125  const Parameters & parameters,
126  const std::string &,
127  const std::string &)
128 {
129  // Initialize with a ball in the middle, which is a segment in 1D, a disk in 2D and a ball in 3D.
130  Point center = parameters.get<Point>("center");
131  Real width = parameters.get<Real>("width");
132  Point pc = p-center;
133  Real r = pc.norm();
134  return (r < width) ? 1.0 : -0.5;
135 }
136 
137 
138 
139 
141  const Parameters & parameters,
142  const std::string &,
143  const std::string &)
144 {
145  // Initialize with a rod in the middle so that we have a z-homogeneous system to model the 2D disk.
146  Point center = parameters.get<Point>("center");
147  Real width = parameters.get<Real>("width");
148  Real r = sqrt((p(0)-center(0))*(p(0)-center(0)) + (p(1)-center(1))*(p(1)-center(1)));
149  return (r < width) ? 1.0 : -0.5;
150 }
151 
152 
153 
154 
155 
157  const Parameters & parameters,
158  const std::string &,
159  const std::string &)
160 {
161  // Initialize with a wide strip in the middle so that we have a yz-homogeneous system to model the 1D.
162  Point center = parameters.get<Point>("center");
163  Real width = parameters.get<Real>("width");
164  Real r = sqrt((p(0)-center(0))*(p(0)-center(0)));
165  return (r < width) ? 1.0 : -0.5;
166 }
167 
168 
169 
170 
172  const Parameters &,
173  const std::string &,
174  const std::string &)
175 {
176  return Gradient(0.0, 0.0, 0.0);
177 }
178 
179 
180 
181 
186 {
187 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
188  if (!R && !J)
189  return;
190 
191  // Declare a performance log. Give it a descriptive
192  // string to identify what part of the code we are
193  // logging, since there may be many PerfLogs in an
194  // application.
195  PerfLog perf_log ("Biharmonic Residual and Jacobian", false);
196 
197  // A reference to the DofMap object for this system. The DofMap
198  // object handles the index translation from node and element numbers
199  // to degree of freedom numbers. We will talk more about the DofMap
200  // in future examples.
201  const DofMap & dof_map = get_dof_map();
202 
203  // Get a constant reference to the Finite Element type
204  // for the first (and only) variable in the system.
205  FEType fe_type = dof_map.variable_type(0);
206 
207  // Build a Finite Element object of the specified type. Since the
208  // FEBase::build() member dynamically creates memory we will
209  // store the object as a UniquePtr<FEBase>. This can be thought
210  // of as a pointer that will clean up after itself.
211  UniquePtr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
212 
213  // Quadrature rule for numerical integration.
214  // With 2D triangles, the Clough quadrature rule puts a Gaussian
215  // quadrature rule on each of the 3 subelements
217 
218  // Tell the finite element object to use our quadrature rule.
219  fe->attach_quadrature_rule (qrule.get());
220 
221  // Here we define some references to element-specific data that
222  // will be used to assemble the linear system.
223  // We begin with the element Jacobian * quadrature weight at each
224  // integration point.
225  const std::vector<Real> & JxW = fe->get_JxW();
226 
227  // The element shape functions evaluated at the quadrature points.
228  const std::vector<std::vector<Real>> & phi = fe->get_phi();
229 
230  // The element shape functions' derivatives evaluated at the quadrature points.
231  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
232 
233  // The element shape functions' second derivatives evaluated at the quadrature points.
234  const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
235 
236  // For efficiency we will compute shape function laplacians n times,
237  // not n^2
238  std::vector<Real> Laplacian_phi_qp;
239 
240  // Define data structures to contain the element matrix
241  // and right-hand-side vector contribution. Following
242  // basic finite element terminology we will denote these
243  // "Je" and "Re". More detail is in example 3.
246 
247  // This vector will hold the degree of freedom indices for
248  // the element. These define where in the global system
249  // the element degrees of freedom get mapped.
250  std::vector<dof_id_type> dof_indices;
251 
252  // Old solution
253  const NumericVector<Number> & u_old = *old_local_solution;
254 
255  // Now we will loop over all the elements in the mesh. We will
256  // compute the element matrix and right-hand-side contribution. See
257  // example 3 for a discussion of the element iterators.
258  for (const auto & elem : _biharmonic._mesh.active_local_element_ptr_range())
259  {
260  // Get the degree of freedom indices for the
261  // current element. These define where in the global
262  // matrix and right-hand-side this element will
263  // contribute to.
264  dof_map.dof_indices (elem, dof_indices);
265 
266  // Compute the element-specific data for the current
267  // element. This involves computing the location of the
268  // quadrature points (q_point) and the shape function
269  // values/derivatives (phi, dphi, d2phi) for the current element.
270  fe->reinit (elem);
271 
272  // Zero the element matrix, the right-hand side and the Laplacian matrix
273  // before summing them.
274  if (J)
275  Je.resize(dof_indices.size(), dof_indices.size());
276 
277  if (R)
278  Re.resize(dof_indices.size());
279 
280  Laplacian_phi_qp.resize(dof_indices.size());
281 
282  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
283  {
284  // AUXILIARY QUANTITIES:
285  // Residual and Jacobian share a few calculations:
286  // at the very least, in the case of interfacial energy only with a constant mobility,
287  // both calculations use Laplacian_phi_qp; more is shared the case of a concentration-dependent
288  // mobility and bulk potentials.
289  Number
290  u_qp = 0.0,
291  u_old_qp = 0.0,
292  Laplacian_u_qp = 0.0,
293  Laplacian_u_old_qp = 0.0;
294 
295  Gradient
296  grad_u_qp(0.0, 0.0, 0.0),
297  grad_u_old_qp(0.0, 0.0, 0.0);
298 
299  Number
300  M_qp = 1.0,
301  M_old_qp = 1.0,
302  M_prime_qp = 0.0,
303  M_prime_old_qp = 0.0;
304 
305  for (std::size_t i=0; i<phi.size(); i++)
306  {
307  Laplacian_phi_qp[i] = d2phi[i][qp](0, 0);
308  grad_u_qp(0) += u(dof_indices[i])*dphi[i][qp](0);
309  grad_u_old_qp(0) += u_old(dof_indices[i])*dphi[i][qp](0);
310 
311  if (_biharmonic._dim > 1)
312  {
313  Laplacian_phi_qp[i] += d2phi[i][qp](1, 1);
314  grad_u_qp(1) += u(dof_indices[i])*dphi[i][qp](1);
315  grad_u_old_qp(1) += u_old(dof_indices[i])*dphi[i][qp](1);
316  }
317  if (_biharmonic._dim > 2)
318  {
319  Laplacian_phi_qp[i] += d2phi[i][qp](2, 2);
320  grad_u_qp(2) += u(dof_indices[i])*dphi[i][qp](2);
321  grad_u_old_qp(2) += u_old(dof_indices[i])*dphi[i][qp](2);
322  }
323  u_qp += phi[i][qp]*u(dof_indices[i]);
324  u_old_qp += phi[i][qp]*u_old(dof_indices[i]);
325  Laplacian_u_qp += Laplacian_phi_qp[i]*u(dof_indices[i]);
326  Laplacian_u_old_qp += Laplacian_phi_qp[i]*u_old(dof_indices[i]);
327  } // for i
328 
330  {
331  M_qp = 1.0 - u_qp*u_qp;
332  M_old_qp = 1.0 - u_old_qp*u_old_qp;
333  M_prime_qp = -2.0*u_qp;
334  M_prime_old_qp = -2.0*u_old_qp;
335  }
336 
337  // ELEMENT RESIDUAL AND JACOBIAN
338  for (std::size_t i=0; i<phi.size(); i++)
339  {
340  // RESIDUAL
341  if (R)
342  {
343  Number ri = 0.0, ri_old = 0.0;
344  ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_u_qp;
345  ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._kappa*Laplacian_u_old_qp;
346 
348  {
349  ri -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp);
350  ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*(_biharmonic._kappa*Laplacian_u_old_qp);
351  }
352 
354  {
356  {
357  ri += Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
358  ri_old += Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
360  {
361  ri += (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
362  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;
363  }
364  }// if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
365 
367  {
368  ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*u_qp;
369  ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*u_old_qp;
371  {
372  ri -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*u_qp;
373  ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*u_old_qp;
374  }
375  } // if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
376 
378  {
380  {
381  case 2:
382  break;
383  case 3:
384  break;
385  default:
386  break;
387  }// switch(_biharmonic._log_truncation)
388  }// if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
389  }// if (_biharmonic._cahn_hillard)
390  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));
391  } // if (R)
392 
393  // JACOBIAN
394  if (J)
395  {
396  Number M_prime_prime_qp = 0.0;
397  if (_biharmonic._degenerate) M_prime_prime_qp = -2.0;
398  for (std::size_t j=0; j<phi.size(); j++)
399  {
400  Number ri_j = 0.0;
401  ri_j -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_phi_qp[j];
403  {
404  ri_j -=
405  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._kappa*Laplacian_u_qp +
406  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp) +
407  (dphi[i][qp]*grad_u_qp)*(M_prime_prime_qp*phi[j][qp])*(_biharmonic._kappa*Laplacian_u_qp) +
408  (dphi[i][qp]*grad_u_qp)*(M_prime_qp)*(_biharmonic._kappa*Laplacian_phi_qp[j]);
409  }
410 
412  {
414  {
415  ri_j +=
416  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
417  Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp] +
418  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
419  (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
420  (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp];
421  }// if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
422 
424  {
425  ri_j -=
426  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*u_qp +
427  Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*phi[j][qp] +
428  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*u_qp +
429  (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*u_qp +
430  (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*phi[j][qp];
431  } // if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
432 
434  {
436  {
437  case 2:
438  break;
439  case 3:
440  break;
441  default:
442  break;
443  }// switch(_biharmonic._log_truncation)
444  }// if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
445  }// if (_biharmonic._cahn_hillard)
446  Je(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] - 0.5*_biharmonic._dt*(2.0-_biharmonic._cnWeight)*ri_j);
447  } // for j
448  } // if (J)
449  } // for i
450  } // for qp
451 
452  // The element matrix and right-hand-side are now built
453  // for this element. Add them to the global matrix and
454  // right-hand-side vector. The SparseMatrix::add_matrix()
455  // and NumericVector::add_vector() members do this for us.
456  // Start logging the insertion of the local (element)
457  // matrix and vector into the global matrix and vector
458  if (R)
459  {
460  // If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
461  dof_map.constrain_element_vector(Re, dof_indices);
462  R->add_vector(Re, dof_indices);
463  }
464 
465  if (J)
466  {
467  // If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
468  dof_map.constrain_element_matrix(Je, dof_indices);
469  J->add_matrix(Je, dof_indices);
470  }
471  } // for el
472 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
473 }
474 
475 
476 
477 
478 
482 {
483  // sys is actually ignored, since it should be the same as *this.
484 
485  // Declare a performance log. Give it a descriptive
486  // string to identify what part of the code we are
487  // logging, since there may be many PerfLogs in an
488  // application.
489  PerfLog perf_log ("Biharmonic bounds", false);
490 
491  // A reference to the DofMap object for this system. The DofMap
492  // object handles the index translation from node and element numbers
493  // to degree of freedom numbers. We will talk more about the DofMap
494  // in future examples.
495  const DofMap & dof_map = get_dof_map();
496 
497  // Get a constant reference to the Finite Element type
498  // for the first (and only) variable in the system.
499  FEType fe_type = dof_map.variable_type(0);
500 
501  // Build a Finite Element object of the specified type. Since the
502  // FEBase::build() member dynamically creates memory we will
503  // store the object as a UniquePtr<FEBase>. This can be thought
504  // of as a pointer that will clean up after itself.
505  UniquePtr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
506 
507  // Define data structures to contain the bound vectors contributions.
508  DenseVector<Number> XLe, XUe;
509 
510  // These vector will hold the degree of freedom indices for
511  // the element. These define where in the global system
512  // the element degrees of freedom get mapped.
513  std::vector<dof_id_type> dof_indices;
514 
515  for (const auto & elem : _biharmonic._mesh.active_local_element_ptr_range())
516  {
517  // Extract the shape function to be evaluated at the nodes
518  const std::vector<std::vector<Real>> & phi = fe->get_phi();
519 
520  // Get the degree of freedom indices for the current element.
521  // They are in 1-1 correspondence with shape functions phi
522  // and define where in the global vector this element will.
523  dof_map.dof_indices (elem, dof_indices);
524 
525  // Resize the local bounds vectors (zeroing them out in the process).
526  XLe.resize(dof_indices.size());
527  XUe.resize(dof_indices.size());
528 
529  // Extract the element node coordinates in the reference frame
530  std::vector<Point> nodes;
531  fe->get_refspace_nodes(elem->type(), nodes);
532 
533  // Evaluate the shape functions at the nodes
534  fe->reinit(elem, &nodes);
535 
536  // Construct the bounds based on the value of the i-th phi at the nodes.
537  // Observe that this doesn't really work in general: we rely on the fact
538  // that for Hermite elements each shape function is nonzero at most at a
539  // single node.
540  // More generally the bounds must be constructed by inspecting a "mass-like"
541  // matrix (m_{ij}) of the shape functions (i) evaluated at their corresponding nodes (j).
542  // The constraints imposed on the dofs (d_i) are then are -1 \leq \sum_i d_i m_{ij} \leq 1,
543  // since \sum_i d_i m_{ij} is the value of the solution at the j-th node.
544  // Auxiliary variables will need to be introduced to reduce this to a "box" constraint.
545  // Additional complications will arise since m might be singular (as is the case for Hermite,
546  // which, however, is easily handled by inspection).
547  for (std::size_t i=0; i<phi.size(); ++i)
548  {
549  // FIXME: should be able to define INF and pass it to the solve
550  Real infinity = 1.0e20;
551  Real bound = infinity;
552  for (std::size_t j = 0; j < nodes.size(); ++j)
553  {
554  if (phi[i][j])
555  {
556  bound = 1.0/std::abs(phi[i][j]);
557  break;
558  }
559  }
560 
561  // The value of the solution at this node must be between 1.0 and -1.0.
562  // Based on the value of phi(i)(i) the nodal coordinate must be between 1.0/phi(i)(i) and its negative.
563  XLe(i) = -bound;
564  XUe(i) = bound;
565  }
566  // The element bound vectors are now built for this element.
567  // Insert them into the global vectors, potentially overwriting
568  // the same dof contributions from other elements: no matter --
569  // the bounds are always -1.0 and 1.0.
570  XL.insert(XLe, dof_indices);
571  XU.insert(XUe, dof_indices);
572  }
573 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
virtual SimpleRange< element_iterator > active_local_element_ptr_range() libmesh_override
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.
double abs(double a)
Biharmonic & _biharmonic
Definition: biharmonic_jr.h:83
This is the EquationSystems class.
void add_periodic_boundary(const PeriodicBoundaryBase &periodic_boundary)
Adds a copy of the specified periodic boundary to the system.
JR(EquationSystems &eqSys, const std::string &name, const unsigned int number)
Constructor.
Definition: biharmonic_jr.C:19
This class provides the ability to map between arbitrary, user-defined strings and several data types...
Definition: parameters.h:63
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
VectorValue< Real > RealVectorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
bool _verbose
Definition: biharmonic.h:99
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
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]...
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:1802
Real _cnWeight
Definition: biharmonic.h:104
Real _kappa
Definition: biharmonic.h:94
Real _theta_c
Definition: biharmonic.h:94
This class provides a specific system class.
const T & get(const std::string &) const
Definition: parameters.h:430
The Biharmonic class encapsulates most of the data structures necessary to calculate the biharmonic r...
Definition: biharmonic.h:46
void residual_and_jacobian(const NumericVector< Number > &u, NumericVector< Number > *R, SparseMatrix< Number > *J, NonlinearImplicitSystem &)
The residual and Jacobian assembly function for the Biharmonic system.
The libMesh namespace provides an interface to certain functionality in the library.
PetscErrorCode Vec Mat Mat pc
The definition of a periodic boundary.
Point _initialCenter
Definition: biharmonic.h:101
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
The PerfLog class allows monitoring of specific events.
Definition: perf_log.h:125
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
bool _cahn_hillard
Definition: biharmonic.h:96
boundary_id_type myboundary
The boundary ID of this boundary and its counterpart.
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.
static Gradient InitialGradientZero(const Point &, const Parameters &, const std::string &, const std::string &)
NumberVectorValue Gradient
This class provides a specific system class.
static Number InitialDensityRod(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
FreeEnergyEnum _energy
Definition: biharmonic.h:97
int _log_truncation
Definition: biharmonic.h:98
unsigned int _dim
Definition: biharmonic.h:93
static Number InitialDensityStrip(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
void bounds(NumericVector< Number > &XL, NumericVector< Number > &XU, NonlinearImplicitSystem &)
Function defining the bounds of the Biharmonic system.
Real _initialWidth
Definition: biharmonic.h:102
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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:1793
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
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
OStreamProxy out
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
Parameters parameters
Data structure holding arbitrary parameters.
ReplicatedMesh & _mesh
Definition: biharmonic.h:112
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
InitialStateEnum _initialState
Definition: biharmonic.h:100
sys get_dof_map()
bool _degenerate
Definition: biharmonic.h:96
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:1917