1 // MFEM Example 19 - Parallel Version
2 //
3 // Compile with: make ex19p
4 //
5 // Sample runs:
6 // mpirun -np 2 ex19p -m ../data/beam-quad.mesh
7 // mpirun -np 2 ex19p -m ../data/beam-tri.mesh
8 // mpirun -np 2 ex19p -m ../data/beam-hex.mesh
9 // mpirun -np 2 ex19p -m ../data/beam-tet.mesh
10 // mpirun -np 2 ex19p -m ../data/beam-wedge.mesh
11 // mpirun -np 2 ex19p -m ../data/beam-quad-amr.mesh
12 //
13 // Description: This examples solves a quasi-static incompressible nonlinear
14 // elasticity problem of the form 0 = H(x), where H is an
15 // incompressible hyperelastic model and x is a block state vector
16 // containing displacement and pressure variables. The geometry of
17 // the domain is assumed to be as follows:
18 //
19 // +---------------------+
20 // boundary --->| |<--- boundary
21 // attribute 1 | | attribute 2
22 // (fixed) +---------------------+ (fixed, nonzero)
23 //
24 // The example demonstrates the use of block nonlinear operators
25 // (the class RubberOperator defining H(x)) as well as a nonlinear
26 // Newton solver for the quasi-static problem. Each Newton step
27 // requires the inversion of a Jacobian matrix, which is done
28 // through a (preconditioned) inner solver. The specialized block
29 // preconditioner is implemented as a user-defined solver.
30 //
31 // We recommend viewing examples 2, 5, and 10 before viewing this
32 // example.
33
34 #include "mfem.hpp"
35 #include <memory>
36 #include <iostream>
37 #include <fstream>
38
39 using namespace std;
40 using namespace mfem;
41
42 class GeneralResidualMonitor : public IterativeSolverMonitor
43 {
44 public:
GeneralResidualMonitor(MPI_Comm comm,const std::string & prefix_,int print_lvl)45 GeneralResidualMonitor(MPI_Comm comm, const std::string& prefix_,
46 int print_lvl)
47 : prefix(prefix_)
48 {
49 #ifndef MFEM_USE_MPI
50 print_level = print_lvl;
51 #else
52 int rank;
53 MPI_Comm_rank(comm, &rank);
54 if (rank == 0)
55 {
56 print_level = print_lvl;
57 }
58 else
59 {
60 print_level = -1;
61 }
62 #endif
63 }
64
65 virtual void MonitorResidual(int it, double norm, const Vector &r, bool final);
66
67 private:
68 const std::string prefix;
69 int print_level;
70 mutable double norm0;
71 };
72
MonitorResidual(int it,double norm,const Vector & r,bool final)73 void GeneralResidualMonitor::MonitorResidual(int it, double norm,
74 const Vector &r, bool final)
75 {
76 if (print_level == 1 || (print_level == 3 && (final || it == 0)))
77 {
78 mfem::out << prefix << " iteration " << setw(2) << it
79 << " : ||r|| = " << norm;
80 if (it > 0)
81 {
82 mfem::out << ", ||r||/||r_0|| = " << norm/norm0;
83 }
84 else
85 {
86 norm0 = norm;
87 }
88 mfem::out << '\n';
89 }
90 }
91
92 // Custom block preconditioner for the Jacobian of the incompressible nonlinear
93 // elasticity operator. It has the form
94 //
95 // P^-1 = [ K^-1 0 ][ I -B^T ][ I 0 ]
96 // [ 0 I ][ 0 I ][ 0 -\gamma S^-1 ]
97 //
98 // where the original Jacobian has the form
99 //
100 // J = [ K B^T ]
101 // [ B 0 ]
102 //
103 // and K^-1 is an approximation of the inverse of the displacement part of the
104 // Jacobian and S^-1 is an approximation of the inverse of the Schur
105 // complement S = B K^-1 B^T. The Schur complement is approximated using
106 // a mass matrix of the pressure variables.
107 class JacobianPreconditioner : public Solver
108 {
109 protected:
110 // Finite element spaces for setting up preconditioner blocks
111 Array<ParFiniteElementSpace *> spaces;
112
113 // Offsets for extracting block vector segments
114 Array<int> &block_trueOffsets;
115
116 // Jacobian for block access
117 BlockOperator *jacobian;
118
119 // Scaling factor for the pressure mass matrix in the block preconditioner
120 double gamma;
121
122 // Objects for the block preconditioner application
123 Operator *pressure_mass;
124 Solver *mass_pcg;
125 Solver *mass_prec;
126 Solver *stiff_pcg;
127 Solver *stiff_prec;
128
129 public:
130 JacobianPreconditioner(Array<ParFiniteElementSpace *> &fes,
131 Operator &mass, Array<int> &offsets);
132
133 virtual void Mult(const Vector &k, Vector &y) const;
134 virtual void SetOperator(const Operator &op);
135
136 virtual ~JacobianPreconditioner();
137 };
138
139 // After spatial discretization, the rubber model can be written as:
140 // 0 = H(x)
141 // where x is the block vector representing the deformation and pressure and
142 // H(x) is the nonlinear incompressible neo-Hookean operator.
143 class RubberOperator : public Operator
144 {
145 protected:
146 // Finite element spaces
147 Array<ParFiniteElementSpace *> spaces;
148
149 // Block nonlinear form
150 ParBlockNonlinearForm *Hform;
151
152 // Pressure mass matrix for the preconditioner
153 Operator *pressure_mass;
154
155 // Newton solver for the hyperelastic operator
156 NewtonSolver newton_solver;
157 GeneralResidualMonitor newton_monitor;
158
159 // Solver for the Jacobian solve in the Newton method
160 Solver *j_solver;
161 GeneralResidualMonitor j_monitor;
162
163 // Preconditioner for the Jacobian
164 Solver *j_prec;
165
166 // Shear modulus coefficient
167 Coefficient μ
168
169 // Block offsets for variable access
170 Array<int> &block_trueOffsets;
171
172 public:
173 RubberOperator(Array<ParFiniteElementSpace *> &fes, Array<Array<int> *>&ess_bdr,
174 Array<int> &block_trueOffsets, double rel_tol, double abs_tol,
175 int iter, Coefficient &mu);
176
177 // Required to use the native newton solver
178 virtual Operator &GetGradient(const Vector &xp) const;
179 virtual void Mult(const Vector &k, Vector &y) const;
180
181 // Driver for the newton solver
182 void Solve(Vector &xp) const;
183
184 virtual ~RubberOperator();
185 };
186
187 // Visualization driver
188 void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
189 ParGridFunction *field, const char *field_name = NULL,
190 bool init_vis = false);
191
192 // Configuration definition functions
193 void ReferenceConfiguration(const Vector &x, Vector &y);
194 void InitialDeformation(const Vector &x, Vector &y);
195
196
main(int argc,char * argv[])197 int main(int argc, char *argv[])
198 {
199 #ifdef HYPRE_USING_CUDA
200 cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this example\n"
201 << "is NOT supported with the CUDA version of hypre.\n\n";
202 return 255;
203 #endif
204
205 // 1. Initialize MPI
206 MPI_Session mpi;
207 const int myid = mpi.WorldRank();
208
209 // 2. Parse command-line options
210 const char *mesh_file = "../data/beam-tet.mesh";
211 int ser_ref_levels = 0;
212 int par_ref_levels = 0;
213 int order = 2;
214 bool visualization = true;
215 double newton_rel_tol = 1e-4;
216 double newton_abs_tol = 1e-6;
217 int newton_iter = 500;
218 double mu = 1.0;
219
220 OptionsParser args(argc, argv);
221 args.AddOption(&mesh_file, "-m", "--mesh",
222 "Mesh file to use.");
223 args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
224 "Number of times to refine the mesh uniformly in serial.");
225 args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
226 "Number of times to refine the mesh uniformly in parallel.");
227 args.AddOption(&order, "-o", "--order",
228 "Order (degree) of the finite elements.");
229 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
230 "--no-visualization",
231 "Enable or disable GLVis visualization.");
232 args.AddOption(&newton_rel_tol, "-rel", "--relative-tolerance",
233 "Relative tolerance for the Newton solve.");
234 args.AddOption(&newton_abs_tol, "-abs", "--absolute-tolerance",
235 "Absolute tolerance for the Newton solve.");
236 args.AddOption(&newton_iter, "-it", "--newton-iterations",
237 "Maximum iterations for the Newton solve.");
238 args.AddOption(&mu, "-mu", "--shear-modulus",
239 "Shear modulus for the neo-Hookean material.");
240 args.Parse();
241 if (!args.Good())
242 {
243 if (myid == 0)
244 {
245 args.PrintUsage(cout);
246 }
247 return 1;
248 }
249 if (myid == 0)
250 {
251 args.PrintOptions(cout);
252 }
253
254 // 3. Read the (serial) mesh from the given mesh file on all processors. We
255 // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
256 // with the same code.
257 Mesh *mesh = new Mesh(mesh_file, 1, 1);
258 int dim = mesh->Dimension();
259
260 // 4. Refine the mesh in serial to increase the resolution. In this example
261 // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
262 // a command-line parameter.
263 for (int lev = 0; lev < ser_ref_levels; lev++)
264 {
265 mesh->UniformRefinement();
266 }
267
268 // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
269 // this mesh further in parallel to increase the resolution. Once the
270 // parallel mesh is defined, the serial mesh can be deleted.
271 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
272 delete mesh;
273 for (int lev = 0; lev < par_ref_levels; lev++)
274 {
275 pmesh->UniformRefinement();
276 }
277
278 // 6. Define the shear modulus for the incompressible Neo-Hookean material
279 ConstantCoefficient c_mu(mu);
280
281 // 7. Define the finite element spaces for displacement and pressure
282 // (Taylor-Hood elements). By default, the displacement (u/x) is a second
283 // order vector field, while the pressure (p) is a linear scalar function.
284 H1_FECollection quad_coll(order, dim);
285 H1_FECollection lin_coll(order-1, dim);
286
287 ParFiniteElementSpace R_space(pmesh, &quad_coll, dim, Ordering::byVDIM);
288 ParFiniteElementSpace W_space(pmesh, &lin_coll);
289
290 Array<ParFiniteElementSpace *> spaces(2);
291 spaces[0] = &R_space;
292 spaces[1] = &W_space;
293
294 HYPRE_BigInt glob_R_size = R_space.GlobalTrueVSize();
295 HYPRE_BigInt glob_W_size = W_space.GlobalTrueVSize();
296
297 // 8. Define the Dirichlet conditions (set to boundary attribute 1 and 2)
298 Array<Array<int> *> ess_bdr(2);
299
300 Array<int> ess_bdr_u(R_space.GetMesh()->bdr_attributes.Max());
301 Array<int> ess_bdr_p(W_space.GetMesh()->bdr_attributes.Max());
302
303 ess_bdr_p = 0;
304 ess_bdr_u = 0;
305 ess_bdr_u[0] = 1;
306 ess_bdr_u[1] = 1;
307
308 ess_bdr[0] = &ess_bdr_u;
309 ess_bdr[1] = &ess_bdr_p;
310
311 // 9. Print the mesh statistics
312 if (myid == 0)
313 {
314 std::cout << "***********************************************************\n";
315 std::cout << "dim(u) = " << glob_R_size << "\n";
316 std::cout << "dim(p) = " << glob_W_size << "\n";
317 std::cout << "dim(u+p) = " << glob_R_size + glob_W_size << "\n";
318 std::cout << "***********************************************************\n";
319 }
320
321 // 10. Define the block structure of the solution vector (u then p)
322 Array<int> block_trueOffsets(3);
323 block_trueOffsets[0] = 0;
324 block_trueOffsets[1] = R_space.TrueVSize();
325 block_trueOffsets[2] = W_space.TrueVSize();
326 block_trueOffsets.PartialSum();
327
328 BlockVector xp(block_trueOffsets);
329
330 // 11. Define grid functions for the current configuration, reference
331 // configuration, final deformation, and pressure
332 ParGridFunction x_gf(&R_space);
333 ParGridFunction x_ref(&R_space);
334 ParGridFunction x_def(&R_space);
335 ParGridFunction p_gf(&W_space);
336
337 VectorFunctionCoefficient deform(dim, InitialDeformation);
338 VectorFunctionCoefficient refconfig(dim, ReferenceConfiguration);
339
340 x_gf.ProjectCoefficient(deform);
341 x_ref.ProjectCoefficient(refconfig);
342 p_gf = 0.0;
343
344 // 12. Set up the block solution vectors
345 x_gf.GetTrueDofs(xp.GetBlock(0));
346 p_gf.GetTrueDofs(xp.GetBlock(1));
347
348 // 13. Initialize the incompressible neo-Hookean operator
349 RubberOperator oper(spaces, ess_bdr, block_trueOffsets,
350 newton_rel_tol, newton_abs_tol, newton_iter, c_mu);
351
352 // 14. Solve the Newton system
353 oper.Solve(xp);
354
355 // 15. Distribute the shared degrees of freedom
356 x_gf.Distribute(xp.GetBlock(0));
357 p_gf.Distribute(xp.GetBlock(1));
358
359 // 16. Compute the final deformation
360 subtract(x_gf, x_ref, x_def);
361
362 // 17. Visualize the results if requested
363 socketstream vis_u, vis_p;
364 if (visualization)
365 {
366 char vishost[] = "localhost";
367 int visport = 19916;
368 vis_u.open(vishost, visport);
369 vis_u.precision(8);
370 visualize(vis_u, pmesh, &x_gf, &x_def, "Deformation", true);
371 // Make sure all ranks have sent their 'u' solution before initiating
372 // another set of GLVis connections (one from each rank):
373 MPI_Barrier(pmesh->GetComm());
374 vis_p.open(vishost, visport);
375 vis_p.precision(8);
376 visualize(vis_p, pmesh, &x_gf, &p_gf, "Pressure", true);
377 }
378
379 // 18. Save the displaced mesh, the final deformation, and the pressure
380 {
381 GridFunction *nodes = &x_gf;
382 int owns_nodes = 0;
383 pmesh->SwapNodes(nodes, owns_nodes);
384
385 ostringstream mesh_name, pressure_name, deformation_name;
386 mesh_name << "mesh." << setfill('0') << setw(6) << myid;
387 pressure_name << "pressure." << setfill('0') << setw(6) << myid;
388 deformation_name << "deformation." << setfill('0') << setw(6) << myid;
389
390 ofstream mesh_ofs(mesh_name.str().c_str());
391 mesh_ofs.precision(8);
392 pmesh->Print(mesh_ofs);
393
394 ofstream pressure_ofs(pressure_name.str().c_str());
395 pressure_ofs.precision(8);
396 p_gf.Save(pressure_ofs);
397
398 ofstream deformation_ofs(deformation_name.str().c_str());
399 deformation_ofs.precision(8);
400 x_def.Save(deformation_ofs);
401 }
402
403 // 19. Free the used memory
404 delete pmesh;
405
406 return 0;
407 }
408
409
JacobianPreconditioner(Array<ParFiniteElementSpace * > & fes,Operator & mass,Array<int> & offsets)410 JacobianPreconditioner::JacobianPreconditioner(Array<ParFiniteElementSpace *>
411 &fes,
412 Operator &mass,
413 Array<int> &offsets)
414 : Solver(offsets[2]), block_trueOffsets(offsets), pressure_mass(&mass)
415 {
416 fes.Copy(spaces);
417
418 gamma = 0.00001;
419
420 // The mass matrix and preconditioner do not change every Newton cycle, so
421 // we only need to define them once
422 HypreBoomerAMG *mass_prec_amg = new HypreBoomerAMG();
423 mass_prec_amg->SetPrintLevel(0);
424
425 mass_prec = mass_prec_amg;
426
427 CGSolver *mass_pcg_iter = new CGSolver(spaces[0]->GetComm());
428 mass_pcg_iter->SetRelTol(1e-12);
429 mass_pcg_iter->SetAbsTol(1e-12);
430 mass_pcg_iter->SetMaxIter(200);
431 mass_pcg_iter->SetPrintLevel(0);
432 mass_pcg_iter->SetPreconditioner(*mass_prec);
433 mass_pcg_iter->SetOperator(*pressure_mass);
434 mass_pcg_iter->iterative_mode = false;
435
436 mass_pcg = mass_pcg_iter;
437
438 // The stiffness matrix does change every Newton cycle, so we will define it
439 // during SetOperator
440 stiff_pcg = NULL;
441 stiff_prec = NULL;
442 }
443
Mult(const Vector & k,Vector & y) const444 void JacobianPreconditioner::Mult(const Vector &k, Vector &y) const
445 {
446 // Extract the blocks from the input and output vectors
447 Vector disp_in;
448 disp_in.MakeRef(const_cast<Vector&>(k), block_trueOffsets[0],
449 block_trueOffsets[1]-block_trueOffsets[0]);
450 Vector pres_in;
451 pres_in.MakeRef(const_cast<Vector&>(k), block_trueOffsets[1],
452 block_trueOffsets[2]-block_trueOffsets[1]);
453
454 Vector disp_out;
455 disp_out.MakeRef(y, block_trueOffsets[0],
456 block_trueOffsets[1]-block_trueOffsets[0]);
457 Vector pres_out;
458 pres_out.MakeRef(y, block_trueOffsets[1],
459 block_trueOffsets[2]-block_trueOffsets[1]);
460
461 Vector temp(block_trueOffsets[1]-block_trueOffsets[0]);
462 Vector temp2(block_trueOffsets[1]-block_trueOffsets[0]);
463
464 // Perform the block elimination for the preconditioner
465 mass_pcg->Mult(pres_in, pres_out);
466 pres_out *= -gamma;
467
468 jacobian->GetBlock(0,1).Mult(pres_out, temp);
469 subtract(disp_in, temp, temp2);
470
471 stiff_pcg->Mult(temp2, disp_out);
472
473 disp_out.SyncAliasMemory(y);
474 pres_out.SyncAliasMemory(y);
475 }
476
SetOperator(const Operator & op)477 void JacobianPreconditioner::SetOperator(const Operator &op)
478 {
479 jacobian = (BlockOperator *) &op;
480
481 // Initialize the stiffness preconditioner and solver
482 if (stiff_prec == NULL)
483 {
484 HypreBoomerAMG *stiff_prec_amg = new HypreBoomerAMG();
485 stiff_prec_amg->SetPrintLevel(0);
486
487 if (!spaces[0]->GetParMesh()->Nonconforming())
488 {
489 #ifndef HYPRE_USING_CUDA
490 // Not available yet when hypre is built with CUDA
491 stiff_prec_amg->SetElasticityOptions(spaces[0]);
492 #endif
493 }
494
495 stiff_prec = stiff_prec_amg;
496
497 GMRESSolver *stiff_pcg_iter = new GMRESSolver(spaces[0]->GetComm());
498 stiff_pcg_iter->SetRelTol(1e-8);
499 stiff_pcg_iter->SetAbsTol(1e-8);
500 stiff_pcg_iter->SetMaxIter(200);
501 stiff_pcg_iter->SetPrintLevel(0);
502 stiff_pcg_iter->SetPreconditioner(*stiff_prec);
503 stiff_pcg_iter->iterative_mode = false;
504
505 stiff_pcg = stiff_pcg_iter;
506 }
507
508 // At each Newton cycle, compute the new stiffness AMG preconditioner by
509 // updating the iterative solver which, in turn, updates its preconditioner
510 stiff_pcg->SetOperator(jacobian->GetBlock(0,0));
511 }
512
~JacobianPreconditioner()513 JacobianPreconditioner::~JacobianPreconditioner()
514 {
515 delete mass_pcg;
516 delete mass_prec;
517 delete stiff_prec;
518 delete stiff_pcg;
519 }
520
521
RubberOperator(Array<ParFiniteElementSpace * > & fes,Array<Array<int> * > & ess_bdr,Array<int> & trueOffsets,double rel_tol,double abs_tol,int iter,Coefficient & c_mu)522 RubberOperator::RubberOperator(Array<ParFiniteElementSpace *> &fes,
523 Array<Array<int> *> &ess_bdr,
524 Array<int> &trueOffsets,
525 double rel_tol,
526 double abs_tol,
527 int iter,
528 Coefficient &c_mu)
529 : Operator(fes[0]->TrueVSize() + fes[1]->TrueVSize()),
530 newton_solver(fes[0]->GetComm()),
531 newton_monitor(fes[0]->GetComm(), "Newton", 1),
532 j_monitor(fes[0]->GetComm(), " GMRES", 3),
533 mu(c_mu), block_trueOffsets(trueOffsets)
534 {
535 Array<Vector *> rhs(2);
536 rhs = NULL; // Set all entries in the array
537
538 fes.Copy(spaces);
539
540 // Define the block nonlinear form
541 Hform = new ParBlockNonlinearForm(spaces);
542
543 // Add the incompressible neo-Hookean integrator
544 Hform->AddDomainIntegrator(new IncompressibleNeoHookeanIntegrator(mu));
545
546 // Set the essential boundary conditions
547 Hform->SetEssentialBC(ess_bdr, rhs);
548
549 // Compute the pressure mass stiffness matrix
550 ParBilinearForm *a = new ParBilinearForm(spaces[1]);
551 ConstantCoefficient one(1.0);
552 OperatorHandle mass(Operator::Hypre_ParCSR);
553 a->AddDomainIntegrator(new MassIntegrator(one));
554 a->Assemble();
555 a->Finalize();
556 a->ParallelAssemble(mass);
557 delete a;
558
559 mass.SetOperatorOwner(false);
560 pressure_mass = mass.Ptr();
561
562 // Initialize the Jacobian preconditioner
563 JacobianPreconditioner *jac_prec =
564 new JacobianPreconditioner(fes, *pressure_mass, block_trueOffsets);
565 j_prec = jac_prec;
566
567 // Set up the Jacobian solver
568 GMRESSolver *j_gmres = new GMRESSolver(spaces[0]->GetComm());
569 j_gmres->iterative_mode = false;
570 j_gmres->SetRelTol(1e-12);
571 j_gmres->SetAbsTol(1e-12);
572 j_gmres->SetMaxIter(300);
573 j_gmres->SetPrintLevel(-1);
574 j_gmres->SetMonitor(j_monitor);
575 j_gmres->SetPreconditioner(*j_prec);
576 j_solver = j_gmres;
577
578 // Set the newton solve parameters
579 newton_solver.iterative_mode = true;
580 newton_solver.SetSolver(*j_solver);
581 newton_solver.SetOperator(*this);
582 newton_solver.SetPrintLevel(-1);
583 newton_solver.SetMonitor(newton_monitor);
584 newton_solver.SetRelTol(rel_tol);
585 newton_solver.SetAbsTol(abs_tol);
586 newton_solver.SetMaxIter(iter);
587 }
588
589 // Solve the Newton system
Solve(Vector & xp) const590 void RubberOperator::Solve(Vector &xp) const
591 {
592 Vector zero;
593 newton_solver.Mult(zero, xp);
594 MFEM_VERIFY(newton_solver.GetConverged(),
595 "Newton Solver did not converge.");
596 }
597
598 // compute: y = H(x,p)
Mult(const Vector & k,Vector & y) const599 void RubberOperator::Mult(const Vector &k, Vector &y) const
600 {
601 Hform->Mult(k, y);
602 }
603
604 // Compute the Jacobian from the nonlinear form
GetGradient(const Vector & xp) const605 Operator &RubberOperator::GetGradient(const Vector &xp) const
606 {
607 return Hform->GetGradient(xp);
608 }
609
~RubberOperator()610 RubberOperator::~RubberOperator()
611 {
612 delete Hform;
613 delete pressure_mass;
614 delete j_solver;
615 delete j_prec;
616 }
617
618
619 // Inline visualization
visualize(ostream & out,ParMesh * mesh,ParGridFunction * deformed_nodes,ParGridFunction * field,const char * field_name,bool init_vis)620 void visualize(ostream &out, ParMesh *mesh, ParGridFunction *deformed_nodes,
621 ParGridFunction *field, const char *field_name, bool init_vis)
622 {
623 if (!out)
624 {
625 return;
626 }
627
628 GridFunction *nodes = deformed_nodes;
629 int owns_nodes = 0;
630
631 mesh->SwapNodes(nodes, owns_nodes);
632
633 out << "parallel " << mesh->GetNRanks() << " " << mesh->GetMyRank() << "\n";
634 out << "solution\n" << *mesh << *field;
635
636 mesh->SwapNodes(nodes, owns_nodes);
637
638 if (init_vis)
639 {
640 out << "window_size 800 800\n";
641 out << "window_title '" << field_name << "'\n";
642 if (mesh->SpaceDimension() == 2)
643 {
644 out << "view 0 0\n"; // view from top
645 out << "keys jlA\n"; // turn off perspective and light, +anti-aliasing
646 }
647 out << "keys cmA\n"; // show colorbar and mesh, +anti-aliasing
648 out << "autoscale value\n"; // update value-range; keep mesh-extents fixed
649 }
650 out << flush;
651 }
652
ReferenceConfiguration(const Vector & x,Vector & y)653 void ReferenceConfiguration(const Vector &x, Vector &y)
654 {
655 // Set the reference, stress free, configuration
656 y = x;
657 }
658
InitialDeformation(const Vector & x,Vector & y)659 void InitialDeformation(const Vector &x, Vector &y)
660 {
661 // Set the initial configuration. Having this different from the reference
662 // configuration can help convergence
663 y = x;
664 y[1] = x[1] + 0.25*x[0];
665 }
666