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
JR(EquationSystems & eqSys,const std::string & name_in,const unsigned int number_in)19 Biharmonic::JR::JR(EquationSystems & eqSys,
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
initialize()94 void Biharmonic::JR::initialize()
95 {
96 if (_biharmonic._verbose)
97 libMesh::out << ">>> Initializing Biharmonic::JR\n";
98
99 Parameters parameters;
100 parameters.set<Point>("center") = _biharmonic._initialCenter;
101 parameters.set<Real>("width") = _biharmonic._initialWidth;
102
103 if (_biharmonic._initialState == Biharmonic::BALL)
104 project_solution(Biharmonic::JR::InitialDensityBall, Biharmonic::JR::InitialGradientZero, parameters);
105
106 if (_biharmonic._initialState == Biharmonic::ROD)
107 project_solution(Biharmonic::JR::InitialDensityRod, Biharmonic::JR::InitialGradientZero, parameters);
108
109 if (_biharmonic._initialState == Biharmonic::STRIP)
110 project_solution(Biharmonic::JR::InitialDensityStrip, Biharmonic::JR::InitialGradientZero, parameters);
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
InitialDensityBall(const Point & p,const Parameters & parameters,const std::string &,const std::string &)124 Number Biharmonic::JR::InitialDensityBall(const Point & p,
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
InitialDensityRod(const Point & p,const Parameters & parameters,const std::string &,const std::string &)140 Number Biharmonic::JR::InitialDensityRod(const Point & p,
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 Point pc = p-center;
149 #if LIBMESH_DIM > 2
150 pc(2) = 0;
151 #endif
152 Real r = pc.norm();
153 return (r < width) ? 1.0 : -0.5;
154 }
155
156
157
158
159
InitialDensityStrip(const Point & p,const Parameters & parameters,const std::string &,const std::string &)160 Number Biharmonic::JR::InitialDensityStrip(const Point & p,
161 const Parameters & parameters,
162 const std::string &,
163 const std::string &)
164 {
165 // Initialize with a wide strip in the middle so that we have a yz-homogeneous system to model the 1D.
166 Point center = parameters.get<Point>("center");
167 Real width = parameters.get<Real>("width");
168 Real r = sqrt((p(0)-center(0))*(p(0)-center(0)));
169 return (r < width) ? 1.0 : -0.5;
170 }
171
172
173
174
InitialGradientZero(const Point &,const Parameters &,const std::string &,const std::string &)175 Gradient Biharmonic::JR::InitialGradientZero(const Point &,
176 const Parameters &,
177 const std::string &,
178 const std::string &)
179 {
180 return Gradient(0.0, 0.0, 0.0);
181 }
182
183
184
185
residual_and_jacobian(const NumericVector<Number> & u,NumericVector<Number> * R,SparseMatrix<Number> * J,NonlinearImplicitSystem &)186 void Biharmonic::JR::residual_and_jacobian(const NumericVector<Number> & u,
187 NumericVector<Number> * R,
188 SparseMatrix<Number> * J,
189 NonlinearImplicitSystem &)
190 {
191 libmesh_ignore(u, R, J); // if we don't use --enable-second
192
193 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
194 if (!R && !J)
195 return;
196
197 // Declare a performance log. Give it a descriptive
198 // string to identify what part of the code we are
199 // logging, since there may be many PerfLogs in an
200 // application.
201 PerfLog perf_log ("Biharmonic Residual and Jacobian", false);
202
203 // A reference to the DofMap object for this system. The DofMap
204 // object handles the index translation from node and element numbers
205 // to degree of freedom numbers. We will talk more about the DofMap
206 // in future examples.
207 const DofMap & dof_map = get_dof_map();
208
209 // Get a constant reference to the Finite Element type
210 // for the first (and only) variable in the system.
211 FEType fe_type = dof_map.variable_type(0);
212
213 // Build a Finite Element object of the specified type. Since the
214 // FEBase::build() member dynamically creates memory we will
215 // store the object as a std::unique_ptr<FEBase>. This can be thought
216 // of as a pointer that will clean up after itself.
217 std::unique_ptr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
218
219 // Quadrature rule for numerical integration.
220 // With 2D triangles, the Clough quadrature rule puts a Gaussian
221 // quadrature rule on each of the 3 subelements
222 std::unique_ptr<QBase> qrule(fe_type.default_quadrature_rule(_biharmonic._dim));
223
224 // Tell the finite element object to use our quadrature rule.
225 fe->attach_quadrature_rule (qrule.get());
226
227 // Here we define some references to element-specific data that
228 // will be used to assemble the linear system.
229 // We begin with the element Jacobian * quadrature weight at each
230 // integration point.
231 const std::vector<Real> & JxW = fe->get_JxW();
232
233 // The element shape functions evaluated at the quadrature points.
234 const std::vector<std::vector<Real>> & phi = fe->get_phi();
235
236 // The element shape functions' derivatives evaluated at the quadrature points.
237 const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
238
239 // The element shape functions' second derivatives evaluated at the quadrature points.
240 const std::vector<std::vector<RealTensor>> & d2phi = fe->get_d2phi();
241
242 // For efficiency we will compute shape function laplacians n times,
243 // not n^2
244 std::vector<Real> Laplacian_phi_qp;
245
246 // Define data structures to contain the element matrix
247 // and right-hand-side vector contribution. Following
248 // basic finite element terminology we will denote these
249 // "Je" and "Re". More detail is in example 3.
250 DenseMatrix<Number> Je;
251 DenseVector<Number> Re;
252
253 // This vector will hold the degree of freedom indices for
254 // the element. These define where in the global system
255 // the element degrees of freedom get mapped.
256 std::vector<dof_id_type> dof_indices;
257
258 // Old solution
259 const NumericVector<Number> & u_old = *old_local_solution;
260
261 // Now we will loop over all the elements in the mesh. We will
262 // compute the element matrix and right-hand-side contribution. See
263 // example 3 for a discussion of the element iterators.
264 for (const auto & elem : _biharmonic._mesh.active_local_element_ptr_range())
265 {
266 // Get the degree of freedom indices for the
267 // current element. These define where in the global
268 // matrix and right-hand-side this element will
269 // contribute to.
270 dof_map.dof_indices (elem, dof_indices);
271
272 const unsigned int n_dofs =
273 cast_int<unsigned int>(dof_indices.size());
274
275 // Compute the element-specific data for the current
276 // element. This involves computing the location of the
277 // quadrature points (q_point) and the shape function
278 // values/derivatives (phi, dphi, d2phi) for the current element.
279 fe->reinit (elem);
280
281 // Zero the element matrix, the right-hand side and the Laplacian matrix
282 // before summing them.
283 if (J)
284 Je.resize(n_dofs, n_dofs);
285
286 if (R)
287 Re.resize(n_dofs);
288
289 Laplacian_phi_qp.resize(n_dofs);
290
291 for (unsigned int qp=0; qp<qrule->n_points(); qp++)
292 {
293 // AUXILIARY QUANTITIES:
294 // Residual and Jacobian share a few calculations:
295 // at the very least, in the case of interfacial energy only with a constant mobility,
296 // both calculations use Laplacian_phi_qp; more is shared the case of a concentration-dependent
297 // mobility and bulk potentials.
298 Number
299 u_qp = 0.0,
300 u_old_qp = 0.0,
301 Laplacian_u_qp = 0.0,
302 Laplacian_u_old_qp = 0.0;
303
304 Gradient
305 grad_u_qp(0.0, 0.0, 0.0),
306 grad_u_old_qp(0.0, 0.0, 0.0);
307
308 Number
309 M_qp = 1.0,
310 M_old_qp = 1.0,
311 M_prime_qp = 0.0,
312 M_prime_old_qp = 0.0;
313
314 for (unsigned int i=0; i<n_dofs; i++)
315 {
316 Laplacian_phi_qp[i] = d2phi[i][qp](0, 0);
317 grad_u_qp(0) += u(dof_indices[i])*dphi[i][qp](0);
318 grad_u_old_qp(0) += u_old(dof_indices[i])*dphi[i][qp](0);
319
320 if (_biharmonic._dim > 1)
321 {
322 Laplacian_phi_qp[i] += d2phi[i][qp](1, 1);
323 grad_u_qp(1) += u(dof_indices[i])*dphi[i][qp](1);
324 grad_u_old_qp(1) += u_old(dof_indices[i])*dphi[i][qp](1);
325 }
326 if (_biharmonic._dim > 2)
327 {
328 Laplacian_phi_qp[i] += d2phi[i][qp](2, 2);
329 grad_u_qp(2) += u(dof_indices[i])*dphi[i][qp](2);
330 grad_u_old_qp(2) += u_old(dof_indices[i])*dphi[i][qp](2);
331 }
332 u_qp += phi[i][qp]*u(dof_indices[i]);
333 u_old_qp += phi[i][qp]*u_old(dof_indices[i]);
334 Laplacian_u_qp += Laplacian_phi_qp[i]*u(dof_indices[i]);
335 Laplacian_u_old_qp += Laplacian_phi_qp[i]*u_old(dof_indices[i]);
336 } // for i
337
338 if (_biharmonic._degenerate)
339 {
340 M_qp = 1.0 - u_qp*u_qp;
341 M_old_qp = 1.0 - u_old_qp*u_old_qp;
342 M_prime_qp = -2.0*u_qp;
343 M_prime_old_qp = -2.0*u_old_qp;
344 }
345
346 // ELEMENT RESIDUAL AND JACOBIAN
347 for (unsigned int i=0; i<n_dofs; i++)
348 {
349 // RESIDUAL
350 if (R)
351 {
352 Number ri = 0.0, ri_old = 0.0;
353 ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_u_qp;
354 ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._kappa*Laplacian_u_old_qp;
355
356 if (_biharmonic._degenerate)
357 {
358 ri -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp);
359 ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*(_biharmonic._kappa*Laplacian_u_old_qp);
360 }
361
362 if (_biharmonic._cahn_hillard)
363 {
364 if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
365 {
366 ri += Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
367 ri_old += Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
368 if (_biharmonic._degenerate)
369 {
370 ri += (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
371 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;
372 }
373 }// if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
374
375 if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
376 {
377 ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*u_qp;
378 ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*u_old_qp;
379 if (_biharmonic._degenerate)
380 {
381 ri -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*u_qp;
382 ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*u_old_qp;
383 }
384 } // if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
385
386 if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
387 {
388 switch(_biharmonic._log_truncation)
389 {
390 case 2:
391 break;
392 case 3:
393 break;
394 default:
395 break;
396 }// switch(_biharmonic._log_truncation)
397 }// if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
398 }// if (_biharmonic._cahn_hillard)
399 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));
400 } // if (R)
401
402 // JACOBIAN
403 if (J)
404 {
405 Number M_prime_prime_qp = 0.0;
406 if (_biharmonic._degenerate) M_prime_prime_qp = -2.0;
407 for (unsigned int j=0; j<n_dofs; j++)
408 {
409 Number ri_j = 0.0;
410 ri_j -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_phi_qp[j];
411 if (_biharmonic._degenerate)
412 {
413 ri_j -=
414 Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._kappa*Laplacian_u_qp +
415 (dphi[i][qp]*dphi[j][qp])*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp) +
416 (dphi[i][qp]*grad_u_qp)*(M_prime_prime_qp*phi[j][qp])*(_biharmonic._kappa*Laplacian_u_qp) +
417 (dphi[i][qp]*grad_u_qp)*(M_prime_qp)*(_biharmonic._kappa*Laplacian_phi_qp[j]);
418 }
419
420 if (_biharmonic._cahn_hillard)
421 {
422 if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
423 {
424 ri_j +=
425 Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
426 Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp] +
427 (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
428 (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
429 (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp];
430 }// if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
431
432 if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
433 {
434 ri_j -=
435 Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*u_qp +
436 Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*phi[j][qp] +
437 (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*u_qp +
438 (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*u_qp +
439 (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*phi[j][qp];
440 } // if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
441
442 if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
443 {
444 switch(_biharmonic._log_truncation)
445 {
446 case 2:
447 break;
448 case 3:
449 break;
450 default:
451 break;
452 }// switch(_biharmonic._log_truncation)
453 }// if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
454 }// if (_biharmonic._cahn_hillard)
455 Je(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] - 0.5*_biharmonic._dt*(2.0-_biharmonic._cnWeight)*ri_j);
456 } // for j
457 } // if (J)
458 } // for i
459 } // for qp
460
461 // The element matrix and right-hand-side are now built
462 // for this element. Add them to the global matrix and
463 // right-hand-side vector. The SparseMatrix::add_matrix()
464 // and NumericVector::add_vector() members do this for us.
465 // Start logging the insertion of the local (element)
466 // matrix and vector into the global matrix and vector
467 if (R)
468 {
469 // If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
470 dof_map.constrain_element_vector(Re, dof_indices);
471 R->add_vector(Re, dof_indices);
472 }
473
474 if (J)
475 {
476 // If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
477 dof_map.constrain_element_matrix(Je, dof_indices);
478 J->add_matrix(Je, dof_indices);
479 }
480 } // for el
481 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
482 }
483
484
485
486
487
bounds(NumericVector<Number> & XL,NumericVector<Number> & XU,NonlinearImplicitSystem &)488 void Biharmonic::JR::bounds(NumericVector<Number> & XL,
489 NumericVector<Number> & XU,
490 NonlinearImplicitSystem &)
491 {
492 // sys is actually ignored, since it should be the same as *this.
493
494 // Declare a performance log. Give it a descriptive
495 // string to identify what part of the code we are
496 // logging, since there may be many PerfLogs in an
497 // application.
498 PerfLog perf_log ("Biharmonic bounds", false);
499
500 // A reference to the DofMap object for this system. The DofMap
501 // object handles the index translation from node and element numbers
502 // to degree of freedom numbers. We will talk more about the DofMap
503 // in future examples.
504 const DofMap & dof_map = get_dof_map();
505
506 // Get a constant reference to the Finite Element type
507 // for the first (and only) variable in the system.
508 FEType fe_type = dof_map.variable_type(0);
509
510 // Build a Finite Element object of the specified type. Since the
511 // FEBase::build() member dynamically creates memory we will
512 // store the object as a std::unique_ptr<FEBase>. This can be thought
513 // of as a pointer that will clean up after itself.
514 std::unique_ptr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
515
516 // Define data structures to contain the bound vectors contributions.
517 DenseVector<Number> XLe, XUe;
518
519 // These vector will hold the degree of freedom indices for
520 // the element. These define where in the global system
521 // the element degrees of freedom get mapped.
522 std::vector<dof_id_type> dof_indices;
523
524 for (const auto & elem : _biharmonic._mesh.active_local_element_ptr_range())
525 {
526 // Extract the shape function to be evaluated at the nodes
527 const std::vector<std::vector<Real>> & phi = fe->get_phi();
528
529 // Get the degree of freedom indices for the current element.
530 // They are in 1-1 correspondence with shape functions phi
531 // and define where in the global vector this element will.
532 dof_map.dof_indices (elem, dof_indices);
533
534 const unsigned int n_dofs =
535 cast_int<unsigned int>(dof_indices.size());
536
537 // Resize the local bounds vectors (zeroing them out in the process).
538 XLe.resize(n_dofs);
539 XUe.resize(n_dofs);
540
541 // Extract the element node coordinates in the reference frame
542 std::vector<Point> nodes;
543 fe->get_refspace_nodes(elem->type(), nodes);
544
545 // Evaluate the shape functions at the nodes
546 fe->reinit(elem, &nodes);
547
548 // Construct the bounds based on the value of the i-th phi at the nodes.
549 // Observe that this doesn't really work in general: we rely on the fact
550 // that for Hermite elements each shape function is nonzero at most at a
551 // single node.
552 // More generally the bounds must be constructed by inspecting a "mass-like"
553 // matrix (m_{ij}) of the shape functions (i) evaluated at their corresponding nodes (j).
554 // The constraints imposed on the dofs (d_i) are then are -1 \leq \sum_i d_i m_{ij} \leq 1,
555 // since \sum_i d_i m_{ij} is the value of the solution at the j-th node.
556 // Auxiliary variables will need to be introduced to reduce this to a "box" constraint.
557 // Additional complications will arise since m might be singular (as is the case for Hermite,
558 // which, however, is easily handled by inspection).
559 for (unsigned int i=0; i<n_dofs; ++i)
560 {
561 // FIXME: should be able to define INF and pass it to the solve
562 Real infinity = 1.0e20;
563 Real bound = infinity;
564 for (unsigned int j = 0; j < nodes.size(); ++j)
565 {
566 if (phi[i][j])
567 {
568 bound = 1.0/std::abs(phi[i][j]);
569 break;
570 }
571 }
572
573 // The value of the solution at this node must be between 1.0 and -1.0.
574 // Based on the value of phi(i)(i) the nodal coordinate must be between 1.0/phi(i)(i) and its negative.
575 XLe(i) = -bound;
576 XUe(i) = bound;
577 }
578 // The element bound vectors are now built for this element.
579 // Insert them into the global vectors, potentially overwriting
580 // the same dof contributions from other elements: no matter --
581 // the bounds are always -1.0 and 1.0.
582 XL.insert(XLe, dof_indices);
583 XU.insert(XUe, dof_indices);
584 }
585 }
586