MFEM  v4.5.1 Finite element discretization library
mesh-optimizer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 -bnd -qt 1 -qo 8
38 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
40 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
41 //
43 // mesh-optimizer -m square01.mesh -o 2 -tid 9 -ni 50 -li 20 -hmid 55 -mid 7 -hr
44 // mesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
45 // mesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
46 //
48 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
49 // (requires GSLIB):
50 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor -mno 1 -ae 1
51 // Adapted discrete size 3D with PA:
52 // mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa
53 // Adapted discrete size 3D with PA on device (requires CUDA):
54 // * mesh-optimizer -m cube.mesh -o 3 -rs 3 -mid 321 -tid 5 -ls 3 -nor -lc 0.1 -pa -d cuda
55 // Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
56 // mesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
58 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
59 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
60 // Adapted discrete size+orientation (requires GSLIB):
61 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -ae 1 -nor
62 // Adapted discrete aspect-ratio+orientation (requires GSLIB):
63 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 8 -ni 10 -bnd -qt 1 -qo 8 -fd -ae 1
64 // Adapted discrete aspect ratio (3D):
65 // mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -bnd -qt 1 -qo 8
66 //
68 // mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
69 // Adaptive limiting through the L-BFGS solver:
70 // mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
71 // Adaptive limiting through FD (requires GSLIB):
72 // * mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
73 //
75 // mesh-optimizer -m square01.mesh -o 3 -rs 1 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 5e4 -rtol 1e-5
76 // mesh-optimizer -m square01-tri.mesh -o 3 -rs 0 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 1e4 -rtol 1e-5
77 // Surface fitting with weight adaptation and termination based on fitting error
78 // mesh-optimizer -m square01.mesh -o 2 -rs 1 -mid 2 -tid 1 -ni 100 -vl 2 -sfc 10 -rtol 1e-20 -st 0 -sfa -sft 1e-5
79 //
81 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
82 // (requires CUDA):
83 // * mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8 -d cuda
84 // Blade shape with FD-based solver:
85 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
87 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
88 // ICF shape and equal size:
89 // mesh-optimizer -o 3 -mid 9 -tid 2 -ni 25 -ls 3 -art 2 -qo 5
90 // ICF shape and initial size:
91 // mesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
92 // ICF shape:
93 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
94 // ICF limited shape:
95 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8 -lc 10
96 // ICF combo shape + size (rings, slow convergence):
97 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 1000 -bnd -qt 1 -qo 8 -cmb 1
98 // Mixed tet / cube / hex mesh with limiting:
99 // mesh-optimizer -m ../../data/fichera-mixed-p2.mesh -o 4 -rs 1 -mid 301 -tid 1 -fix-bnd -qo 6 -nor -lc 0.25
100 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
101 // * mesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -mid 303 -tid 1 -ni 20 -li 500 -fix-bnd
102 // 2D non-conforming shape and equal size:
103 // mesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -bnd -qt 1 -qo 8
104 //
105 // 2D untangling:
106 // mesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
107 // 2D untangling with shifted barrier metric:
108 // mesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
109 // 3D untangling (the mesh is in the mfem/data GitHub repository):
110 // * mesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
111
112 #include "mfem.hpp"
113 #include "../common/mfem-common.hpp"
114 #include <fstream>
115 #include <iostream>
116 #include "mesh-optimizer.hpp"
117
118 using namespace mfem;
119 using namespace std;
120
121 int main(int argc, char *argv[])
122 {
123  // 0. Set the method's default parameters.
124  const char *mesh_file = "icf.mesh";
125  int mesh_poly_deg = 1;
126  int rs_levels = 0;
127  double jitter = 0.0;
128  int metric_id = 1;
129  int target_id = 1;
130  double lim_const = 0.0;
132  double surface_fit_const = 0.0;
135  int solver_type = 0;
136  int solver_iter = 20;
137  double solver_rtol = 1e-10;
138  int solver_art_type = 0;
139  int lin_solver = 2;
140  int max_lin_iter = 100;
141  bool move_bnd = true;
142  int combomet = 0;
144  int h_metric_id = -1;
145  bool normalization = false;
146  bool visualization = true;
147  int verbosity_level = 0;
148  bool fdscheme = false;
150  bool exactaction = false;
151  const char *devopt = "cpu";
152  bool pa = false;
153  int n_hr_iter = 5;
154  int n_h_iter = 1;
156  double surface_fit_threshold = -10;
157  int mesh_node_ordering = 0;
158  int barrier_type = 0;
159  int worst_case_type = 0;
160
161  // 1. Parse command-line options.
162  OptionsParser args(argc, argv);
164  "Mesh file to use.");
166  "Polynomial degree of mesh finite element space.");
168  "Number of times to refine the mesh uniformly in serial.");
170  "Random perturbation scaling factor.");
172  "Mesh optimization metric:\n\t"
173  "T-metrics\n\t"
174  "1 : |T|^2 -- 2D no type\n\t"
175  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
176  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
177  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
178  "14 : |T-I|^2 -- 2D shape+size+orientation\n\t"
179  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
180  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
181  "55 : (tau-1)^2 -- 2D size\n\t"
182  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
183  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
184  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
185  "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t"
186  "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t"
187  "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t"
188  // "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
189  // "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
190  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
191  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
192  "303: (|T|^2)/3/tau^(2/3)-1 -- 3D shape\n\t"
193  "304: (|T|^3)/3^{3/2}/tau-1 -- 3D shape\n\t"
194  //"311: (tau-1)^2-tau+sqrt(tau^2+eps)-- 3D untangling\n\t"
195  "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t"
196  "315: (tau-1)^2 -- 3D no type\n\t"
197  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D no type\n\t"
198  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
199  "322: |T-adjT^-t|^2 -- 3D shape+size\n\t"
200  "323: |J|^3-3sqrt(3)ln(det(J))-3sqrt(3) -- 3D shape+size\n\t"
201  "328: (1-gamma) mu_301 + gamma mu_316 -- 3D shape+size\n\t"
202  "332: (1-gamma) mu_302 + gamma mu_315 -- 3D shape+size\n\t"
203  "333: (1-gamma) mu_302 + gamma mu_316 -- 3D shape+size\n\t"
204  "334: (1-gamma) mu_303 + gamma mu_316 -- 3D shape+size\n\t"
205  "347: (1-gamma) mu_304 + gamma mu_316 -- 3D shape+size\n\t"
206  // "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling\n\t"
207  "360: (|T|^3)/3^{3/2}-tau -- 3D shape\n\t"
208  "A-metrics\n\t"
209  "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t"
210  "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t"
211  "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t"
212  "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t"
213  );
215  "Target (ideal element) type:\n\t"
216  "1: Ideal shape, unit size\n\t"
217  "2: Ideal shape, equal size\n\t"
218  "3: Ideal shape, initial size\n\t"
219  "4: Given full analytic Jacobian (in physical space)\n\t"
220  "5: Ideal shape, given size (in physical space)");
221  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
225  "Surface preservation constant.");
228  "1: Gauss-Lobatto\n\t"
229  "2: Gauss-Legendre\n\t"
230  "3: Closed uniform points");
232  "Order of the quadrature rule.");
234  " Type of solver: (default) 0: Newton, 1: LBFGS");
236  "Maximum number of Newton iterations.");
238  "Relative tolerance for the Newton solver.");
240  "Type of adaptive relative linear solver tolerance:\n\t"
241  "0: None (default)\n\t"
242  "1: Eisenstat-Walker type 1\n\t"
243  "2: Eisenstat-Walker type 2");
245  "Linear solver:\n\t"
246  "0: l1-Jacobi\n\t"
247  "1: CG\n\t"
248  "2: MINRES\n\t"
249  "3: MINRES + Jacobi preconditioner\n\t"
250  "4: MINRES + l1-Jacobi preconditioner");
252  "Maximum number of iterations in the linear solve.");
254  "--fix-boundary",
255  "Enable motion along horizontal and vertical boundaries.");
257  "Combination of metrics options:\n\t"
258  "0: Use single metric\n\t"
259  "1: Shape + space-dependent size given analytically\n\t"
260  "2: Shape + adapted size given discretely; shared target");
265  "Same options as metric_id. Used to determine refinement"
266  " type for each element if h-adaptivity is enabled.");
268  "--no-normalization",
269  "Make all terms in the optimization functional unitless.");
271  "-no-fd", "--no-fd-approx",
272  "Enable finite difference based derivative computations.");
274  "-no-ex", "--no-exact-action",
275  "Enable exact action of TMOP_Integrator.");
277  "--no-visualization",
278  "Enable or disable GLVis visualization.");
280  "Set the verbosity level - 0, 1, or 2.");
282  "0 - Advection based (DEFAULT), 1 - GSLIB.");
284  "Device configuration string, see Device::Configure().");
286  "--no-partial-assembly", "Enable Partial Assembly.");
291  "iteration.");
294  "Enable or disable adaptive surface fitting.");
296  "Set threshold for surface fitting. TMOP solver will"
297  "terminate when max surface fitting error is below this limit");
299  "Ordering of mesh nodes."
300  "0 (default): byNodes, 1: byVDIM");
302  "0 - None,"
303  "1 - Shifted Barrier,"
304  "2 - Pseudo Barrier.");
306  "0 - None,"
307  "1 - Beta,"
308  "2 - PMean.");
309  args.Parse();
310  if (!args.Good())
311  {
312  args.PrintUsage(cout);
313  return 1;
314  }
315  args.PrintOptions(cout);
316
317  if (h_metric_id < 0) { h_metric_id = metric_id; }
318
320  {
321  MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
322  " supported on cpus.");
323  }
324  Device device(devopt);
325  device.Print();
326
327  // 2. Initialize and refine the starting mesh.
328  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
329  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
330  const int dim = mesh->Dimension();
331
332  if (hradaptivity) { mesh->EnsureNCMesh(); }
333
334  // 3. Define a finite element space on the mesh-> Here we use vector finite
335  // elements which are tensor products of quadratic finite elements. The
336  // number of components in the vector finite element space is specified by
337  // the last parameter of the FiniteElementSpace constructor.
339  if (mesh_poly_deg <= 0)
340  {
342  mesh_poly_deg = 2;
343  }
344  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
345  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim,
346  mesh_node_ordering);
347
348  // 4. Make the mesh curved based on the above finite element space. This
349  // means that we define the mesh elements through a fespace-based
350  // transformation of the reference element.
351  mesh->SetNodalFESpace(fespace);
352
353  // 5. Set up an empty right-hand side vector b, which is equivalent to b=0.
354  Vector b(0);
355
356  // 6. Get the mesh nodes (vertices and other degrees of freedom in the finite
357  // element space) as a finite element grid function in fespace. Note that
358  // changing x automatically changes the shapes of the mesh elements.
359  GridFunction x(fespace);
360  mesh->SetNodalGridFunction(&x);
361
362  // 7. Define a vector representing the minimal local mesh size in the mesh
363  // nodes. We index the nodes using the scalar version of the degrees of
364  // freedom in fespace. Note: this is partition-dependent.
365  //
366  // In addition, compute average mesh size and total volume.
367  Vector h0(fespace->GetNDofs());
368  h0 = infinity();
369  double mesh_volume = 0.0;
370  Array<int> dofs;
371  for (int i = 0; i < mesh->GetNE(); i++)
372  {
373  // Get the local scalar element degrees of freedom in dofs.
374  fespace->GetElementDofs(i, dofs);
375  // Adjust the value of h0 in dofs based on the local mesh size.
376  const double hi = mesh->GetElementSize(i);
377  for (int j = 0; j < dofs.Size(); j++)
378  {
379  h0(dofs[j]) = min(h0(dofs[j]), hi);
380  }
381  mesh_volume += mesh->GetElementVolume(i);
382  }
383  const double small_phys_size = pow(mesh_volume, 1.0 / dim) / 100.0;
384
385  // 8. Add a random perturbation to the nodes in the interior of the domain.
386  // We define a random grid function of fespace and make sure that it is
387  // zero on the boundary and its values are locally of the order of h0.
388  // The latter is based on the DofToVDof() method which maps the scalar to
389  // the vector degrees of freedom in fespace.
390  GridFunction rdm(fespace);
391  rdm.Randomize();
392  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
393  rdm *= jitter;
395  // Scale the random values to be of order of the local mesh size.
396  for (int i = 0; i < fespace->GetNDofs(); i++)
397  {
398  for (int d = 0; d < dim; d++)
399  {
400  rdm(fespace->DofToVDof(i,d)) *= h0(i);
401  }
402  }
403  Array<int> vdofs;
404  for (int i = 0; i < fespace->GetNBE(); i++)
405  {
406  // Get the vector degrees of freedom in the boundary element.
407  fespace->GetBdrElementVDofs(i, vdofs);
408  // Set the boundary values to zero.
409  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
410  }
411  x -= rdm;
412  x.SetTrueVector();
413  x.SetFromTrueVector();
414
415  // 9. Save the starting (prior to the optimization) mesh to a file. This
416  // output can be viewed later using GLVis: "glvis -m perturbed.mesh".
417  {
418  ofstream mesh_ofs("perturbed.mesh");
419  mesh->Print(mesh_ofs);
420  }
421
422  // 10. Store the starting (prior to the optimization) positions.
423  GridFunction x0(fespace);
424  x0 = x;
425
426  // 11. Form the integrator that uses the chosen metric and target.
427  double min_detJ = -0.1;
428  TMOP_QualityMetric *metric = NULL;
429  switch (metric_id)
430  {
431  // T-metrics
432  case 1: metric = new TMOP_Metric_001; break;
433  case 2: metric = new TMOP_Metric_002; break;
434  case 4: metric = new TMOP_Metric_004; break;
435  case 7: metric = new TMOP_Metric_007; break;
436  case 9: metric = new TMOP_Metric_009; break;
437  case 14: metric = new TMOP_Metric_014; break;
438  case 22: metric = new TMOP_Metric_022(min_detJ); break;
439  case 50: metric = new TMOP_Metric_050; break;
440  case 55: metric = new TMOP_Metric_055; break;
441  case 56: metric = new TMOP_Metric_056; break;
442  case 58: metric = new TMOP_Metric_058; break;
443  case 66: metric = new TMOP_Metric_066(0.5); break;
444  case 77: metric = new TMOP_Metric_077; break;
445  case 80: metric = new TMOP_Metric_080(0.5); break;
446  case 85: metric = new TMOP_Metric_085; break;
447  case 98: metric = new TMOP_Metric_098; break;
448  // case 211: metric = new TMOP_Metric_211; break;
449  // case 252: metric = new TMOP_Metric_252(min_detJ); break;
450  case 301: metric = new TMOP_Metric_301; break;
451  case 302: metric = new TMOP_Metric_302; break;
452  case 303: metric = new TMOP_Metric_303; break;
453  case 304: metric = new TMOP_Metric_304; break;
454  // case 311: metric = new TMOP_Metric_311; break;
455  case 313: metric = new TMOP_Metric_313(min_detJ); break;
456  case 315: metric = new TMOP_Metric_315; break;
457  case 316: metric = new TMOP_Metric_316; break;
458  case 321: metric = new TMOP_Metric_321; break;
459  case 322: metric = new TMOP_Metric_322; break;
460  case 323: metric = new TMOP_Metric_323; break;
461  case 328: metric = new TMOP_Metric_328(0.5); break;
462  case 332: metric = new TMOP_Metric_332(0.5); break;
463  case 333: metric = new TMOP_Metric_333(0.5); break;
464  case 334: metric = new TMOP_Metric_334(0.5); break;
465  case 347: metric = new TMOP_Metric_347(0.5); break;
466  // case 352: metric = new TMOP_Metric_352(min_detJ); break;
467  case 360: metric = new TMOP_Metric_360; break;
468  // A-metrics
469  case 11: metric = new TMOP_AMetric_011; break;
470  case 36: metric = new TMOP_AMetric_036; break;
471  case 107: metric = new TMOP_AMetric_107a; break;
472  case 126: metric = new TMOP_AMetric_126(0.9); break;
473  default:
474  cout << "Unknown metric_id: " << metric_id << endl;
475  return 3;
476  }
477  TMOP_QualityMetric *h_metric = NULL;
479  {
480  switch (h_metric_id)
481  {
482  case 1: h_metric = new TMOP_Metric_001; break;
483  case 2: h_metric = new TMOP_Metric_002; break;
484  case 7: h_metric = new TMOP_Metric_007; break;
485  case 9: h_metric = new TMOP_Metric_009; break;
486  case 55: h_metric = new TMOP_Metric_055; break;
487  case 56: h_metric = new TMOP_Metric_056; break;
488  case 58: h_metric = new TMOP_Metric_058; break;
489  case 77: h_metric = new TMOP_Metric_077; break;
490  case 315: h_metric = new TMOP_Metric_315; break;
491  case 316: h_metric = new TMOP_Metric_316; break;
492  case 321: h_metric = new TMOP_Metric_321; break;
493  default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
494  endl;
495  return 3;
496  }
497  }
498
500  switch (barrier_type)
501  {
503  break;
505  break;
507  break;
508  default: cout << "barrier_type not supported: " << barrier_type << endl;
509  return 3;
510  }
511
513  switch (worst_case_type)
514  {
516  break;
518  break;
520  break;
521  default: cout << "worst_case_type not supported: " << worst_case_type << endl;
522  return 3;
523  }
524
525  TMOP_QualityMetric *untangler_metric = NULL;
526  if (barrier_type > 0 || worst_case_type > 0)
527  {
528  if (barrier_type > 0)
529  {
530  MFEM_VERIFY(metric_id == 4 || metric_id == 14 || metric_id == 66,
531  "Metric not supported for shifted/pseudo barriers.");
532  }
533  untangler_metric = new TMOP_WorstCaseUntangleOptimizer_Metric(*metric,
534  2,
535  1.5,
536  0.001, //0.01 for pseudo barrier
537  0.001,
538  btype,
539  wctype);
540  }
541
542  if (metric_id < 300 || h_metric_id < 300)
543  {
544  MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
545  }
546  if (metric_id >= 300 || h_metric_id >= 300)
547  {
548  MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
549  }
550
552  TargetConstructor *target_c = NULL;
555  H1_FECollection ind_fec(mesh_poly_deg, dim);
556  FiniteElementSpace ind_fes(mesh, &ind_fec);
557  FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
558  GridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
559  GridFunction aspr3d(&ind_fesv);
560
561  const AssemblyLevel al =
563
564  switch (target_id)
565  {
566  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
567  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
568  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
569  case 4: // Analytic
570  {
573  adapt_coeff = new HessianCoefficient(dim, metric_id);
575  target_c = tc;
576  break;
577  }
578  case 5: // Discrete size 2D or 3D
579  {
583  {
585  }
586  else
587  {
588 #ifdef MFEM_USE_GSLIB
590 #else
591  MFEM_ABORT("MFEM is not built with GSLIB.");
592 #endif
593  }
594  if (dim == 2)
595  {
597  size.ProjectCoefficient(size_coeff);
598  }
599  else if (dim == 3)
600  {
602  size.ProjectCoefficient(size_coeff);
603  }
604  tc->SetSerialDiscreteTargetSize(size);
605  target_c = tc;
606  break;
607  }
608  case 6: // Discrete size + aspect ratio - 2D
609  {
610  GridFunction d_x(&ind_fes), d_y(&ind_fes), disc(&ind_fes);
611
615  disc.ProjectCoefficient(mat_coeff);
617  {
619  }
620  else
621  {
622 #ifdef MFEM_USE_GSLIB
624 #else
625  MFEM_ABORT("MFEM is not built with GSLIB.");
626 #endif
627  }
628
629  // Diffuse the interface
630  DiffuseField(disc,2);
631
632  // Get partials with respect to x and y of the grid function
633  disc.GetDerivative(1,0,d_x);
634  disc.GetDerivative(1,1,d_y);
635
636  // Compute the squared magnitude of the gradient
637  for (int i = 0; i < size.Size(); i++)
638  {
639  size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
640  }
641  const double max = size.Max();
642
643  for (int i = 0; i < d_x.Size(); i++)
644  {
645  d_x(i) = std::abs(d_x(i));
646  d_y(i) = std::abs(d_y(i));
647  }
648  const double eps = 0.01;
649  const double aspr_ratio = 20.0;
650  const double size_ratio = 40.0;
651
652  for (int i = 0; i < size.Size(); i++)
653  {
654  size(i) = (size(i)/max);
655  aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
656  aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
657  if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
658  if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
659  }
660  Vector vals;
661  const int NE = mesh->GetNE();
662  double volume = 0.0, volume_ind = 0.0;
663
664  for (int i = 0; i < NE; i++)
665  {
667  const IntegrationRule &ir =
668  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr->OrderJ());
669  size.GetValues(i, ir, vals);
670  for (int j = 0; j < ir.GetNPoints(); j++)
671  {
672  const IntegrationPoint &ip = ir.IntPoint(j);
673  Tr->SetIntPoint(&ip);
674  volume += ip.weight * Tr->Weight();
675  volume_ind += vals(j) * ip.weight * Tr->Weight();
676  }
677  }
678
679  const double avg_zone_size = volume / NE;
680
681  const double small_avg_ratio = (volume_ind + (volume - volume_ind) /
682  size_ratio) /
683  volume;
684
685  const double small_zone_size = small_avg_ratio * avg_zone_size;
686  const double big_zone_size = size_ratio * small_zone_size;
687
688  for (int i = 0; i < size.Size(); i++)
689  {
690  const double val = size(i);
691  const double a = (big_zone_size - small_zone_size) / small_zone_size;
692  size(i) = big_zone_size / (1.0+a*val);
693  }
694
695  DiffuseField(size, 2);
696  DiffuseField(aspr, 2);
697
698  tc->SetSerialDiscreteTargetSize(size);
700  target_c = tc;
701  break;
702  }
703  case 7: // Discrete aspect ratio 3D
704  {
708  {
710  }
711  else
712  {
713 #ifdef MFEM_USE_GSLIB
715 #else
716  MFEM_ABORT("MFEM is not built with GSLIB.");
717 #endif
718  }
720  aspr3d.ProjectCoefficient(fd_aspr3d);
721
723  target_c = tc;
724  break;
725  }
726  case 8: // shape/size + orientation 2D
727  {
731  {
733  }
734  else
735  {
736 #ifdef MFEM_USE_GSLIB
738 #else
739  MFEM_ABORT("MFEM is not built with GSLIB.");
740 #endif
741  }
742
743  if (metric_id == 14 || metric_id == 36)
744  {
745  ConstantCoefficient size_coeff(0.1*0.1);
746  size.ProjectCoefficient(size_coeff);
747  tc->SetSerialDiscreteTargetSize(size);
748  }
749
750  if (metric_id == 85)
751  {
753  aspr.ProjectCoefficient(aspr_coeff);
754  DiffuseField(aspr,2);
756  }
757
759  ori.ProjectCoefficient(ori_coeff);
761  target_c = tc;
762  break;
763  }
764  // Targets used for hr-adaptivity tests.
765  case 9: // size target in an annular region.
766  case 10: // size+aspect-ratio in an annular region.
767  case 11: // size+aspect-ratio target for a rotate sine wave
768  {
771  hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
773  target_c = tc;
774  break;
775  }
776  default: cout << "Unknown target_id: " << target_id << endl; return 3;
777  }
778  if (target_c == NULL)
779  {
780  target_c = new TargetConstructor(target_t);
781  }
782  target_c->SetNodes(x0);
783  TMOP_QualityMetric *metric_to_use = barrier_type > 0 || worst_case_type > 0
784  ? untangler_metric
785  : metric;
786  TMOP_Integrator *tmop_integ = new TMOP_Integrator(metric_to_use, target_c,
787  h_metric);
788  if (barrier_type > 0 || worst_case_type > 0)
789  {
790  tmop_integ->ComputeUntangleMetricQuantiles(x, *fespace);
791  }
792
793
794  // Finite differences for computations of derivatives.
795  if (fdscheme)
796  {
797  MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
798  tmop_integ->EnableFiniteDifferences(x);
799  }
800  tmop_integ->SetExactActionFlag(exactaction);
801
802  // Setup the quadrature rules for the TMOP integrator.
803  IntegrationRules *irules = NULL;
805  {
806  case 1: irules = &IntRulesLo; break;
807  case 2: irules = &IntRules; break;
808  case 3: irules = &IntRulesCU; break;
809  default: cout << "Unknown quad_type: " << quad_type << endl; return 3;
810  }
812  if (dim == 2)
813  {
814  cout << "Triangle quadrature points: "
817  << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
818  }
819  if (dim == 3)
820  {
821  cout << "Tetrahedron quadrature points: "
823  << "\nHexahedron quadrature points: "
825  << "\nPrism quadrature points: "
826  << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
827  }
828
829  // Limit the node movement.
830  // The limiting distances can be given by a general function of space.
831  FiniteElementSpace dist_fespace(mesh, fec); // scalar space
832  GridFunction dist(&dist_fespace);
833  dist = 1.0;
834  // The small_phys_size is relevant only with proper normalization.
835  if (normalization) { dist = small_phys_size; }
836  ConstantCoefficient lim_coeff(lim_const);
837  if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
838
844  {
845  MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
846
849
851  else if (adapt_eval == 1)
852  {
853 #ifdef MFEM_USE_GSLIB
855 #else
856  MFEM_ABORT("MFEM is not built with GSLIB support!");
857 #endif
858  }
859  else { MFEM_ABORT("Bad interpolation option."); }
860
863  if (visualization)
864  {
865  socketstream vis1;
866  common::VisualizeField(vis1, "localhost", 19916, adapt_lim_gf0, "Zeta 0",
867  300, 600, 300, 300);
868  }
869  }
870
871  // Surface fitting.
872  L2_FECollection mat_coll(0, dim);
873  H1_FECollection surf_fit_fec(mesh_poly_deg, dim);
874  FiniteElementSpace surf_fit_fes(mesh, &surf_fit_fec);
875  FiniteElementSpace mat_fes(mesh, &mat_coll);
876  GridFunction mat(&mat_fes);
877  GridFunction surf_fit_mat_gf(&surf_fit_fes);
878  GridFunction surf_fit_gf0(&surf_fit_fes);
879  Array<bool> surf_fit_marker(surf_fit_gf0.Size());
880  ConstantCoefficient surf_fit_coeff(surface_fit_const);
882  if (surface_fit_const > 0.0)
883  {
885  "Surface fitting with HR is not implemented yet.");
886  MFEM_VERIFY(pa == false,
887  "Surface fitting with PA is not implemented yet.");
888
890  surf_fit_gf0.ProjectCoefficient(ls_coeff);
891
892  for (int i = 0; i < mesh->GetNE(); i++)
893  {
894  mat(i) = material_id(i, surf_fit_gf0);
895  mesh->SetAttribute(i, static_cast<int>(mat(i) + 1));
896  }
897
898  GridFunctionCoefficient mat_coeff(&mat);
899  surf_fit_mat_gf.ProjectDiscCoefficient(mat_coeff, GridFunction::ARITHMETIC);
900  for (int j = 0; j < surf_fit_marker.Size(); j++)
901  {
902  if (surf_fit_mat_gf(j) > 0.1 && surf_fit_mat_gf(j) < 0.9)
903  {
904  surf_fit_marker[j] = true;
905  surf_fit_mat_gf(j) = 1.0;
906  }
907  else
908  {
909  surf_fit_marker[j] = false;
910  surf_fit_mat_gf(j) = 0.0;
911  }
912  }
913
915  else if (adapt_eval == 1)
916  {
917 #ifdef MFEM_USE_GSLIB
919 #else
920  MFEM_ABORT("MFEM is not built with GSLIB support!");
921 #endif
922  }
923  else { MFEM_ABORT("Bad interpolation option."); }
924
925  tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker,
927  if (visualization)
928  {
929  socketstream vis1, vis2, vis3;
930  common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0, "Level Set 0",
931  300, 600, 300, 300);
932  common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
933  600, 600, 300, 300);
934  common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
935  "Dofs to Move",
936  900, 600, 300, 300);
937  }
938  }
939
940  // Has to be after the enabling of the limiting / alignment, as it computes
941  // normalization factors for these terms as well.
942  if (normalization) { tmop_integ->EnableNormalization(x0); }
943
944  // 12. Setup the final NonlinearForm (which defines the integral of interest,
945  // its first and second derivatives). Here we can use a combination of
946  // metrics, i.e., optimize the sum of two integrals, where both are
947  // scaled by used-defined space-dependent weights. Note that there are no
948  // command-line options for the weights and the type of the second
949  // metric; one should update those in the code.
950  NonlinearForm a(fespace);
952  ConstantCoefficient *metric_coeff1 = NULL;
953  TMOP_QualityMetric *metric2 = NULL;
954  TargetConstructor *target_c2 = NULL;
955  FunctionCoefficient metric_coeff2(weight_fun);
956
957  // Explicit combination of metrics.
958  if (combomet > 0)
959  {
960  // First metric.
961  metric_coeff1 = new ConstantCoefficient(1.0);
962  tmop_integ->SetCoefficient(*metric_coeff1);
963
964  // Second metric.
965  if (dim == 2) { metric2 = new TMOP_Metric_077; }
966  else { metric2 = new TMOP_Metric_315; }
967  TMOP_Integrator *tmop_integ2 = NULL;
968  if (combomet == 1)
969  {
970  target_c2 = new TargetConstructor(
972  target_c2->SetVolumeScale(0.01);
973  target_c2->SetNodes(x0);
974  tmop_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
975  tmop_integ2->SetCoefficient(metric_coeff2);
976  }
977  else { tmop_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
979  if (fdscheme) { tmop_integ2->EnableFiniteDifferences(x); }
980  tmop_integ2->SetExactActionFlag(exactaction);
981
985  if (normalization) { combo->EnableNormalization(x0); }
986  if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
987
989  }
990  else
991  {
993  }
994
995  if (pa) { a.Setup(); }
996
997  // Compute the minimum det(J) of the starting mesh.
998  min_detJ = infinity();
999  const int NE = mesh->GetNE();
1000  for (int i = 0; i < NE; i++)
1001  {
1002  const IntegrationRule &ir =
1005  for (int j = 0; j < ir.GetNPoints(); j++)
1006  {
1007  transf->SetIntPoint(&ir.IntPoint(j));
1008  min_detJ = min(min_detJ, transf->Jacobian().Det());
1009  }
1010  }
1011  cout << "Minimum det(J) of the original mesh is " << min_detJ << endl;
1012
1013  if (min_detJ < 0.0 && barrier_type == 0
1014  && metric_id != 22 && metric_id != 211 && metric_id != 252
1015  && metric_id != 311 && metric_id != 313 && metric_id != 352)
1016  {
1017  MFEM_ABORT("The input mesh is inverted! Try an untangling metric.");
1018  }
1019  if (min_detJ < 0.0)
1020  {
1021  MFEM_VERIFY(target_t == TargetConstructor::IDEAL_SHAPE_UNIT_SIZE,
1022  "Untangling is supported only for ideal targets.");
1023
1024  const DenseMatrix &Wideal =
1026  min_detJ /= Wideal.Det();
1027
1028  // Slightly below minJ0 to avoid div by 0.
1029  min_detJ -= 0.01 * h0.Min();
1030  }
1031
1032  // For HR tests, the energy is normalized by the number of elements.
1033  const double init_energy = a.GetGridFunctionEnergy(x) /
1034  (hradaptivity ? mesh->GetNE() : 1);
1035  double init_metric_energy = init_energy;
1036  if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
1037  {
1038  lim_coeff.constant = 0.0;
1040  surf_fit_coeff.constant = 0.0;
1041  init_metric_energy = a.GetGridFunctionEnergy(x) /
1042  (hradaptivity ? mesh->GetNE() : 1);
1043  lim_coeff.constant = lim_const;
1045  surf_fit_coeff.constant = surface_fit_const;
1046  }
1047
1048  // Visualize the starting mesh and metric values.
1049  // Note that for combinations of metrics, this only shows the first metric.
1050  if (visualization)
1051  {
1052  char title[] = "Initial metric values";
1053  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
1054  }
1055
1056  // 13. Fix all boundary nodes, or fix only a given component depending on the
1057  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
1058  // fixed x/y/z components of the node. Attribute 4 corresponds to an
1059  // entirely fixed node. Other boundary attributes do not affect the node
1060  // movement boundary conditions.
1061  if (move_bnd == false)
1062  {
1063  Array<int> ess_bdr(mesh->bdr_attributes.Max());
1064  ess_bdr = 1;
1065  a.SetEssentialBC(ess_bdr);
1066  }
1067  else
1068  {
1069  int n = 0;
1070  for (int i = 0; i < mesh->GetNBE(); i++)
1071  {
1072  const int nd = fespace->GetBE(i)->GetDof();
1073  const int attr = mesh->GetBdrElement(i)->GetAttribute();
1074  MFEM_VERIFY(!(dim == 2 && attr == 3),
1075  "Boundary attribute 3 must be used only for 3D meshes. "
1076  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
1078  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
1079  if (attr == 4) { n += nd * dim; }
1080  }
1081  Array<int> ess_vdofs(n);
1082  n = 0;
1083  for (int i = 0; i < mesh->GetNBE(); i++)
1084  {
1085  const int nd = fespace->GetBE(i)->GetDof();
1086  const int attr = mesh->GetBdrElement(i)->GetAttribute();
1087  fespace->GetBdrElementVDofs(i, vdofs);
1088  if (attr == 1) // Fix x components.
1089  {
1090  for (int j = 0; j < nd; j++)
1091  { ess_vdofs[n++] = vdofs[j]; }
1092  }
1093  else if (attr == 2) // Fix y components.
1094  {
1095  for (int j = 0; j < nd; j++)
1096  { ess_vdofs[n++] = vdofs[j+nd]; }
1097  }
1098  else if (attr == 3) // Fix z components.
1099  {
1100  for (int j = 0; j < nd; j++)
1101  { ess_vdofs[n++] = vdofs[j+2*nd]; }
1102  }
1103  else if (attr == 4) // Fix all components.
1104  {
1105  for (int j = 0; j < vdofs.Size(); j++)
1106  { ess_vdofs[n++] = vdofs[j]; }
1107  }
1108  }
1109  a.SetEssentialVDofs(ess_vdofs);
1110  }
1111
1112  // 14. As we use the Newton method to solve the resulting nonlinear system,
1113  // here we setup the linear solver for the system's Jacobian.
1114  Solver *S = NULL, *S_prec = NULL;
1115  const double linsol_rtol = 1e-12;
1116  if (lin_solver == 0)
1117  {
1118  S = new DSmoother(1, 1.0, max_lin_iter);
1119  }
1120  else if (lin_solver == 1)
1121  {
1122  CGSolver *cg = new CGSolver;
1123  cg->SetMaxIter(max_lin_iter);
1124  cg->SetRelTol(linsol_rtol);
1125  cg->SetAbsTol(0.0);
1126  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
1127  S = cg;
1128  }
1129  else
1130  {
1131  MINRESSolver *minres = new MINRESSolver;
1132  minres->SetMaxIter(max_lin_iter);
1133  minres->SetRelTol(linsol_rtol);
1134  minres->SetAbsTol(0.0);
1135  if (verbosity_level > 2) { minres->SetPrintLevel(1); }
1136  minres->SetPrintLevel(verbosity_level == 2 ? 3 : -1);
1137  if (lin_solver == 3 || lin_solver == 4)
1138  {
1139  if (pa)
1140  {
1141  MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
1142  auto js = new OperatorJacobiSmoother;
1143  js->SetPositiveDiagonal(true);
1144  S_prec = js;
1145  }
1146  else
1147  {
1148  auto ds = new DSmoother((lin_solver == 3) ? 0 : 1, 1.0, 1);
1149  ds->SetPositiveDiagonal(true);
1150  S_prec = ds;
1151  }
1152  minres->SetPreconditioner(*S_prec);
1153  }
1154  S = minres;
1155  }
1156
1157  // Perform the nonlinear optimization.
1158  const IntegrationRule &ir =
1160  TMOPNewtonSolver solver(ir, solver_type);
1162  if (surface_fit_threshold > 0)
1163  {
1164  solver.SetTerminationWithMaxSurfaceFittingError(surface_fit_threshold);
1165  }
1166  // Provide all integration rules in case of a mixed mesh.
1168  if (solver_type == 0)
1169  {
1170  // Specify linear solver when we use a Newton-based solver.
1171  solver.SetPreconditioner(*S);
1172  }
1173  // For untangling, the solver will update the min det(T) values.
1174  solver.SetMinDetPtr(&min_detJ);
1175  solver.SetMaxIter(solver_iter);
1176  solver.SetRelTol(solver_rtol);
1177  solver.SetAbsTol(0.0);
1178  if (solver_art_type > 0)
1179  {
1181  }
1182  solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
1183
1185  // If hr-adaptivity is disabled, r-adaptivity is done once using the
1186  // TMOPNewtonSolver.
1187  // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
1189  // The solver terminates if an h-adaptivity iteration does not modify
1190  // any element in the mesh.
1191  TMOPHRSolver hr_solver(*mesh, a, solver,
1193  mesh_poly_deg, h_metric_id,
1194  n_hr_iter, n_h_iter);
1197  {
1200  }
1201  hr_solver.Mult();
1202
1203  // 15. Save the optimized mesh to a file. This output can be viewed later
1204  // using GLVis: "glvis -m optimized.mesh".
1205  {
1206  ofstream mesh_ofs("optimized.mesh");
1207  mesh_ofs.precision(14);
1208  mesh->Print(mesh_ofs);
1209  }
1210
1211  const double fin_energy = a.GetGridFunctionEnergy(x) /
1212  (hradaptivity ? mesh->GetNE() : 1);
1213  double fin_metric_energy = fin_energy;
1214  if (lim_const > 0.0 || adapt_lim_const > 0.0)
1215  {
1216  lim_coeff.constant = 0.0;
1218  surf_fit_coeff.constant = 0.0;
1219  fin_metric_energy = a.GetGridFunctionEnergy(x) /
1220  (hradaptivity ? mesh->GetNE() : 1);
1221  lim_coeff.constant = lim_const;
1223  surf_fit_coeff.constant = surface_fit_const;
1224  }
1225  std::cout << std::scientific << std::setprecision(4);
1226  cout << "Initial strain energy: " << init_energy
1227  << " = metrics: " << init_metric_energy
1228  << " + extra terms: " << init_energy - init_metric_energy << endl;
1229  cout << " Final strain energy: " << fin_energy
1230  << " = metrics: " << fin_metric_energy
1231  << " + extra terms: " << fin_energy - fin_metric_energy << endl;
1232  cout << "The strain energy decreased by: "
1233  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
1234
1235  // 16. Visualize the final mesh and metric values.
1236  if (visualization)
1237  {
1238  char title[] = "Final metric values";
1239  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
1240  }
1241
1242  if (adapt_lim_const > 0.0 && visualization)
1243  {
1244  socketstream vis0;
1245  common::VisualizeField(vis0, "localhost", 19916, adapt_lim_gf0, "Xi 0",
1246  600, 600, 300, 300);
1247  }
1248
1249  if (surface_fit_const > 0.0)
1250  {
1251  if (visualization)
1252  {
1253  socketstream vis2, vis3;
1254  common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
1255  600, 900, 300, 300);
1256  common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf, "Surface dof",
1257  900, 900, 300, 300);
1258  }
1259  double err_avg, err_max;
1260  tmop_integ->GetSurfaceFittingErrors(err_avg, err_max);
1261  std::cout << "Avg fitting error: " << err_avg << std::endl
1262  << "Max fitting error: " << err_max << std::endl;
1263  }
1264
1265  // 17. Visualize the mesh displacement.
1266  if (visualization)
1267  {
1268  osockstream sock(19916, "localhost");
1269  sock << "solution\n";
1270  mesh->Print(sock);
1271  x0 -= x;
1272  x0.Save(sock);
1273  sock.send();
1274  sock << "window_title 'Displacements'\n"
1275  << "window_geometry "
1276  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
1277  << "keys jRmclA" << endl;
1278  }
1279
1280  delete S;
1281  delete S_prec;
1282  delete target_c2;
1283  delete metric2;
1284  delete metric_coeff1;
1287  delete target_c;
1290  delete h_metric;
1291  delete metric;
1292  delete untangler_metric;
1293  delete fespace;
1294  delete fec;
1295  delete mesh;
1296
1297  return 0;
1298 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
void ComputeUntangleMetricQuantiles(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:4125
void discrete_aspr_3d(const Vector &x, Vector &v)
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Definition: solvers.hpp:465
Definition: ex25.cpp:148
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
3D non-barrier Shape (S) metric.
Definition: tmop.hpp:939
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:3993
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:251
double discrete_size_2d(const Vector &x)
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:744
2D non-barrier Shape+Size (VS) metric.
Definition: tmop.hpp:452
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:923
void SetPositiveDiagonal(bool pos_diag=true)
Replace diagonal entries with their absolute values.
Definition: solvers.hpp:311
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:197
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:152
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:872
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
Definition: tmop.cpp:1802
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
Definition: tmop.cpp:3829
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:926
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:625
Definition: tmop.hpp:1495
virtual void Setup()
Setup the NonlinearForm: based on the current AssemblyLevel and the current mesh, optionally...
double Det() const
Definition: densemat.cpp:449
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2573
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:493
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
double surface_level_set(const Vector &x)
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Container class for integration rules.
Definition: intrules.hpp:311
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:604
void GetSurfaceFittingErrors(double &err_avg, double &err_max)
Definition: tmop.cpp:2786
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:1252
2D barrier shape (S) metric (polyconvex).
Definition: tmop.hpp:275
Restriction of the node positions to certain regions.
Definition: tmop.cpp:2715
3D non-barrier metric without a type.
Definition: tmop.hpp:707
2D barrier size (V) metric (polyconvex).
Definition: tmop.hpp:475
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop.hpp:1823
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
Definition: tmop.cpp:4264
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:1033
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:786
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
Definition: tmop_tools.cpp:927
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
Definition: tmop.hpp:1978
MINRES method.
Definition: solvers.hpp:575
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:785
virtual int OrderJ() const =0
Return the order of the elements of the Jacobian of the transformation.
void SetMinDetPtr(double *md_ptr)
Definition: tmop_tools.hpp:203
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:314
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1069
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:765
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:4165
2D Shifted barrier form of shape metric (mu_2).
Definition: tmop.hpp:363
virtual void SetSerialDiscreteTargetAspectRatio(const GridFunction &tspec_)
Definition: tmop.cpp:1792
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:297
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:433
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:646
void SetTerminationWithMaxSurfaceFittingError(double max_error)
Definition: tmop_tools.hpp:224
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3676
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
2D non-barrier shape (S) metric.
Definition: tmop.hpp:296
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
Geometry Geometries
Definition: fe.cpp:49
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void GetDerivative(int comp, int der_comp, GridFunction &der)
Compute a certain derivative of a function&#39;s component. Derivatives of the function are computed at t...
Definition: gridfunc.cpp:1424
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:146
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:588
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:76
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Limiting of the mesh displacements (general version).
Definition: tmop.cpp:2690
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:623
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Definition: tmop.hpp:1835
2D non-barrier size (V) metric (not polyconvex).
Definition: tmop.hpp:398
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9473
3D barrier metric without a type.
Definition: tmop.hpp:725
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:1246
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:896
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:720
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:131
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:531
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5285
Definition: tmop_amr.hpp:253
int Dimension() const
Definition: mesh.hpp:1006
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
Definition: tmop_tools.hpp:220
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:1015
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2678
A general vector function coefficient.
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:807
void SetAbsTol(double atol)
Definition: solvers.hpp:199
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
Definition: tmop_amr.hpp:252
void SetRelTol(double rtol)
Definition: solvers.hpp:198
2D barrier (not a shape) metric (polyconvex).
Definition: tmop.hpp:382
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: tmop_tools.hpp:254
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
double weight_fun(const Vector &x)
2D barrier size (V) metric (polyconvex).
Definition: tmop.hpp:415
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
void SetAdaptiveLinRtol(const int type=2, const double rtol0=0.5, const double rtol_max=0.9, const double alpha=0.5 *(1.0+sqrt(5.0)), const double gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
Definition: solvers.cpp:1922
Adds a new TMOP_Integrator to the combination.
Definition: tmop.hpp:2008
double material_indicator_2d(const Vector &x)
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:332
double a
Definition: lissajous.cpp:41
double discrete_size_3d(const Vector &x)
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 energy corresponding to the state x.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1671
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:997
3D Shape (S) metric, untangling version of 303.
Definition: tmop.hpp:686
int dim
Definition: ex24.cpp:53
virtual 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:2781
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2399
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
A general function coefficient.
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1483
Vector data type.
Definition: vector.hpp:60
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:92
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:851
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:1170
void EnableSurfaceFitting(const GridFunction &s0, const Array< bool > &smarker, Coefficient &coeff, AdaptivityEvaluator &ae)
Fitting of certain DOFs to the zero level set of a function.
Definition: tmop.cpp:2750
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Specify essential boundary conditions.
Base class for solvers.
Definition: operator.hpp:655
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:3100
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:5332
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1166
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, TMOPMatrixCoefficient *mspec)
Definition: tmop.cpp:1543
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:828
int material_id(int el_id, const GridFunction &g)
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
Definition: tmop.cpp:1772
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:516
Definition: vector.hpp:469
int main()
double GetElementVolume(int i)
Definition: mesh.cpp:112
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1044
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:585
2D non-barrier metric without a type.
Definition: tmop.hpp:197
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1565
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:348