MFEM  v3.3.2
Finite element discretization library
pmesh-optimizer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 //
12 // ---------------------------------------------------------------------
13 // Mesh Optimizer Miniapp: Optimize high-order meshes - Parallel Version
14 // ---------------------------------------------------------------------
15 //
16 // This miniapp performs mesh optimization using the Target-Matrix Optimization
17 // Paradigm (TMOP) by P.Knupp et al., and a global variational minimization
18 // approach. It minimizes the quantity sum_T int_T mu(J(x)), where T are the
19 // target (ideal) elements, J is the Jacobian of the transformation from the
20 // target to the physical element, and mu is the mesh quality metric. This
21 // metric can measure shape, size or alignment of the region around each
22 // quadrature point. The combination of targets & quality metrics is used to
23 // optimize the physical node positions, i.e., they must be as close as possible
24 // to the shape / size / alignment of their targets. This code also demonstrates
25 // a possible use of nonlinear operators (the class TMOP_QualityMetric, defining
26 // mu(J), and the class TMOP_Integrator, defining int mu(J)), as well
27 // as their coupling to Newton methods for solving minimization problems. Note
28 // that the utilized Newton methods are oriented towards avoiding invalid meshes
29 // with negative Jacobian determinants. Each Newton step requires the inversion
30 // of a Jacobian matrix, which is done through an inner linear solver.
31 //
32 // Compile with: make pmesh-optimizer
33 //
34 // Sample runs:
35 // Blade shape:
36 // mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
37 // ICF shape and equal size:
38 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
39 // ICF shape and initial size:
40 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 9 -tid 3 -ni 5000 -ls 2 -li 100 -bnd -qt 1 -qo 8
41 // ICF shape:
42 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8
43 // ICF limited shape:
44 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -lc 6.67
45 // ICF combo shape + size (rings, slow convergence):
46 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 1 -tid 2 -ni 1000 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb
47 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
48 // * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -rs 0 -mid 303 -tid 1 -ni 20 -ls 2 -li 500 -fix-bnd
49 
50 #include "mfem.hpp"
51 #include <fstream>
52 #include <iostream>
53 
54 using namespace mfem;
55 using namespace std;
56 
57 double weight_fun(const Vector &x);
58 
59 // Metric values are visualized by creating an L2 finite element functions and
60 // computing the metric values at the nodes.
61 void vis_metric(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc,
62  ParMesh &pmesh, char *title, int position)
63 {
65  ParFiniteElementSpace fes(&pmesh, &fec, 1);
66  ParGridFunction metric(&fes);
67  InterpolateTMOP_QualityMetric(qm, tc, pmesh, metric);
68  socketstream sock;
69  if (pmesh.GetMyRank() == 0)
70  {
71  sock.open("localhost", 19916);
72  sock << "solution\n";
73  }
74  pmesh.PrintAsOne(sock);
75  metric.SaveAsOne(sock);
76  if (pmesh.GetMyRank() == 0)
77  {
78  sock << "window_title '"<< title << "'\n"
79  << "window_geometry "
80  << position << " " << 0 << " " << 600 << " " << 600 << "\n"
81  << "keys jRmclA" << endl;
82  }
83 }
84 
85 class RelaxedNewtonSolver : public NewtonSolver
86 {
87 private:
88  // Quadrature points that are checked for negative Jacobians etc.
89  const IntegrationRule &ir;
90 
91 public:
92  RelaxedNewtonSolver(const IntegrationRule &irule) : ir(irule) { }
93 
94 #ifdef MFEM_USE_MPI
95  RelaxedNewtonSolver(const IntegrationRule &irule, MPI_Comm _comm)
96  : NewtonSolver(_comm), ir(irule) { }
97 #endif
98 
99  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
100 };
101 
102 double RelaxedNewtonSolver::ComputeScalingFactor(const Vector &x,
103  const Vector &b) const
104 {
105  const ParNonlinearForm *nlf = dynamic_cast<const ParNonlinearForm *>(oper);
106  MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
107  const bool have_b = (b.Size() == Height());
108  ParFiniteElementSpace *pfes = nlf->ParFESpace();
109 
110  const int NE = pfes->GetParMesh()->GetNE(), dim = pfes->GetFE(0)->GetDim(),
111  dof = pfes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
112  Array<int> xdofs(dof * dim);
113  DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
114  Vector posV(pos.Data(), dof * dim);
115 
116  Vector x_out(x.Size());
117  ParGridFunction x_out_gf(pfes);
118  bool x_out_ok = false;
119  const double energy_in = nlf->GetEnergy(x);
120  double scale = 1.0, energy_out;
121  double norm0 = Norm(r);
122 
123  // Decreases the scaling of the update until the new mesh is valid.
124  for (int i = 0; i < 12; i++)
125  {
126  add(x, -scale, c, x_out);
127  x_out_gf.Distribute(x_out);
128 
129  energy_out = nlf->GetEnergy(x_out_gf);
130  if (energy_out > 1.2*energy_in || isnan(energy_out) != 0)
131  {
132  if (print_level >= 0)
133  { cout << "Scale = " << scale << " Increasing energy." << endl; }
134  scale *= 0.5; continue;
135  }
136 
137  int jac_ok = 1;
138  for (int i = 0; i < NE; i++)
139  {
140  pfes->GetElementVDofs(i, xdofs);
141  x_out_gf.GetSubVector(xdofs, posV);
142  for (int j = 0; j < nsp; j++)
143  {
144  pfes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
145  MultAtB(pos, dshape, Jpr);
146  if (Jpr.Det() <= 0.0) { jac_ok = 0; goto break2; }
147  }
148  }
149  break2:
150  int jac_ok_all;
151  MPI_Allreduce(&jac_ok, &jac_ok_all, 1, MPI_INT, MPI_LAND,
152  pfes->GetComm());
153 
154  if (jac_ok_all == 0)
155  {
156  if (print_level >= 0)
157  { cout << "Scale = " << scale << " Neg det(J) found." << endl; }
158  scale *= 0.5; continue;
159  }
160 
161  oper->Mult(x_out, r);
162  if (have_b) { r -= b; }
163  double norm = Norm(r);
164 
165  if (norm > 1.2*norm0)
166  {
167  if (print_level >= 0)
168  { cout << "Scale = " << scale << " Norm increased." << endl; }
169  scale *= 0.5; continue;
170  }
171  else { x_out_ok = true; break; }
172  }
173 
174  if (print_level >= 0)
175  {
176  cout << "Energy decrease: "
177  << (energy_in - energy_out) / energy_in * 100.0
178  << "% with " << scale << " scaling." << endl;
179  }
180 
181  if (x_out_ok == false) { scale = 0.0; }
182 
183  return scale;
184 }
185 
186 // Allows negative Jacobians. Used in untangling metrics.
187 class DescentNewtonSolver : public NewtonSolver
188 {
189 private:
190  // Quadrature points that are checked for negative Jacobians etc.
191  const IntegrationRule &ir;
192 
193 public:
194  DescentNewtonSolver(const IntegrationRule &irule) : ir(irule) { }
195 
196 #ifdef MFEM_USE_MPI
197  DescentNewtonSolver(const IntegrationRule &irule, MPI_Comm _comm)
198  : NewtonSolver(_comm), ir(irule) { }
199 #endif
200 
201  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
202 };
203 
204 double DescentNewtonSolver::ComputeScalingFactor(const Vector &x,
205  const Vector &b) const
206 {
207  const ParNonlinearForm *nlf = dynamic_cast<const ParNonlinearForm *>(oper);
208  MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
209  ParFiniteElementSpace *pfes = nlf->ParFESpace();
210 
211  const int NE = pfes->GetParMesh()->GetNE(), dim = pfes->GetFE(0)->GetDim(),
212  dof = pfes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
213  Array<int> xdofs(dof * dim);
214  DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
215  Vector posV(pos.Data(), dof * dim);
216 
217  ParGridFunction x_gf(pfes);
218  x_gf.Distribute(x);
219 
220  double min_detJ = numeric_limits<double>::infinity();
221  for (int i = 0; i < NE; i++)
222  {
223  pfes->GetElementVDofs(i, xdofs);
224  x_gf.GetSubVector(xdofs, posV);
225  for (int j = 0; j < nsp; j++)
226  {
227  pfes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
228  MultAtB(pos, dshape, Jpr);
229  min_detJ = min(min_detJ, Jpr.Det());
230  }
231  }
232  double min_detJ_all;
233  MPI_Allreduce(&min_detJ, &min_detJ_all, 1, MPI_DOUBLE, MPI_MIN,
234  pfes->GetComm());
235  if (print_level >= 0)
236  { cout << "Minimum det(J) = " << min_detJ_all << endl; }
237 
238  Vector x_out(x.Size());
239  bool x_out_ok = false;
240  const double energy_in = nlf->GetEnergy(x_gf);
241  double scale = 1.0, energy_out;
242 
243  for (int i = 0; i < 7; i++)
244  {
245  add(x, -scale, c, x_out);
246 
247  energy_out = nlf->GetEnergy(x_out);
248  if (energy_out > energy_in || isnan(energy_out) != 0)
249  {
250  scale *= 0.5;
251  }
252  else { x_out_ok = true; break; }
253  }
254 
255  if (print_level >= 0)
256  {
257  cout << "Energy decrease: "
258  << (energy_in - energy_out) / energy_in * 100.0
259  << "% with " << scale << " scaling." << endl;
260  }
261 
262  if (x_out_ok == false) { return 0.0; }
263 
264  return scale;
265 }
266 
267 // Additional IntegrationRules that can be used with the --quad-type option.
270 
271 
272 int main (int argc, char *argv[])
273 {
274  // 0. Initialize MPI.
275  int num_procs, myid;
276  MPI_Init(&argc, &argv);
277  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
278  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
279 
280  // 1. Set the method's default parameters.
281  const char *mesh_file = "icf.mesh";
282  int mesh_poly_deg = 1;
283  int rs_levels = 0;
284  int rp_levels = 0;
285  double jitter = 0.0;
286  int metric_id = 1;
287  int target_id = 1;
288  double lim_const = 0.0;
289  int quad_type = 1;
290  int quad_order = 8;
291  int newton_iter = 10;
292  double newton_rtol = 1e-12;
293  int lin_solver = 2;
294  int max_lin_iter = 100;
295  bool move_bnd = true;
296  bool combomet = 0;
297  bool visualization = true;
298  int verbosity_level = 0;
299 
300  // 2. Parse command-line options.
301  OptionsParser args(argc, argv);
302  args.AddOption(&mesh_file, "-m", "--mesh",
303  "Mesh file to use.");
304  args.AddOption(&mesh_poly_deg, "-o", "--order",
305  "Polynomial degree of mesh finite element space.");
306  args.AddOption(&rs_levels, "-rs", "--refine-serial",
307  "Number of times to refine the mesh uniformly in serial.");
308  args.AddOption(&rp_levels, "-rp", "--refine-parallel",
309  "Number of times to refine the mesh uniformly in parallel.");
310  args.AddOption(&jitter, "-ji", "--jitter",
311  "Random perturbation scaling factor.");
312  args.AddOption(&metric_id, "-mid", "--metric-id",
313  "Mesh optimization metric:\n\t"
314  "1 : |T|^2 -- 2D shape\n\t"
315  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
316  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
317  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
318  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
319  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
320  "55 : (tau-1)^2 -- 2D size\n\t"
321  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
322  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
323  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
324  "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
325  "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
326  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
327  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
328  "303: (|T|^2)/3*tau^(2/3)-1 -- 3D shape\n\t"
329  "315: (tau-1)^2 -- 3D size\n\t"
330  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D size\n\t"
331  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
332  "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling");
333  args.AddOption(&target_id, "-tid", "--target-id",
334  "Target (ideal element) type:\n\t"
335  "1: Ideal shape, unit size\n\t"
336  "2: Ideal shape, equal size\n\t"
337  "3: Ideal shape, initial size");
338  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
339  args.AddOption(&quad_type, "-qt", "--quad-type",
340  "Quadrature rule type:\n\t"
341  "1: Gauss-Lobatto\n\t"
342  "2: Gauss-Legendre\n\t"
343  "3: Closed uniform points");
344  args.AddOption(&quad_order, "-qo", "--quad_order",
345  "Order of the quadrature rule.");
346  args.AddOption(&newton_iter, "-ni", "--newton-iters",
347  "Maximum number of Newton iterations.");
348  args.AddOption(&newton_rtol, "-rtol", "--newton-rel-tolerance",
349  "Relative tolerance for the Newton solver.");
350  args.AddOption(&lin_solver, "-ls", "--lin-solver",
351  "Linear solver: 0 - l1-Jacobi, 1 - CG, 2 - MINRES.");
352  args.AddOption(&max_lin_iter, "-li", "--lin-iter",
353  "Maximum number of iterations in the linear solve.");
354  args.AddOption(&move_bnd, "-bnd", "--move-boundary", "-fix-bnd",
355  "--fix-boundary",
356  "Enable motion along horizontal and vertical boundaries.");
357  args.AddOption(&combomet, "-cmb", "--combo-met", "-no-cmb", "--no-combo-met",
358  "Combination of metrics.");
359  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
360  "--no-visualization",
361  "Enable or disable GLVis visualization.");
362  args.AddOption(&verbosity_level, "-vl", "--verbosity-level",
363  "Set the verbosity level - 0, 1, or 2.");
364  args.Parse();
365  if (!args.Good())
366  {
367  if (myid == 0) { args.PrintUsage(cout); }
368  return 1;
369  }
370  if (myid == 0) { args.PrintOptions(cout); }
371 
372  // 3. Initialize and refine the starting mesh.
373  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
374  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
375  const int dim = mesh->Dimension();
376  if (myid == 0)
377  {
378  cout << "Mesh curvature: ";
379  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
380  else { cout << "(NONE)"; }
381  cout << endl;
382  }
383  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
384  delete mesh;
385  for (int lev = 0; lev < rp_levels; lev++) { pmesh->UniformRefinement(); }
386 
387  // 4. Define a finite element space on the mesh. Here we use vector finite
388  // elements which are tensor products of quadratic finite elements. The
389  // number of components in the vector finite element space is specified by
390  // the last parameter of the FiniteElementSpace constructor.
392  if (mesh_poly_deg <= 0)
393  {
394  fec = new QuadraticPosFECollection;
395  mesh_poly_deg = 2;
396  }
397  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
398  ParFiniteElementSpace *pfespace = new ParFiniteElementSpace(pmesh, fec, dim);
399 
400  // 5. Make the mesh curved based on the above finite element space. This
401  // means that we define the mesh elements through a fespace-based
402  // transformation of the reference element.
403  pmesh->SetNodalFESpace(pfespace);
404 
405  // 6. Set up an empty right-hand side vector b, which is equivalent to b=0.
406  Vector b(0);
407 
408  // 7. Get the mesh nodes (vertices and other degrees of freedom in the finite
409  // element space) as a finite element grid function in fespace. Note that
410  // changing x automatically changes the shapes of the mesh elements.
411  ParGridFunction x(pfespace);
412  pmesh->SetNodalGridFunction(&x);
413 
414  // 8. Define a vector representing the minimal local mesh size in the mesh
415  // nodes. We index the nodes using the scalar version of the degrees of
416  // freedom in pfespace.
417  Vector h0(pfespace->GetNDofs());
418  h0 = numeric_limits<double>::infinity();
419  Array<int> dofs;
420  for (int i = 0; i < pmesh->GetNE(); i++)
421  {
422  // Get the local scalar element degrees of freedom in dofs.
423  pfespace->GetElementDofs(i, dofs);
424  // Adjust the value of h0 in dofs based on the local mesh size.
425  for (int j = 0; j < dofs.Size(); j++)
426  {
427  h0(dofs[j]) = min(h0(dofs[j]), pmesh->GetElementSize(i));
428  }
429  }
430 
431  // 9. Add a random perturbation to the nodes in the interior of the domain.
432  // We define a random grid function of fespace and make sure that it is
433  // zero on the boundary and its values are locally of the order of h0.
434  // The latter is based on the DofToVDof() method which maps the scalar to
435  // the vector degrees of freedom in fespace.
436  ParGridFunction rdm(pfespace);
437  rdm.Randomize();
438  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
439  rdm *= jitter;
440  // Scale the random values to be of order of the local mesh size.
441  for (int i = 0; i < pfespace->GetNDofs(); i++)
442  {
443  for (int d = 0; d < dim; d++)
444  {
445  rdm(pfespace->DofToVDof(i,d)) *= h0(i);
446  }
447  }
448  Array<int> vdofs;
449  for (int i = 0; i < pfespace->GetNBE(); i++)
450  {
451  // Get the vector degrees of freedom in the boundary element.
452  pfespace->GetBdrElementVDofs(i, vdofs);
453  // Set the boundary values to zero.
454  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
455  }
456  // Average the perturbation of the overlapping nodes.
457  HypreParVector *trueF = rdm.ParallelAverage();
458  rdm = *trueF;
459  x -= rdm;
460  delete trueF;
461 
462  // 10. Save the starting (prior to the optimization) mesh to a file. This
463  // output can be viewed later using GLVis: "glvis -m perturbed -np
464  // num_mpi_tasks".
465  {
466  ostringstream mesh_name;
467  mesh_name << "perturbed." << setfill('0') << setw(6) << myid;
468  ofstream mesh_ofs(mesh_name.str().c_str());
469  mesh_ofs.precision(8);
470  pmesh->Print(mesh_ofs);
471  }
472 
473  // 11. Store the starting (prior to the optimization) positions.
474  ParGridFunction x0(pfespace);
475  x0 = x;
476 
477  // 12. Form the integrator that uses the chosen metric and target.
478  double tauval = -0.1;
479  TMOP_QualityMetric *metric = NULL;
480  switch (metric_id)
481  {
482  case 1: metric = new TMOP_Metric_001; break;
483  case 2: metric = new TMOP_Metric_002; break;
484  case 7: metric = new TMOP_Metric_007; break;
485  case 9: metric = new TMOP_Metric_009; break;
486  case 22: metric = new TMOP_Metric_022(tauval); break;
487  case 50: metric = new TMOP_Metric_050; break;
488  case 55: metric = new TMOP_Metric_055; break;
489  case 56: metric = new TMOP_Metric_056; break;
490  case 58: metric = new TMOP_Metric_058; break;
491  case 77: metric = new TMOP_Metric_077; break;
492  case 211: metric = new TMOP_Metric_211; break;
493  case 252: metric = new TMOP_Metric_252(tauval); break;
494  case 301: metric = new TMOP_Metric_301; break;
495  case 302: metric = new TMOP_Metric_302; break;
496  case 303: metric = new TMOP_Metric_303; break;
497  case 315: metric = new TMOP_Metric_315; break;
498  case 316: metric = new TMOP_Metric_316; break;
499  case 321: metric = new TMOP_Metric_321; break;
500  case 352: metric = new TMOP_Metric_352(tauval); break;
501  default:
502  if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
503  return 3;
504  }
506  switch (target_id)
507  {
508  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
509  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
510  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
511  default:
512  if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
513  return 3;
514  }
515  TargetConstructor *target_c;
516  target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
517  target_c->SetNodes(x);
518  TMOP_Integrator *he_nlf_integ;
519  he_nlf_integ = new TMOP_Integrator(metric, target_c);
520 
521  // 13. Setup the quadrature rule for the non-linear form integrator.
522  const IntegrationRule *ir = NULL;
523  const int geom_type = pfespace->GetFE(0)->GetGeomType();
524  switch (quad_type)
525  {
526  case 1: ir = &IntRulesLo.Get(geom_type, quad_order); break;
527  case 2: ir = &IntRules.Get(geom_type, quad_order); break;
528  case 3: ir = &IntRulesCU.Get(geom_type, quad_order); break;
529  default:
530  if (myid == 0) { cout << "Unknown quad_type: " << quad_type << endl; }
531  return 3;
532  }
533  if (myid == 0)
534  { cout << "Quadrature points per cell: " << ir->GetNPoints() << endl; }
535  he_nlf_integ->SetIntegrationRule(*ir);
536 
537  // 14. Limit the node movement.
538  ConstantCoefficient lim_coeff(lim_const);
539  if (lim_const != 0.0) { he_nlf_integ->EnableLimiting(x0, lim_coeff); }
540 
541  // 15. Setup the final NonlinearForm (which defines the integral of interest,
542  // its first and second derivatives). Here we can use a combination of
543  // metrics, i.e., optimize the sum of two integrals, where both are
544  // scaled by used-defined space-dependent weights. Note that there are
545  // no command-line options for the weights and the type of the second
546  // metric; one should update those in the code.
547  ParNonlinearForm a(pfespace);
548  Coefficient *coeff1 = NULL;
549  TMOP_QualityMetric *metric2 = NULL;
550  TargetConstructor *target_c2 = NULL;
552 
553  if (combomet)
554  {
555  // Weight of the original metric.
556  coeff1 = new ConstantCoefficient(1.0);
557  he_nlf_integ->SetCoefficient(*coeff1);
558  a.AddDomainIntegrator(he_nlf_integ);
559 
560  metric2 = new TMOP_Metric_077;
561  target_c2 = new TargetConstructor(target_t, MPI_COMM_WORLD);
562  target_c2->SetVolumeScale(0.01);
563  target_c2->SetNodes(x);
564  TMOP_Integrator *he_nlf_integ2;
565  he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
566  he_nlf_integ2->SetIntegrationRule(*ir);
567 
568  // Weight of metric2.
569  he_nlf_integ2->SetCoefficient(coeff2);
570  a.AddDomainIntegrator(he_nlf_integ2);
571  }
572  else { a.AddDomainIntegrator(he_nlf_integ); }
573  const double init_en = a.GetEnergy(x);
574  if (myid == 0) { cout << "Initial strain energy: " << init_en << endl; }
575 
576  // 16. Visualize the starting mesh and metric values.
577  if (visualization)
578  {
579  char title[] = "Initial metric values";
580  vis_metric(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
581  }
582 
583  // 17. Fix all boundary nodes, or fix only a given component depending on the
584  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
585  // fixed x/y/z components of the node. Attribute 4 corresponds to an
586  // entirely fixed node. Other boundary attributes do not affect the node
587  // movement boundary conditions.
588  if (move_bnd == false)
589  {
590  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
591  ess_bdr = 1;
592  a.SetEssentialBC(ess_bdr);
593  }
594  else
595  {
596  const int nd = pfespace->GetBE(0)->GetDof();
597  int n = 0;
598  for (int i = 0; i < pmesh->GetNBE(); i++)
599  {
600  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
601  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
602  if (attr == 4) { n += nd * dim; }
603  }
604  Array<int> ess_vdofs(n), vdofs;
605  n = 0;
606  for (int i = 0; i < pmesh->GetNBE(); i++)
607  {
608  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
609  pfespace->GetBdrElementVDofs(i, vdofs);
610  if (attr == 1) // Fix x components.
611  {
612  for (int j = 0; j < nd; j++)
613  { ess_vdofs[n++] = vdofs[j]; }
614  }
615  else if (attr == 2) // Fix y components.
616  {
617  for (int j = 0; j < nd; j++)
618  { ess_vdofs[n++] = vdofs[j+nd]; }
619  }
620  else if (attr == 3) // Fix z components.
621  {
622  for (int j = 0; j < nd; j++)
623  { ess_vdofs[n++] = vdofs[j+2*nd]; }
624  }
625  else if (attr == 4) // Fix all components.
626  {
627  for (int j = 0; j < vdofs.Size(); j++)
628  { ess_vdofs[n++] = vdofs[j]; }
629  }
630  }
631  a.SetEssentialVDofs(ess_vdofs);
632  }
633 
634  // 18. As we use the Newton method to solve the resulting nonlinear system,
635  // here we setup the linear solver for the system's Jacobian.
636  Solver *S = NULL;
637  const double linsol_rtol = 1e-12;
638  if (lin_solver == 0)
639  {
640  S = new DSmoother(1, 1.0, max_lin_iter);
641  }
642  else if (lin_solver == 1)
643  {
644  CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
645  cg->SetMaxIter(max_lin_iter);
646  cg->SetRelTol(linsol_rtol);
647  cg->SetAbsTol(0.0);
648  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
649  S = cg;
650  }
651  else
652  {
653  MINRESSolver *minres = new MINRESSolver(MPI_COMM_WORLD);
654  minres->SetMaxIter(max_lin_iter);
655  minres->SetRelTol(linsol_rtol);
656  minres->SetAbsTol(0.0);
657  minres->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
658  S = minres;
659  }
660 
661  // 19. Compute the minimum det(J) of the starting mesh.
662  tauval = numeric_limits<double>::infinity();
663  const int NE = pmesh->GetNE();
664  for (int i = 0; i < NE; i++)
665  {
667  for (int j = 0; j < ir->GetNPoints(); j++)
668  {
669  transf->SetIntPoint(&ir->IntPoint(j));
670  tauval = min(tauval, transf->Jacobian().Det());
671  }
672  }
673  double minJ0;
674  MPI_Allreduce(&tauval, &minJ0, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
675  tauval = minJ0;
676  if (myid == 0)
677  { cout << "Minimum det(J) of the original mesh is " << tauval << endl; }
678 
679  // 20. Finally, perform the nonlinear optimization.
680  NewtonSolver *newton = NULL;
681  if (tauval > 0.0)
682  {
683  tauval = 0.0;
684  newton = new RelaxedNewtonSolver(*ir, MPI_COMM_WORLD);
685  if (myid == 0)
686  { cout << "RelaxedNewtonSolver is used (as all det(J) > 0)." << endl; }
687  }
688  else
689  {
690  if ( (dim == 2 && metric_id != 22 && metric_id != 252) ||
691  (dim == 3 && metric_id != 352) )
692  {
693  if (myid == 0)
694  { cout << "The mesh is inverted. Use an untangling metric." << endl; }
695  return 3;
696  }
697  double h0min = h0.Min(), h0min_all;
698  MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
699  tauval -= 0.01 * h0min_all; // Slightly below minJ0 to avoid div by 0.
700  newton = new DescentNewtonSolver(*ir, MPI_COMM_WORLD);
701  if (myid == 0)
702  { cout << "DescentNewtonSolver is used (as some det(J) < 0)." << endl; }
703  }
704  newton->SetPreconditioner(*S);
705  newton->SetMaxIter(newton_iter);
706  newton->SetRelTol(newton_rtol);
707  newton->SetAbsTol(0.0);
708  newton->SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
709  newton->SetOperator(a);
710  Vector X(pfespace->TrueVSize());
711  pfespace->GetRestrictionMatrix()->Mult(x, X);
712  newton->Mult(b, X);
713  if (myid == 0 && newton->GetConverged() == false)
714  {
715  cout << "NewtonIteration: rtol = " << newton_rtol << " not achieved."
716  << endl;
717  }
718  pfespace->Dof_TrueDof_Matrix()->Mult(X, x);
719  delete newton;
720 
721  // 21. Save the optimized mesh to a file. This output can be viewed later
722  // using GLVis: "glvis -m optimized -np num_mpi_tasks".
723  {
724  ostringstream mesh_name;
725  mesh_name << "optimized." << setfill('0') << setw(6) << myid;
726  ofstream mesh_ofs(mesh_name.str().c_str());
727  mesh_ofs.precision(8);
728  pmesh->Print(mesh_ofs);
729  }
730 
731  // 22. Compute the amount of energy decrease.
732  const double fin_en = a.GetEnergy(x);
733  if (myid == 0)
734  {
735  cout << "Final strain energy : " << fin_en << endl;
736  cout << "The strain energy decreased by: " << setprecision(12)
737  << (init_en - fin_en) * 100.0 / init_en << " %." << endl;
738  }
739 
740  // 23. Visualize the final mesh and metric values.
741  if (visualization)
742  {
743  char title[] = "Final metric values";
744  vis_metric(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
745  }
746 
747  // 23. Visualize the mesh displacement.
748  if (visualization)
749  {
750  x0 -= x;
751  socketstream sock;
752  if (myid == 0)
753  {
754  sock.open("localhost", 19916);
755  sock << "solution\n";
756  }
757  pmesh->PrintAsOne(sock);
758  x0.SaveAsOne(sock);
759  if (myid == 0)
760  {
761  sock << "window_title 'Displacements'\n"
762  << "window_geometry "
763  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
764  << "keys jRmclA" << endl;
765  }
766  }
767 
768  // 24. Free the used memory.
769  delete S;
770  delete target_c2;
771  delete metric2;
772  delete coeff1;
773  delete target_c;
774  delete metric;
775  delete pfespace;
776  delete fec;
777  delete pmesh;
778 
779  MPI_Finalize();
780  return 0;
781 }
782 
783 // Defined with respect to the icf mesh.
784 double weight_fun(const Vector &x)
785 {
786  const double r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
787  const double den = 0.002;
788  double l2 = 0.2 + 0.5 * (std::tanh((r-0.16)/den) - std::tanh((r-0.17)/den)
789  + std::tanh((r-0.23)/den) - std::tanh((r-0.24)/den));
790  return l2;
791 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
int Size() const
Logical size of the array.
Definition: array.hpp:110
Shifted barrier form of metric 56 (area, ideal barrier metric), 2D.
Definition: tmop.hpp:265
Shifted barrier form of 3D metric 16 (volume, ideal barrier metric), 3D.
Definition: tmop.hpp:383
Conjugate gradient method.
Definition: solvers.hpp:111
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:116
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:161
ParFiniteElementSpace * ParFESpace() const
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:104
Shape & volume, ideal barrier metric, 3D.
Definition: tmop.hpp:367
Data type for scaled Jacobi-type smoother of sparse matrix.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
Subclass constant coefficient.
Definition: coefficient.hpp:57
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:617
void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric, const TargetConstructor &tc, const Mesh &mesh, GridFunction &metric_gf)
Interpolates the metric's values at the nodes of metric_gf.
Definition: tmop.cpp:986
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
Definition: fespace.cpp:133
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:317
double Det() const
Definition: densemat.cpp:435
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
void SaveAsOne(std::ostream &out=mfem::out)
Merge the local grid functions.
Definition: pgridfunc.cpp:376
Container class for integration rules.
Definition: intrules.hpp:279
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:301
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:473
Parallel non-linear operator on the true dofs.
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:92
Volume metric, 3D.
Definition: tmop.hpp:333
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:229
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:614
Abstract parallel finite element space.
Definition: pfespace.hpp:31
virtual const SparseMatrix * GetRestrictionMatrix()
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:250
STL namespace.
void vis_metric(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
MINRES method.
Definition: solvers.hpp:220
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:651
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0)
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:947
Shape & area, ideal barrier metric, 2D.
Definition: tmop.hpp:108
Shifted barrier form of metric 2 (shape, ideal barrier metric), 2D.
Definition: tmop.hpp:140
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:264
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:211
virtual double GetEnergy(const ParGridFunction &x) const
Compute the energy of a ParGridFunction.
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:24
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:225
double weight_fun(const Vector &x)
int dim
Definition: ex3.cpp:47
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:72
MPI_Comm GetComm() const
Definition: pfespace.hpp:166
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:190
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Definition: tmop.hpp:528
Area metric, 2D.
Definition: tmop.hpp:175
Volume, ideal barrier metric, 3D.
Definition: tmop.hpp:349
void PrintAsOne(std::ostream &out=mfem::out)
Definition: pmesh.cpp:3527
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:6718
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:67
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:470
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:198
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:108
Newton's method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:258
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:3399
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:423
int main(int argc, char *argv[])
int GetAttribute() const
Return element's attribute.
Definition: element.hpp:50
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3298
int GetConverged() const
Definition: solvers.hpp:67
int Dimension() const
Definition: mesh.hpp:641
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1241
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:119
double GetElementSize(int i, int type=0)
Definition: mesh.cpp:51
Wrapper for hypre's parallel vector class.
Definition: hypre.hpp:73
void SetAbsTol(double atol)
Definition: solvers.hpp:62
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:172
int GetMyRank() const
Definition: pmesh.hpp:119
void SetRelTol(double rtol)
Definition: solvers.hpp:61
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:159
Untangling metric, 2D.
Definition: tmop.hpp:246
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:122
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:192
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:475
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Definition: optparser.hpp:74
void EnableLimiting(const GridFunction &n0, Coefficient &w0)
Adds a limiting term to the integrator.
Definition: tmop.hpp:537
Shape & area metric, 2D.
Definition: tmop.hpp:124
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
Definition: pfespace.cpp:266
virtual void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
aka closed Newton-Cotes
Definition: intrules.hpp:267
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:255
void SetIntegrationRule(const IntegrationRule &irule)
Prescribe a fixed IntegrationRule to use.
Definition: nonlininteg.hpp:40
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i'th element.
Definition: fespace.cpp:1171
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
int open(const char hostname[], int port)
class for C-function coefficient
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:3637
Vector data type.
Definition: vector.hpp:41
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5401
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:93
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:414
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Base class for solvers.
Definition: operator.hpp:259
Class for parallel grid function.
Definition: pgridfunc.hpp:32
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i'th boundary element.
Definition: fespace.cpp:1406
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i'th boundary element.
Definition: fespace.cpp:139
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:3304
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:410
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Class for parallel meshes.
Definition: pmesh.hpp:29
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:343
void ParallelAverage(Vector &tv) const
Returns the vector averaged on the true dofs.
Definition: pgridfunc.cpp:133
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:676
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:285
Metric without a type, 2D.
Definition: tmop.hpp:76
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1230
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:195
bool Good() const
Definition: optparser.hpp:120
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:490