MFEM  v4.2.0 Finite element discretization library
mesh-optimizer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // --------------------------------------------------
13 // Mesh Optimizer Miniapp: Optimize high-order meshes
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 as their
27 // coupling to Newton methods for solving minimization problems. Note that the
28 // utilized Newton methods are oriented towards avoiding invalid meshes with
29 // negative Jacobian determinants. Each Newton step requires the inversion of a
30 // Jacobian matrix, which is done through an inner linear solver.
31 //
32 // Compile with: make mesh-optimizer
33 //
34 // Sample runs:
36 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
38 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
40 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
41 //
43 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
44 // Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
45 // mesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
47 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
48 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
49 // Adapted discrete size+orientation (requires GSLIB):
50 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 8 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
51 // Adapted discrete aspect-ratio+orientation (requires GSLIB):
52 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 8 -ni 10 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd -ae 1
53 // Adapted discrete aspect ratio (3D):
54 // mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -ls 2 -li 100 -bnd -qt 1 -qo 8
55 //
57 // mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
58 // Adaptive limiting through the L-BFGS solver:
59 // mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
60 // Adaptive limiting through FD (requires GSLIB):
61 // * mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
62 //
64 // mesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
65 // Blade shape with FD-based solver:
66 // mesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8 -fd
68 // mesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8 -lc 5000
69 // ICF shape and equal size:
70 // mesh-optimizer -o 3 -rs 0 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
71 // ICF shape and initial size:
72 // mesh-optimizer -o 3 -rs 0 -mid 9 -tid 3 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8
73 // ICF shape:
74 // mesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8
75 // ICF limited shape:
76 // mesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -lc 10
77 // ICF combo shape + size (rings, slow convergence):
78 // mesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 1000 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb 1
79 // Mixed tet / cube / hex mesh with limiting:
80 // mesh-optimizer -m ../../data/fichera-mixed-p2.mesh -o 4 -rs 1 -mid 301 -tid 1 -fix-bnd -qo 6 -nor -lc 0.25
81 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
82 // * mesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -rs 0 -mid 303 -tid 1 -ni 20 -ls 2 -li 500 -fix-bnd
83 // 2D non-conforming shape and equal size:
84 // mesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
85
86
87 #include "mfem.hpp"
88 #include "../common/mfem-common.hpp"
89 #include <fstream>
90 #include <iostream>
91 #include "mesh-optimizer.hpp"
92
93 using namespace mfem;
94 using namespace std;
95
96 int main(int argc, char *argv[])
97 {
98  // 0. Set the method's default parameters.
99  const char *mesh_file = "icf.mesh";
100  int mesh_poly_deg = 1;
101  int rs_levels = 0;
102  double jitter = 0.0;
103  int metric_id = 1;
104  int target_id = 1;
105  double lim_const = 0.0;
109  int solver_type = 0;
110  int solver_iter = 10;
111  double solver_rtol = 1e-10;
112  int lin_solver = 2;
113  int max_lin_iter = 100;
114  bool move_bnd = true;
115  int combomet = 0;
116  bool normalization = false;
117  bool visualization = true;
118  int verbosity_level = 0;
119  bool fdscheme = false;
121  bool exactaction = false;
122
123  // 1. Parse command-line options.
124  OptionsParser args(argc, argv);
126  "Mesh file to use.");
128  "Polynomial degree of mesh finite element space.");
130  "Number of times to refine the mesh uniformly in serial.");
132  "Random perturbation scaling factor.");
134  "Mesh optimization metric:\n\t"
135  "1 : |T|^2 -- 2D shape\n\t"
136  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
137  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
138  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
139  "14: 0.5*(1-cos(theta_A - theta_W) -- 2D Sh+Sz+Alignment\n\t"
140  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
141  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
142  "55 : (tau-1)^2 -- 2D size\n\t"
143  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
144  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
145  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
146  "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
147  "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
148  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
149  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
150  "303: (|T|^2)/3*tau^(2/3)-1 -- 3D shape\n\t"
151  "315: (tau-1)^2 -- 3D size\n\t"
152  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D size\n\t"
153  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
154  "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling");
156  "Target (ideal element) type:\n\t"
157  "1: Ideal shape, unit size\n\t"
158  "2: Ideal shape, equal size\n\t"
159  "3: Ideal shape, initial size\n\t"
160  "4: Given full analytic Jacobian (in physical space)\n\t"
161  "5: Ideal shape, given size (in physical space)");
162  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
167  "1: Gauss-Lobatto\n\t"
168  "2: Gauss-Legendre\n\t"
169  "3: Closed uniform points");
171  "Order of the quadrature rule.");
173  " Type of solver: (default) 0: Newton, 1: LBFGS");
175  "Maximum number of Newton iterations.");
177  "Relative tolerance for the Newton solver.");
179  "Linear solver: 0 - l1-Jacobi, 1 - CG, 2 - MINRES.");
181  "Maximum number of iterations in the linear solve.");
183  "--fix-boundary",
184  "Enable motion along horizontal and vertical boundaries.");
186  "Combination of metrics options:"
187  "0: Use single metric\n\t"
188  "1: Shape + space-dependent size given analytically\n\t"
189  "2: Shape + adapted size given discretely; shared target");
191  "--no-normalization",
192  "Make all terms in the optimization functional unitless.");
194  "-no-fd", "--no-fd-approx",
195  "Enable finite difference based derivative computations.");
197  "-no-ex", "--no-exact-action",
198  "Enable exact action of TMOP_Integrator.");
200  "--no-visualization",
201  "Enable or disable GLVis visualization.");
203  "Set the verbosity level - 0, 1, or 2.");
205  "0 - Advection based (DEFAULT), 1 - GSLIB.");
206  args.Parse();
207  if (!args.Good())
208  {
209  args.PrintUsage(cout);
210  return 1;
211  }
212  args.PrintOptions(cout);
213
214  // 2. Initialize and refine the starting mesh.
215  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
216  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
217  const int dim = mesh->Dimension();
218  cout << "Mesh curvature: ";
219  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
220  else { cout << "(NONE)"; }
221  cout << endl;
222
223  // 3. Define a finite element space on the mesh. Here we use vector finite
224  // elements which are tensor products of quadratic finite elements. The
225  // number of components in the vector finite element space is specified by
226  // the last parameter of the FiniteElementSpace constructor.
228  if (mesh_poly_deg <= 0)
229  {
231  mesh_poly_deg = 2;
232  }
233  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
234  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim);
235
236  // 4. Make the mesh curved based on the above finite element space. This
237  // means that we define the mesh elements through a fespace-based
238  // transformation of the reference element.
239  mesh->SetNodalFESpace(fespace);
240
241  // 5. Set up an empty right-hand side vector b, which is equivalent to b=0.
242  Vector b(0);
243
244  // 6. Get the mesh nodes (vertices and other degrees of freedom in the finite
245  // element space) as a finite element grid function in fespace. Note that
246  // changing x automatically changes the shapes of the mesh elements.
247  GridFunction x(fespace);
248  mesh->SetNodalGridFunction(&x);
249
250  // 7. Define a vector representing the minimal local mesh size in the mesh
251  // nodes. We index the nodes using the scalar version of the degrees of
252  // freedom in fespace. Note: this is partition-dependent.
253  //
254  // In addition, compute average mesh size and total volume.
255  Vector h0(fespace->GetNDofs());
256  h0 = infinity();
257  double volume = 0.0;
258  Array<int> dofs;
259  for (int i = 0; i < mesh->GetNE(); i++)
260  {
261  // Get the local scalar element degrees of freedom in dofs.
262  fespace->GetElementDofs(i, dofs);
263  // Adjust the value of h0 in dofs based on the local mesh size.
264  const double hi = mesh->GetElementSize(i);
265  for (int j = 0; j < dofs.Size(); j++)
266  {
267  h0(dofs[j]) = min(h0(dofs[j]), hi);
268  }
269  volume += mesh->GetElementVolume(i);
270  }
271  const double small_phys_size = pow(volume, 1.0 / dim) / 100.0;
272
273  // 8. Add a random perturbation to the nodes in the interior of the domain.
274  // We define a random grid function of fespace and make sure that it is
275  // zero on the boundary and its values are locally of the order of h0.
276  // The latter is based on the DofToVDof() method which maps the scalar to
277  // the vector degrees of freedom in fespace.
278  GridFunction rdm(fespace);
279  rdm.Randomize();
280  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
281  rdm *= jitter;
282  // Scale the random values to be of order of the local mesh size.
283  for (int i = 0; i < fespace->GetNDofs(); i++)
284  {
285  for (int d = 0; d < dim; d++)
286  {
287  rdm(fespace->DofToVDof(i,d)) *= h0(i);
288  }
289  }
290  Array<int> vdofs;
291  for (int i = 0; i < fespace->GetNBE(); i++)
292  {
293  // Get the vector degrees of freedom in the boundary element.
294  fespace->GetBdrElementVDofs(i, vdofs);
295  // Set the boundary values to zero.
296  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
297  }
298  x -= rdm;
299  x.SetTrueVector();
300  x.SetFromTrueVector();
301
302  // 9. Save the starting (prior to the optimization) mesh to a file. This
303  // output can be viewed later using GLVis: "glvis -m perturbed.mesh".
304  {
305  ofstream mesh_ofs("perturbed.mesh");
306  mesh->Print(mesh_ofs);
307  }
308
309  // 10. Store the starting (prior to the optimization) positions.
310  GridFunction x0(fespace);
311  x0 = x;
312
313  // 11. Form the integrator that uses the chosen metric and target.
314  double tauval = -0.1;
315  TMOP_QualityMetric *metric = NULL;
316  switch (metric_id)
317  {
318  case 1: metric = new TMOP_Metric_001; break;
319  case 2: metric = new TMOP_Metric_002; break;
320  case 7: metric = new TMOP_Metric_007; break;
321  case 9: metric = new TMOP_Metric_009; break;
322  case 14: metric = new TMOP_Metric_SSA2D; break;
323  case 22: metric = new TMOP_Metric_022(tauval); break;
324  case 50: metric = new TMOP_Metric_050; break;
325  case 55: metric = new TMOP_Metric_055; break;
326  case 56: metric = new TMOP_Metric_056; break;
327  case 58: metric = new TMOP_Metric_058; break;
328  case 77: metric = new TMOP_Metric_077; break;
329  case 85: metric = new TMOP_Metric_085; break;
330  case 211: metric = new TMOP_Metric_211; break;
331  case 252: metric = new TMOP_Metric_252(tauval); break;
332  case 301: metric = new TMOP_Metric_301; break;
333  case 302: metric = new TMOP_Metric_302; break;
334  case 303: metric = new TMOP_Metric_303; break;
335  case 315: metric = new TMOP_Metric_315; break;
336  case 316: metric = new TMOP_Metric_316; break;
337  case 321: metric = new TMOP_Metric_321; break;
338  case 352: metric = new TMOP_Metric_352(tauval); break;
339  default: cout << "Unknown metric_id: " << metric_id << endl; return 3;
340  }
342  TargetConstructor *target_c = NULL;
344  H1_FECollection ind_fec(mesh_poly_deg, dim);
345  FiniteElementSpace ind_fes(mesh, &ind_fec);
346  FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
347  GridFunction size(&ind_fes), aspr(&ind_fes), disc(&ind_fes), ori(&ind_fes);
348  GridFunction aspr3d(&ind_fesv), size3d(&ind_fesv);
349  switch (target_id)
350  {
351  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
352  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
353  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
354  case 4: // Analytic
355  {
358  adapt_coeff = new HessianCoefficient(dim, metric_id);
360  target_c = tc;
361  break;
362  }
363  case 5: // Discrete size 2D
364  {
368  {
370  }
371  else
372  {
373 #ifdef MFEM_USE_GSLIB
375 #else
376  MFEM_ABORT("MFEM is not built with GSLIB.");
377 #endif
378  }
380  size.ProjectCoefficient(ind_coeff);
381  tc->SetSerialDiscreteTargetSize(size);
382  target_c = tc;
383  break;
384  }
385  case 6: // Discrete size + aspect ratio - 2D
386  {
387  GridFunction d_x(&ind_fes), d_y(&ind_fes);
388
392  disc.ProjectCoefficient(ind_coeff);
394  {
396  }
397  else
398  {
399 #ifdef MFEM_USE_GSLIB
401 #else
402  MFEM_ABORT("MFEM is not built with GSLIB.");
403 #endif
404  }
405
406  //Diffuse the interface
407  DiffuseField(disc,2);
408
409  //Get partials with respect to x and y of the grid function
410  disc.GetDerivative(1,0,d_x);
411  disc.GetDerivative(1,1,d_y);
412
413  //Compute the squared magnitude of the gradient
414  for (int i = 0; i < size.Size(); i++)
415  {
416  size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
417  }
418  const double max = size.Max();
419
420  for (int i = 0; i < d_x.Size(); i++)
421  {
422  d_x(i) = std::abs(d_x(i));
423  d_y(i) = std::abs(d_y(i));
424  }
425  const double eps = 0.01;
426  const double aspr_ratio = 20.0;
427  const double size_ratio = 40.0;
428
429  for (int i = 0; i < size.Size(); i++)
430  {
431  size(i) = (size(i)/max);
432  aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
433  aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
434  if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
435  if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
436  }
437  Vector vals;
438  const int NE = mesh->GetNE();
439  double volume = 0.0, volume_ind = 0.0;
440
441  for (int i = 0; i < NE; i++)
442  {
444  const IntegrationRule &ir =
445  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr->OrderJ());
446  size.GetValues(i, ir, vals);
447  for (int j = 0; j < ir.GetNPoints(); j++)
448  {
449  const IntegrationPoint &ip = ir.IntPoint(j);
450  Tr->SetIntPoint(&ip);
451  volume += ip.weight * Tr->Weight();
452  volume_ind += vals(j) * ip.weight * Tr->Weight();
453  }
454  }
455
456  const double avg_zone_size = volume / NE;
457
458  const double small_avg_ratio = (volume_ind + (volume - volume_ind) /
459  size_ratio) /
460  volume;
461
462  const double small_zone_size = small_avg_ratio * avg_zone_size;
463  const double big_zone_size = size_ratio * small_zone_size;
464
465  for (int i = 0; i < size.Size(); i++)
466  {
467  const double val = size(i);
468  const double a = (big_zone_size - small_zone_size) / small_zone_size;
469  size(i) = big_zone_size / (1.0+a*val);
470  }
471
472  DiffuseField(size, 2);
473  DiffuseField(aspr, 2);
474
475  tc->SetSerialDiscreteTargetSize(size);
477  target_c = tc;
478  break;
479  }
480  case 7: // Discrete aspect ratio 3D
481  {
485  {
487  }
488  else
489  {
490 #ifdef MFEM_USE_GSLIB
492 #else
493  MFEM_ABORT("MFEM is not built with GSLIB.");
494 #endif
495  }
497  aspr3d.ProjectCoefficient(fd_aspr3d);
498
500  target_c = tc;
501  break;
502  }
503  case 8: // shape/size + orientation 2D
504  {
508  {
510  }
511  else
512  {
513 #ifdef MFEM_USE_GSLIB
515 #else
516  MFEM_ABORT("MFEM is not built with GSLIB.");
517 #endif
518  }
519
520  if (metric_id == 14)
521  {
522  ConstantCoefficient ind_coeff(0.1*0.1);
523  size.ProjectCoefficient(ind_coeff);
524  tc->SetSerialDiscreteTargetSize(size);
525  }
526
527  if (metric_id == 85)
528  {
530  aspr.ProjectCoefficient(aspr_coeff);
531  DiffuseField(aspr,2);
533  }
534
536  ori.ProjectCoefficient(ori_coeff);
538  target_c = tc;
539  break;
540  }
541  default: cout << "Unknown target_id: " << target_id << endl; return 3;
542  }
543  if (target_c == NULL)
544  {
545  target_c = new TargetConstructor(target_t);
546  }
547  target_c->SetNodes(x0);
548  TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c);
549  if (fdscheme) { he_nlf_integ->EnableFiniteDifferences(x); }
550  he_nlf_integ->SetExactActionFlag(exactaction);
551
552  // Setup the quadrature rules for the TMOP integrator.
553  IntegrationRules *irules = NULL;
555  {
556  case 1: irules = &IntRulesLo; break;
557  case 2: irules = &IntRules; break;
558  case 3: irules = &IntRulesCU; break;
559  default: cout << "Unknown quad_type: " << quad_type << endl; return 3;
560  }
562  if (dim == 2)
563  {
564  cout << "Triangle quadrature points: "
567  << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
568  }
569  if (dim == 3)
570  {
571  cout << "Tetrahedron quadrature points: "
573  << "\nHexahedron quadrature points: "
575  << "\nPrism quadrature points: "
576  << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
577  }
578
579  if (normalization) { he_nlf_integ->EnableNormalization(x0); }
580
581  // Limit the node movement.
582  // The limiting distances can be given by a general function of space.
583  GridFunction dist(fespace);
584  dist = 1.0;
585  // The small_phys_size is relevant only with proper normalization.
586  if (normalization) { dist = small_phys_size; }
587  ConstantCoefficient lim_coeff(lim_const);
588  if (lim_const != 0.0) { he_nlf_integ->EnableLimiting(x0, dist, lim_coeff); }
589
591  GridFunction zeta_0(&ind_fes);
595  {
597  zeta_0.ProjectCoefficient(alim_coeff);
598
600  else if (adapt_eval == 1)
601  {
602 #ifdef MFEM_USE_GSLIB
604 #else
605  MFEM_ABORT("MFEM is not built with GSLIB support!");
606 #endif
607  }
608  else { MFEM_ABORT("Bad interpolation option."); }
609
611  if (visualization)
612  {
613  socketstream vis1;
614  common::VisualizeField(vis1, "localhost", 19916, zeta_0, "Zeta 0",
615  300, 600, 300, 300);
616  }
617  }
618
619  // 12. Setup the final NonlinearForm (which defines the integral of interest,
620  // its first and second derivatives). Here we can use a combination of
621  // metrics, i.e., optimize the sum of two integrals, where both are
622  // scaled by used-defined space-dependent weights. Note that there are no
623  // command-line options for the weights and the type of the second
624  // metric; one should update those in the code.
625  NonlinearForm a(fespace);
626  ConstantCoefficient *coeff1 = NULL;
627  TMOP_QualityMetric *metric2 = NULL;
628  TargetConstructor *target_c2 = NULL;
630
631  // Explicit combination of metrics.
632  if (combomet > 0)
633  {
634  // First metric.
635  coeff1 = new ConstantCoefficient(1.0);
636  he_nlf_integ->SetCoefficient(*coeff1);
637
638  // Second metric.
639  metric2 = new TMOP_Metric_077;
640  TMOP_Integrator *he_nlf_integ2 = NULL;
641  if (combomet == 1)
642  {
643  target_c2 = new TargetConstructor(
645  target_c2->SetVolumeScale(0.01);
646  target_c2->SetNodes(x0);
647  he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
648  he_nlf_integ2->SetCoefficient(coeff2);
649  }
650  else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c); }
652  if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
653  he_nlf_integ2->SetExactActionFlag(exactaction);
654
658  if (normalization) { combo->EnableNormalization(x0); }
659  if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
660
662  }
664
665  const double init_energy = a.GetGridFunctionEnergy(x);
666
667  // Visualize the starting mesh and metric values.
668  // Note that for combinations of metrics, this only shows the first metric.
669  if (visualization)
670  {
671  char title[] = "Initial metric values";
672  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
673  }
674
675  // 13. Fix all boundary nodes, or fix only a given component depending on the
676  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
677  // fixed x/y/z components of the node. Attribute 4 corresponds to an
678  // entirely fixed node. Other boundary attributes do not affect the node
679  // movement boundary conditions.
680  if (move_bnd == false)
681  {
682  Array<int> ess_bdr(mesh->bdr_attributes.Max());
683  ess_bdr = 1;
684  a.SetEssentialBC(ess_bdr);
685  }
686  else
687  {
688  int n = 0;
689  for (int i = 0; i < mesh->GetNBE(); i++)
690  {
691  const int nd = fespace->GetBE(i)->GetDof();
692  const int attr = mesh->GetBdrElement(i)->GetAttribute();
693  MFEM_VERIFY(!(dim == 2 && attr == 3),
694  "Boundary attribute 3 must be used only for 3D meshes. "
695  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
697  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
698  if (attr == 4) { n += nd * dim; }
699  }
700  Array<int> ess_vdofs(n), vdofs;
701  n = 0;
702  for (int i = 0; i < mesh->GetNBE(); i++)
703  {
704  const int nd = fespace->GetBE(i)->GetDof();
705  const int attr = mesh->GetBdrElement(i)->GetAttribute();
706  fespace->GetBdrElementVDofs(i, vdofs);
707  if (attr == 1) // Fix x components.
708  {
709  for (int j = 0; j < nd; j++)
710  { ess_vdofs[n++] = vdofs[j]; }
711  }
712  else if (attr == 2) // Fix y components.
713  {
714  for (int j = 0; j < nd; j++)
715  { ess_vdofs[n++] = vdofs[j+nd]; }
716  }
717  else if (attr == 3) // Fix z components.
718  {
719  for (int j = 0; j < nd; j++)
720  { ess_vdofs[n++] = vdofs[j+2*nd]; }
721  }
722  else if (attr == 4) // Fix all components.
723  {
724  for (int j = 0; j < vdofs.Size(); j++)
725  { ess_vdofs[n++] = vdofs[j]; }
726  }
727  }
728  a.SetEssentialVDofs(ess_vdofs);
729  }
730
731  // 14. As we use the Newton method to solve the resulting nonlinear system,
732  // here we setup the linear solver for the system's Jacobian.
733  Solver *S = NULL;
734  const double linsol_rtol = 1e-12;
735  if (lin_solver == 0)
736  {
737  S = new DSmoother(1, 1.0, max_lin_iter);
738  }
739  else if (lin_solver == 1)
740  {
741  CGSolver *cg = new CGSolver;
742  cg->SetMaxIter(max_lin_iter);
743  cg->SetRelTol(linsol_rtol);
744  cg->SetAbsTol(0.0);
745  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
746  S = cg;
747  }
748  else
749  {
750  MINRESSolver *minres = new MINRESSolver;
751  minres->SetMaxIter(max_lin_iter);
752  minres->SetRelTol(linsol_rtol);
753  minres->SetAbsTol(0.0);
754  minres->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
755  S = minres;
756  }
757
758  // Compute the minimum det(J) of the starting mesh.
759  tauval = infinity();
760  const int NE = mesh->GetNE();
761  for (int i = 0; i < NE; i++)
762  {
763  const IntegrationRule &ir =
766  for (int j = 0; j < ir.GetNPoints(); j++)
767  {
768  transf->SetIntPoint(&ir.IntPoint(j));
769  tauval = min(tauval, transf->Jacobian().Det());
770  }
771  }
772  cout << "Minimum det(J) of the original mesh is " << tauval << endl;
773  tauval -= 0.01 * h0.Min(); // Slightly below minJ0 to avoid div by 0.
774
775  // Perform the nonlinear optimization.
776  const IntegrationRule &ir =
778  TMOPNewtonSolver solver(ir, solver_type);
779  // Provide all integration rules in case of a mixed mesh.
781  if (solver_type == 0)
782  {
783  // Specify linear solver when we use a Newton-based solver.
784  solver.SetPreconditioner(*S);
785  }
786  solver.SetMaxIter(solver_iter);
787  solver.SetRelTol(solver_rtol);
788  solver.SetAbsTol(0.0);
789  solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
790  solver.SetOperator(a);
791  solver.Mult(b, x.GetTrueVector());
792  x.SetFromTrueVector();
793  if (solver.GetConverged() == false)
794  {
795  cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
796  }
797
798  // 15. Save the optimized mesh to a file. This output can be viewed later
799  // using GLVis: "glvis -m optimized.mesh".
800  {
801  ofstream mesh_ofs("optimized.mesh");
802  mesh_ofs.precision(14);
803  mesh->Print(mesh_ofs);
804  }
805
806  // 16. Compute the amount of energy decrease.
807  const double fin_energy = a.GetGridFunctionEnergy(x);
808  double metric_part = fin_energy;
809  if (lim_const > 0.0 || adapt_lim_const > 0.0)
810  {
811  lim_coeff.constant = 0.0;
812  coef_zeta.constant = 0.0;
813  metric_part = a.GetGridFunctionEnergy(x);
814  lim_coeff.constant = lim_const;
816  }
817  cout << "Initial strain energy: " << init_energy
818  << " = metrics: " << init_energy
819  << " + limiting term: " << 0.0 << endl;
820  cout << " Final strain energy: " << fin_energy
821  << " = metrics: " << metric_part
822  << " + limiting term: " << fin_energy - metric_part << endl;
823  cout << "The strain energy decreased by: " << setprecision(12)
824  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
825
826  // 17. Visualize the final mesh and metric values.
827  if (visualization)
828  {
829  char title[] = "Final metric values";
830  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
831  }
832
833  if (adapt_lim_const > 0.0 && visualization)
834  {
835  socketstream vis0;
836  common::VisualizeField(vis0, "localhost", 19916, zeta_0, "Xi 0",
837  600, 600, 300, 300);
838  }
839
840  // 18. Visualize the mesh displacement.
841  if (visualization)
842  {
843  x0 -= x;
844  osockstream sock(19916, "localhost");
845  sock << "solution\n";
846  mesh->Print(sock);
847  x0.Save(sock);
848  sock.send();
849  sock << "window_title 'Displacements'\n"
850  << "window_geometry "
851  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
852  << "keys jRmclA" << endl;
853  }
854
855  // 19. Free the used memory.
856  delete S;
857  delete target_c2;
858  delete metric2;
859  delete coeff1;
861  delete target_c;
863  delete metric;
864  delete fespace;
865  delete fec;
866  delete mesh;
867
868  return 0;
869 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
void discrete_aspr_3d(const Vector &x, Vector &v)
int Size() const
Return the logical size of the array.
Definition: array.hpp:124
Shifted barrier form of metric 56 (area, ideal barrier metric), 2D.
Definition: tmop.hpp:354
Shifted barrier form of 3D metric 16 (volume, ideal barrier metric), 3D.
Definition: tmop.hpp:472
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1234
Definition: solvers.hpp:258
Definition: ex25.cpp:144
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:397
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:2747
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:144
double discrete_size_2d(const Vector &x)
Shape &amp; volume, ideal barrier metric, 3D.
Definition: tmop.hpp:456
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:915
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:147
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:142
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
Definition: tmop.cpp:1156
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
Definition: tmop.cpp:2619
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:740
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:406
Definition: tmop.hpp:849
double Det() const
Definition: densemat.cpp:451
void DiffuseField(GridFunction &field, int smooth_steps)
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
Container class for integration rules.
Definition: intrules.hpp:309
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:390
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:668
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:166
Restriction of the node positions to certain regions.
Definition: tmop.cpp:1921
Volume metric, 3D.
Definition: tmop.hpp:422
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:303
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop.hpp:1032
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
Definition: tmop.cpp:2848
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
Definition: tmop_tools.cpp:651
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
Definition: tmop.hpp:1123
MINRES method.
Definition: solvers.hpp:368
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:725
int main(int argc, char *argv[])
Definition: ex1.cpp:66
virtual int OrderJ() const =0
Return the order of the elements of the Jacobian of the transformation.
Shape &amp; area, ideal barrier metric, 2D.
Definition: tmop.hpp:182
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:834
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Adds the limiting term to the first integrator. Disables it for the rest.
Definition: tmop.cpp:2775
Shifted barrier form of metric 2 (shape, ideal barrier metric), 2D.
Definition: tmop.hpp:214
virtual void SetSerialDiscreteTargetAspectRatio(const GridFunction &tspec_)
Definition: tmop.cpp:1147
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:285
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3417
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:312
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:136
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Limiting of the mesh displacements (general version).
Definition: tmop.cpp:1896
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:436
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Definition: tmop.hpp:1044
Area metric, 2D.
Definition: tmop.hpp:249
Volume, ideal barrier metric, 3D.
Definition: tmop.hpp:438
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:111
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:665
void SetMaxIter(int max_it)
Definition: solvers.hpp:98
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:1696
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:483
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
double discrete_ori_2d(const Vector &x)
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4126
int GetConverged() const
Definition: solvers.hpp:102
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
const Vector & GetTrueVector() const
Definition: gridfunc.hpp:124
A general vector function coefficient.
double GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Definition: mesh.cpp:74
void SetAbsTol(double atol)
Definition: solvers.hpp:97
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:201
void SetRelTol(double rtol)
Definition: solvers.hpp:96
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:233
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: tmop_tools.hpp:182
Untangling metric, 2D.
Definition: tmop.hpp:335
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:315
double weight_fun(const Vector &x)
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:266
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
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)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
Adds a new TMOP_Integrator to the combination.
Definition: tmop.hpp:1141
double material_indicator_2d(const Vector &x)
Shape &amp; area metric, 2D.
Definition: tmop.hpp:198
double a
Definition: lissajous.cpp:41
double discrete_aspr_2d(const Vector &x)
Class for integration point with weight.
Definition: intrules.hpp:25
double GetGridFunctionEnergy(const Vector &x) const
Compute the enery corresponding to the state x.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:336
int dim
Definition: ex24.cpp:53
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:1798
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:45
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:2252
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
A general function coefficient.
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: tmop_tools.hpp:157
Vector data type.
Definition: vector.hpp:51
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
Definition: fem_extras.cpp:94
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6603
Shape+Size+Orientation metric, 2D.
Definition: tmop.hpp:151
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:159
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:607
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Specify essential boundary conditions.
Base class for solvers.
Definition: operator.hpp:634
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:2047
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:179
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:4154
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:603
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, TMOPMatrixCoefficient *mspec)
Definition: tmop.cpp:941
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
Definition: tmop.cpp:1130
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:378
Shape &amp; orientation metric, 2D.
Definition: tmop.hpp:320
double GetElementVolume(int i)
Definition: mesh.cpp:101
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:823
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:374
Metric without a type, 2D.
Definition: tmop.hpp:75
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1552
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:145
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:884