MFEM  v4.5.2 Finite element discretization library
pmesh-optimizer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 - Parallel Version
14 // ---------------------------------------------------------------------
15 //
16 // This miniapp performs mesh optimization using the Target-Matrix Optimization
17 // Paradigm (TMOP) by P.Knupp et al., and a global variational minimization
18 // approach. It minimizes the quantity sum_T int_T mu(J(x)), where T are the
19 // target (ideal) elements, J is the Jacobian of the transformation from the
20 // target to the physical element, and mu is the mesh quality metric. This
21 // metric can measure shape, size or alignment of the region around each
22 // quadrature point. The combination of targets & quality metrics is used to
23 // optimize the physical node positions, i.e., they must be as close as possible
24 // to the shape / size / alignment of their targets. This code also demonstrates
25 // a possible use of nonlinear operators (the class TMOP_QualityMetric, defining
26 // mu(J), and the class TMOP_Integrator, defining int mu(J)), as well 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 pmesh-optimizer
33 //
34 // Sample runs:
36 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -bnd -qt 1 -qo 8
38 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 200 -bnd -qt 1 -qo 8 -fd
40 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
41 //
43 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 9 -ni 50 -li 20 -hmid 55 -mid 7 -hr
44 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
45 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
46 //
48 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
49 // (requires GSLIB):
50 // * mpirun -np 4 pmesh-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 // mpirun -np 4 pmesh-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 // * mpirun -n 4 pmesh-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 // mpirun -np 4 pmesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
58 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
59 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
61 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -nor
62 // Adapted discrete aspect ratio (3D):
63 // mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -bnd -qt 1 -qo 8
64 //
66 // mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
67 // Adaptive limiting through the L-BFGS solver:
68 // mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
69 // Adaptive limiting through FD (requires GSLIB):
70 // * mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
71 //
73 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 3 -rs 1 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 5e4 -rtol 1e-5
74 // mpirun -np 4 pmesh-optimizer -m square01-tri.mesh -o 3 -rs 0 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 1e4 -rtol 1e-5
75 // Surface fitting with weight adaptation and termination based on fitting error
76 // mpirun -np 4 pmesh-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
77 //
79 // mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
80 // (requires CUDA):
81 // * mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8 -d cuda
82 // Blade shape with FD-based solver:
83 // mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
85 // mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
86 // ICF shape and equal size:
87 // mpirun -np 4 pmesh-optimizer -o 3 -mid 80 -bec -tid 2 -ni 25 -ls 3 -art 2 -qo 5
88 // ICF shape and initial size:
89 // mpirun -np 4 pmesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
90 // ICF shape:
91 // mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
92 // ICF limited shape:
93 // mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8 -lc 10
94 // ICF combo shape + size (rings, slow convergence):
95 // mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 1000 -bnd -qt 1 -qo 8 -cmb 1
96 // Mixed tet / cube / hex mesh with limiting:
97 // mpirun -np 4 pmesh-optimizer -m ../../data/fichera-mixed-p2.mesh -o 4 -rs 1 -mid 301 -tid 1 -fix-bnd -qo 6 -nor -lc 0.25
98 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
99 // * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -mid 303 -tid 1 -ni 20 -li 500 -fix-bnd
100 // 2D non-conforming shape and equal size:
101 // mpirun -np 4 pmesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -bnd -qt 1 -qo 8
102 //
103 // 2D untangling:
104 // mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
105 // 2D untangling with shifted barrier metric:
106 // mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
107 // 3D untangling (the mesh is in the mfem/data GitHub repository):
108 // * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
109 // Shape optimization for a Kershaw transformed mesh using partial assembly:
110 // Mesh for Kershaw transformation must be a Cartesian mesh with nx % 6 = ny % 2 = nz % 2 = 0.
111 // Kershaw transformation can be imposed using the transformation ('t') feature in the mesh-explorer miniapp.
112 // * mpirun - np 6 pmesh-optimizer -m kershaw-24x24x24.mesh -mid 303 -tid 1 -bnd -ni 100 -art 1 -ls 3 -qo 8 -li 40 -o 2 -qo 8 -ker -pa
113
114 #include "mfem.hpp"
115 #include "../common/mfem-common.hpp"
116 #include <iostream>
117 #include <fstream>
118 #include "mesh-optimizer.hpp"
119
120 using namespace mfem;
121 using namespace std;
122
123 int main (int argc, char *argv[])
124 {
125  // 0. Initialize MPI and HYPRE.
126  Mpi::Init(argc, argv);
127  int myid = Mpi::WorldRank();
128  Hypre::Init();
129
130  // 1. Set the method's default parameters.
131  const char *mesh_file = "icf.mesh";
132  int mesh_poly_deg = 1;
133  int rs_levels = 0;
134  int rp_levels = 0;
135  double jitter = 0.0;
136  int metric_id = 1;
137  int target_id = 1;
138  double lim_const = 0.0;
140  double surface_fit_const = 0.0;
143  int solver_type = 0;
144  int solver_iter = 20;
145  double solver_rtol = 1e-10;
146  int solver_art_type = 0;
147  int lin_solver = 2;
148  int max_lin_iter = 100;
149  bool move_bnd = true;
150  int combomet = 0;
151  bool bal_expl_combo = false;
153  int h_metric_id = -1;
154  bool normalization = false;
155  bool visualization = true;
156  int verbosity_level = 0;
157  bool fdscheme = false;
159  bool exactaction = false;
160  const char *devopt = "cpu";
161  bool pa = false;
162  int n_hr_iter = 5;
163  int n_h_iter = 1;
165  double surface_fit_threshold = -10;
166  int mesh_node_ordering = 0;
167  int barrier_type = 0;
168  int worst_case_type = 0;
169
170  // 2. Parse command-line options.
171  OptionsParser args(argc, argv);
173  "Mesh file to use.");
175  "Polynomial degree of mesh finite element space.");
177  "Number of times to refine the mesh uniformly in serial.");
179  "Number of times to refine the mesh uniformly in parallel.");
181  "Random perturbation scaling factor.");
183  "Mesh optimization metric:\n\t"
184  "T-metrics\n\t"
185  "1 : |T|^2 -- 2D no type\n\t"
186  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
187  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
188  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
189  "14 : |T-I|^2 -- 2D shape+size+orientation\n\t"
190  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
191  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
192  "55 : (tau-1)^2 -- 2D size\n\t"
193  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
194  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
195  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
196  "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t"
197  "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t"
198  "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t"
199  // "211: (tau-1)^2-tau+sqrt(tau^2+eps) -- 2D untangling\n\t"
200  // "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
201  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
202  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
203  "303: (|T|^2)/3/tau^(2/3)-1 -- 3D shape\n\t"
204  "304: (|T|^3)/3^{3/2}/tau-1 -- 3D shape\n\t"
205  // "311: (tau-1)^2-tau+sqrt(tau^2+eps)-- 3D untangling\n\t"
206  "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t"
207  "315: (tau-1)^2 -- 3D no type\n\t"
208  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D no type\n\t"
209  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
210  "322: |T-adjT^-t|^2 -- 3D shape+size\n\t"
211  "323: |J|^3-3sqrt(3)ln(det(J))-3sqrt(3) -- 3D shape+size\n\t"
212  "328: (1-gamma) mu_301 + gamma mu_316 -- 3D shape+size\n\t"
213  "332: (1-gamma) mu_302 + gamma mu_315 -- 3D shape+size\n\t"
214  "333: (1-gamma) mu_302 + gamma mu_316 -- 3D shape+size\n\t"
215  "334: (1-gamma) mu_303 + gamma mu_316 -- 3D shape+size\n\t"
216  "347: (1-gamma) mu_304 + gamma mu_316 -- 3D shape+size\n\t"
217  // "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling\n\t"
218  "360: (|T|^3)/3^{3/2}-tau -- 3D shape\n\t"
219  "A-metrics\n\t"
220  "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t"
221  "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t"
222  "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t"
223  "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t"
224  );
226  "Target (ideal element) type:\n\t"
227  "1: Ideal shape, unit size\n\t"
228  "2: Ideal shape, equal size\n\t"
229  "3: Ideal shape, initial size\n\t"
230  "4: Given full analytic Jacobian (in physical space)\n\t"
231  "5: Ideal shape, given size (in physical space)");
232  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
236  "Surface preservation constant.");
239  "1: Gauss-Lobatto\n\t"
240  "2: Gauss-Legendre\n\t"
241  "3: Closed uniform points");
243  "Order of the quadrature rule.");
245  " Type of solver: (default) 0: Newton, 1: LBFGS");
247  "Maximum number of Newton iterations.");
249  "Relative tolerance for the Newton solver.");
251  "Type of adaptive relative linear solver tolerance:\n\t"
252  "0: None (default)\n\t"
253  "1: Eisenstat-Walker type 1\n\t"
254  "2: Eisenstat-Walker type 2");
256  "Linear solver:\n\t"
257  "0: l1-Jacobi\n\t"
258  "1: CG\n\t"
259  "2: MINRES\n\t"
260  "3: MINRES + Jacobi preconditioner\n\t"
261  "4: MINRES + l1-Jacobi preconditioner");
263  "Maximum number of iterations in the linear solve.");
265  "--fix-boundary",
266  "Enable motion along horizontal and vertical boundaries.");
268  "Combination of metrics options:\n\t"
269  "0: Use single metric\n\t"
270  "1: Shape + space-dependent size given analytically\n\t"
271  "2: Shape + adapted size given discretely; shared target");
273  "-no-bec", "--balance-explicit-combo",
274  "Automatic balancing of explicit combo metrics.");
279  "Same options as metric_id. Used to determine refinement"
280  " type for each element if h-adaptivity is enabled.");
282  "--no-normalization",
283  "Make all terms in the optimization functional unitless.");
285  "-no-fd", "--no-fd-approx",
286  "Enable finite difference based derivative computations.");
288  "-no-ex", "--no-exact-action",
289  "Enable exact action of TMOP_Integrator.");
291  "--no-visualization",
292  "Enable or disable GLVis visualization.");
294  "Set the verbosity level - 0, 1, or 2.");
296  "0 - Advection based (DEFAULT), 1 - GSLIB.");
298  "Device configuration string, see Device::Configure().");
300  "--no-partial-assembly", "Enable Partial Assembly.");
305  "iteration.");
308  "Enable or disable adaptive surface fitting.");
310  "Set threshold for surface fitting. TMOP solver will"
311  "terminate when max surface fitting error is below this limit");
313  "Ordering of mesh nodes."
314  "0 (default): byNodes, 1: byVDIM");
316  "0 - None,"
317  "1 - Shifted Barrier,"
318  "2 - Pseudo Barrier.");
320  "0 - None,"
321  "1 - Beta,"
322  "2 - PMean.");
323
324  args.Parse();
325  if (!args.Good())
326  {
327  if (myid == 0) { args.PrintUsage(cout); }
328  return 1;
329  }
330  if (myid == 0) { args.PrintOptions(cout); }
331  if (h_metric_id < 0) { h_metric_id = metric_id; }
332
334  {
335  MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
336  " supported on cpus.");
337  }
338  Device device(devopt);
339  if (myid == 0) { device.Print();}
340
341  // 3. Initialize and refine the starting mesh.
342  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
343  for (int lev = 0; lev < rs_levels; lev++)
344  {
345  mesh->UniformRefinement();
346  }
347  const int dim = mesh->Dimension();
348
349  if (hradaptivity) { mesh->EnsureNCMesh(); }
350  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
351
352  delete mesh;
353  for (int lev = 0; lev < rp_levels; lev++)
354  {
355  pmesh->UniformRefinement();
356  }
357
358  // 4. Define a finite element space on the mesh. Here we use vector finite
359  // elements which are tensor products of quadratic finite elements. The
360  // number of components in the vector finite element space is specified by
361  // the last parameter of the FiniteElementSpace constructor.
363  if (mesh_poly_deg <= 0)
364  {
366  mesh_poly_deg = 2;
367  }
368  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
369  ParFiniteElementSpace *pfespace = new ParFiniteElementSpace(pmesh, fec, dim,
370  mesh_node_ordering);
371
372  // 5. Make the mesh curved based on the above finite element space. This
373  // means that we define the mesh elements through a fespace-based
374  // transformation of the reference element.
375  pmesh->SetNodalFESpace(pfespace);
376
377  // 7. Get the mesh nodes (vertices and other degrees of freedom in the finite
378  // element space) as a finite element grid function in fespace. Note that
379  // changing x automatically changes the shapes of the mesh elements.
380  ParGridFunction x(pfespace);
381  pmesh->SetNodalGridFunction(&x);
382
383  // 8. Define a vector representing the minimal local mesh size in the mesh
384  // nodes. We index the nodes using the scalar version of the degrees of
385  // freedom in pfespace. Note: this is partition-dependent.
386  //
387  // In addition, compute average mesh size and total volume.
388  Vector h0(pfespace->GetNDofs());
389  h0 = infinity();
390  double vol_loc = 0.0;
391  Array<int> dofs;
392  for (int i = 0; i < pmesh->GetNE(); i++)
393  {
394  // Get the local scalar element degrees of freedom in dofs.
395  pfespace->GetElementDofs(i, dofs);
396  // Adjust the value of h0 in dofs based on the local mesh size.
397  const double hi = pmesh->GetElementSize(i);
398  for (int j = 0; j < dofs.Size(); j++)
399  {
400  h0(dofs[j]) = min(h0(dofs[j]), hi);
401  }
402  vol_loc += pmesh->GetElementVolume(i);
403  }
404  double vol_glb;
405  MPI_Allreduce(&vol_loc, &vol_glb, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
406  const double small_phys_size = pow(vol_glb, 1.0 / dim) / 100.0;
407
408  // 9. Add a random perturbation to the nodes in the interior of the domain.
409  // We define a random grid function of fespace and make sure that it is
410  // zero on the boundary and its values are locally of the order of h0.
411  // The latter is based on the DofToVDof() method which maps the scalar to
412  // the vector degrees of freedom in pfespace.
413  ParGridFunction rdm(pfespace);
414  rdm.Randomize();
415  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
416  rdm *= jitter;
418  // Scale the random values to be of order of the local mesh size.
419  for (int i = 0; i < pfespace->GetNDofs(); i++)
420  {
421  for (int d = 0; d < dim; d++)
422  {
423  rdm(pfespace->DofToVDof(i,d)) *= h0(i);
424  }
425  }
426  Array<int> vdofs;
427  for (int i = 0; i < pfespace->GetNBE(); i++)
428  {
429  // Get the vector degrees of freedom in the boundary element.
430  pfespace->GetBdrElementVDofs(i, vdofs);
431  // Set the boundary values to zero.
432  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
433  }
434  x -= rdm;
435  // Set the perturbation of all nodes from the true nodes.
436  x.SetTrueVector();
437  x.SetFromTrueVector();
438
439  // 10. Save the starting (prior to the optimization) mesh to a file. This
440  // output can be viewed later using GLVis: "glvis -m perturbed -np
442  {
443  ostringstream mesh_name;
444  mesh_name << "perturbed.mesh";
445  ofstream mesh_ofs(mesh_name.str().c_str());
446  mesh_ofs.precision(8);
447  pmesh->PrintAsOne(mesh_ofs);
448  }
449
450  // 11. Store the starting (prior to the optimization) positions.
451  ParGridFunction x0(pfespace);
452  x0 = x;
453
454  // 12. Form the integrator that uses the chosen metric and target.
455  double min_detJ = -0.1;
456  TMOP_QualityMetric *metric = NULL;
457  switch (metric_id)
458  {
459  // T-metrics
460  case 1: metric = new TMOP_Metric_001; break;
461  case 2: metric = new TMOP_Metric_002; break;
462  case 4: metric = new TMOP_Metric_004; break;
463  case 7: metric = new TMOP_Metric_007; break;
464  case 9: metric = new TMOP_Metric_009; break;
465  case 14: metric = new TMOP_Metric_014; break;
466  case 22: metric = new TMOP_Metric_022(min_detJ); break;
467  case 50: metric = new TMOP_Metric_050; break;
468  case 55: metric = new TMOP_Metric_055; break;
469  case 56: metric = new TMOP_Metric_056; break;
470  case 58: metric = new TMOP_Metric_058; break;
471  case 66: metric = new TMOP_Metric_066(0.5); break;
472  case 77: metric = new TMOP_Metric_077; break;
473  case 80: metric = new TMOP_Metric_080(0.5); break;
474  case 85: metric = new TMOP_Metric_085; break;
475  case 98: metric = new TMOP_Metric_098; break;
476  // case 211: metric = new TMOP_Metric_211; break;
477  // case 252: metric = new TMOP_Metric_252(min_detJ); break;
478  case 301: metric = new TMOP_Metric_301; break;
479  case 302: metric = new TMOP_Metric_302; break;
480  case 303: metric = new TMOP_Metric_303; break;
481  case 304: metric = new TMOP_Metric_304; break;
482  // case 311: metric = new TMOP_Metric_311; break;
483  case 313: metric = new TMOP_Metric_313(min_detJ); break;
484  case 315: metric = new TMOP_Metric_315; break;
485  case 316: metric = new TMOP_Metric_316; break;
486  case 321: metric = new TMOP_Metric_321; break;
487  case 322: metric = new TMOP_Metric_322; break;
488  case 323: metric = new TMOP_Metric_323; break;
489  case 328: metric = new TMOP_Metric_328(0.5); break;
490  case 332: metric = new TMOP_Metric_332(0.5); break;
491  case 333: metric = new TMOP_Metric_333(0.5); break;
492  case 334: metric = new TMOP_Metric_334(0.5); break;
493  case 347: metric = new TMOP_Metric_347(0.5); break;
494  // case 352: metric = new TMOP_Metric_352(min_detJ); break;
495  case 360: metric = new TMOP_Metric_360; break;
496  // A-metrics
497  case 11: metric = new TMOP_AMetric_011; break;
498  case 36: metric = new TMOP_AMetric_036; break;
499  case 107: metric = new TMOP_AMetric_107a; break;
500  case 126: metric = new TMOP_AMetric_126(0.9); break;
501  default:
502  if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
503  return 3;
504  }
505  TMOP_QualityMetric *h_metric = NULL;
507  {
508  switch (h_metric_id)
509  {
510  case 1: h_metric = new TMOP_Metric_001; break;
511  case 2: h_metric = new TMOP_Metric_002; break;
512  case 7: h_metric = new TMOP_Metric_007; break;
513  case 9: h_metric = new TMOP_Metric_009; break;
514  case 55: h_metric = new TMOP_Metric_055; break;
515  case 56: h_metric = new TMOP_Metric_056; break;
516  case 58: h_metric = new TMOP_Metric_058; break;
517  case 77: h_metric = new TMOP_Metric_077; break;
518  case 315: h_metric = new TMOP_Metric_315; break;
519  case 316: h_metric = new TMOP_Metric_316; break;
520  case 321: h_metric = new TMOP_Metric_321; break;
521  default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
522  endl;
523  return 3;
524  }
525  }
526
528  switch (barrier_type)
529  {
531  break;
533  break;
535  break;
536  default: cout << "barrier_type not supported: " << barrier_type << endl;
537  return 3;
538  }
539
541  switch (worst_case_type)
542  {
544  break;
546  break;
548  break;
549  default: cout << "worst_case_type not supported: " << worst_case_type << endl;
550  return 3;
551  }
552
553  TMOP_QualityMetric *untangler_metric = NULL;
554  if (barrier_type > 0 || worst_case_type > 0)
555  {
556  if (barrier_type > 0)
557  {
558  MFEM_VERIFY(metric_id == 4 || metric_id == 14 || metric_id == 66,
559  "Metric not supported for shifted/pseudo barriers.");
560  }
561  untangler_metric = new TMOP_WorstCaseUntangleOptimizer_Metric(*metric,
562  2,
563  1.5,
564  0.001,//0.01 for pseudo barrier
565  0.001,
566  btype,
567  wctype);
568  }
569
570  if (metric_id < 300 || h_metric_id < 300)
571  {
572  MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
573  }
574  if (metric_id >= 300 || h_metric_id >= 300)
575  {
576  MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
577  }
578
580  TargetConstructor *target_c = NULL;
583  H1_FECollection ind_fec(mesh_poly_deg, dim);
584  ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
585  ParFiniteElementSpace ind_fesv(pmesh, &ind_fec, dim);
586  ParGridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
587  ParGridFunction aspr3d(&ind_fesv);
588
589  const AssemblyLevel al =
591
592  switch (target_id)
593  {
594  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
595  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
596  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
597  case 4:
598  {
601  adapt_coeff = new HessianCoefficient(dim, metric_id);
603  target_c = tc;
604  break;
605  }
606  case 5: // Discrete size 2D or 3D
607  {
611  {
613  }
614  else
615  {
616 #ifdef MFEM_USE_GSLIB
618 #else
619  MFEM_ABORT("MFEM is not built with GSLIB.");
620 #endif
621  }
622  if (dim == 2)
623  {
625  size.ProjectCoefficient(size_coeff);
626  }
627  else if (dim == 3)
628  {
630  size.ProjectCoefficient(size_coeff);
631  }
632  tc->SetParDiscreteTargetSize(size);
633  target_c = tc;
634  break;
635  }
636  case 6: // material indicator 2D
637  {
638  ParGridFunction d_x(&ind_fes), d_y(&ind_fes), disc(&ind_fes);
639
643  disc.ProjectCoefficient(mat_coeff);
645  {
647  }
648  else
649  {
650 #ifdef MFEM_USE_GSLIB
652 #else
653  MFEM_ABORT("MFEM is not built with GSLIB.");
654 #endif
655  }
656  // Diffuse the interface
657  DiffuseField(disc,2);
658
659  // Get partials with respect to x and y of the grid function
660  disc.GetDerivative(1,0,d_x);
661  disc.GetDerivative(1,1,d_y);
662
663  // Compute the squared magnitude of the gradient
664  for (int i = 0; i < size.Size(); i++)
665  {
666  size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
667  }
668  const double max = size.Max();
669  double max_all;
670  MPI_Allreduce(&max, &max_all, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
671
672  for (int i = 0; i < d_x.Size(); i++)
673  {
674  d_x(i) = std::abs(d_x(i));
675  d_y(i) = std::abs(d_y(i));
676  }
677  const double eps = 0.01;
678  const double aspr_ratio = 20.0;
679  const double size_ratio = 40.0;
680
681  for (int i = 0; i < size.Size(); i++)
682  {
683  size(i) = (size(i)/max_all);
684  aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
685  aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
686  if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
687  if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
688  }
689  Vector vals;
690  const int NE = pmesh->GetNE();
691  double volume = 0.0, volume_ind = 0.0;
692
693  for (int i = 0; i < NE; i++)
694  {
696  const IntegrationRule &ir =
697  IntRules.Get(pmesh->GetElementBaseGeometry(i), Tr->OrderJ());
698  size.GetValues(i, ir, vals);
699  for (int j = 0; j < ir.GetNPoints(); j++)
700  {
701  const IntegrationPoint &ip = ir.IntPoint(j);
702  Tr->SetIntPoint(&ip);
703  volume += ip.weight * Tr->Weight();
704  volume_ind += vals(j) * ip.weight * Tr->Weight();
705  }
706  }
707  double volume_all, volume_ind_all;
708  MPI_Allreduce(&volume, &volume_all, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
709  MPI_Allreduce(&volume_ind, &volume_ind_all, 1, MPI_DOUBLE, MPI_SUM,
710  MPI_COMM_WORLD);
711  const int NE_ALL = pmesh->GetGlobalNE();
712
713  const double avg_zone_size = volume_all / NE_ALL;
714
715  const double small_avg_ratio =
716  (volume_ind_all + (volume_all - volume_ind_all) / size_ratio)
717  / volume_all;
718
719  const double small_zone_size = small_avg_ratio * avg_zone_size;
720  const double big_zone_size = size_ratio * small_zone_size;
721
722  for (int i = 0; i < size.Size(); i++)
723  {
724  const double val = size(i);
725  const double a = (big_zone_size - small_zone_size) / small_zone_size;
726  size(i) = big_zone_size / (1.0+a*val);
727  }
728
729  DiffuseField(size, 2);
730  DiffuseField(aspr, 2);
731
732  tc->SetParDiscreteTargetSize(size);
734  target_c = tc;
735  break;
736  }
737  case 7: // Discrete aspect ratio 3D
738  {
742  {
744  }
745  else
746  {
747 #ifdef MFEM_USE_GSLIB
749 #else
750  MFEM_ABORT("MFEM is not built with GSLIB.");
751 #endif
752  }
754  aspr3d.ProjectCoefficient(fd_aspr3d);
756  target_c = tc;
757  break;
758  }
759  case 8: // shape/size + orientation 2D
760  {
764  {
766  }
767  else
768  {
769 #ifdef MFEM_USE_GSLIB
771 #else
772  MFEM_ABORT("MFEM is not built with GSLIB.");
773 #endif
774  }
775
776  ConstantCoefficient size_coeff(0.1*0.1);
777  size.ProjectCoefficient(size_coeff);
778  tc->SetParDiscreteTargetSize(size);
779
781  ori.ProjectCoefficient(ori_coeff);
783  target_c = tc;
784  break;
785  }
786  // Targets used for hr-adaptivity tests.
787  case 9: // size target in an annular region.
788  case 10: // size+aspect-ratio in an annular region.
789  case 11: // size+aspect-ratio target for a rotate sine wave
790  {
793  hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
795  target_c = tc;
796  break;
797  }
798  default:
799  if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
800  return 3;
801  }
802  if (target_c == NULL)
803  {
804  target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
805  }
806  target_c->SetNodes(x0);
807
808  // Automatically balanced gamma in composite metrics.
809  auto metric_combo = dynamic_cast<TMOP_Combo_QualityMetric *>(metric);
810  if (metric_combo && bal_expl_combo)
811  {
812  Vector bal_weights;
813  metric_combo->ComputeBalancedWeights(x, *target_c, bal_weights);
814  metric_combo->SetWeights(bal_weights);
815  }
816
817  TMOP_QualityMetric *metric_to_use = barrier_type > 0 || worst_case_type > 0
818  ? untangler_metric
819  : metric;
820  TMOP_Integrator *tmop_integ = new TMOP_Integrator(metric_to_use, target_c,
821  h_metric);
822  if (barrier_type > 0 || worst_case_type > 0)
823  {
824  tmop_integ->ComputeUntangleMetricQuantiles(x, *pfespace);
825  }
826
827  // Finite differences for computations of derivatives.
828  if (fdscheme)
829  {
830  MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
831  tmop_integ->EnableFiniteDifferences(x);
832  }
833  tmop_integ->SetExactActionFlag(exactaction);
834
835  // Setup the quadrature rules for the TMOP integrator.
836  IntegrationRules *irules = NULL;
838  {
839  case 1: irules = &IntRulesLo; break;
840  case 2: irules = &IntRules; break;
841  case 3: irules = &IntRulesCU; break;
842  default:
843  if (myid == 0) { cout << "Unknown quad_type: " << quad_type << endl; }
844  return 3;
845  }
847  if (myid == 0 && dim == 2)
848  {
849  cout << "Triangle quadrature points: "
852  << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
853  }
854  if (myid == 0 && dim == 3)
855  {
856  cout << "Tetrahedron quadrature points: "
858  << "\nHexahedron quadrature points: "
860  << "\nPrism quadrature points: "
861  << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
862  }
863
864  // Limit the node movement.
865  // The limiting distances can be given by a general function of space.
866  ParFiniteElementSpace dist_pfespace(pmesh, fec); // scalar space
867  ParGridFunction dist(&dist_pfespace);
868  dist = 1.0;
869  // The small_phys_size is relevant only with proper normalization.
870  if (normalization) { dist = small_phys_size; }
871  ConstantCoefficient lim_coeff(lim_const);
872  if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
873
879  {
880  MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
881
884
886  else if (adapt_eval == 1)
887  {
888 #ifdef MFEM_USE_GSLIB
890 #else
891  MFEM_ABORT("MFEM is not built with GSLIB support!");
892 #endif
893  }
894  else { MFEM_ABORT("Bad interpolation option."); }
895
898  if (visualization)
899  {
900  socketstream vis1;
901  common::VisualizeField(vis1, "localhost", 19916, adapt_lim_gf0, "Zeta 0",
902  300, 600, 300, 300);
903  }
904  }
905
906  // Surface fitting.
907  L2_FECollection mat_coll(0, dim);
908  H1_FECollection surf_fit_fec(mesh_poly_deg, dim);
909  ParFiniteElementSpace surf_fit_fes(pmesh, &surf_fit_fec);
910  ParFiniteElementSpace mat_fes(pmesh, &mat_coll);
911  ParGridFunction mat(&mat_fes);
912  ParGridFunction surf_fit_mat_gf(&surf_fit_fes);
913  ParGridFunction surf_fit_gf0(&surf_fit_fes);
914  Array<bool> surf_fit_marker(surf_fit_gf0.Size());
915  ConstantCoefficient surf_fit_coeff(surface_fit_const);
917  if (surface_fit_const > 0.0)
918  {
920  "Surface fitting with HR is not implemented yet.");
921  MFEM_VERIFY(pa == false,
922  "Surface fitting with PA is not implemented yet.");
923
925  surf_fit_gf0.ProjectCoefficient(ls_coeff);
926
927  for (int i = 0; i < pmesh->GetNE(); i++)
928  {
929  mat(i) = material_id(i, surf_fit_gf0);
930  pmesh->SetAttribute(i, static_cast<int>(mat(i) + 1));
931  }
932
933  GridFunctionCoefficient coeff_mat(&mat);
934  surf_fit_mat_gf.ProjectDiscCoefficient(coeff_mat, GridFunction::ARITHMETIC);
935  for (int j = 0; j < surf_fit_marker.Size(); j++)
936  {
937  if (surf_fit_mat_gf(j) > 0.1 && surf_fit_mat_gf(j) < 0.9)
938  {
939  surf_fit_marker[j] = true;
940  surf_fit_mat_gf(j) = 1.0;
941  }
942  else
943  {
944  surf_fit_marker[j] = false;
945  surf_fit_mat_gf(j) = 0.0;
946  }
947  }
948
950  else if (adapt_eval == 1)
951  {
952 #ifdef MFEM_USE_GSLIB
954 #else
955  MFEM_ABORT("MFEM is not built with GSLIB support!");
956 #endif
957  }
958  else { MFEM_ABORT("Bad interpolation option."); }
959
960  tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker, surf_fit_coeff,
962  if (visualization)
963  {
964  socketstream vis1, vis2, vis3;
965  common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0, "Level Set 0",
966  300, 600, 300, 300);
967  common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
968  600, 600, 300, 300);
969  common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
970  "Dofs to Move",
971  900, 600, 300, 300);
972  }
973  }
974
975  // Has to be after the enabling of the limiting / alignment, as it computes
976  // normalization factors for these terms as well.
977  if (normalization) { tmop_integ->ParEnableNormalization(x0); }
978
979  // 13. Setup the final NonlinearForm (which defines the integral of interest,
980  // its first and second derivatives). Here we can use a combination of
981  // metrics, i.e., optimize the sum of two integrals, where both are
982  // scaled by used-defined space-dependent weights. Note that there are
983  // no command-line options for the weights and the type of the second
984  // metric; one should update those in the code.
985  ParNonlinearForm a(pfespace);
986  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
987  ConstantCoefficient *metric_coeff1 = NULL;
988  TMOP_QualityMetric *metric2 = NULL;
989  TargetConstructor *target_c2 = NULL;
990  FunctionCoefficient metric_coeff2(weight_fun);
991
992  // Explicit combination of metrics.
993  if (combomet > 0)
994  {
995  // First metric.
996  metric_coeff1 = new ConstantCoefficient(1.0);
997  tmop_integ->SetCoefficient(*metric_coeff1);
998
999  // Second metric.
1000  if (dim == 2) { metric2 = new TMOP_Metric_077; }
1001  else { metric2 = new TMOP_Metric_315; }
1002  TMOP_Integrator *tmop_integ2 = NULL;
1003  if (combomet == 1)
1004  {
1005  target_c2 = new TargetConstructor(
1007  target_c2->SetVolumeScale(0.01);
1008  target_c2->SetNodes(x0);
1009  tmop_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
1010  tmop_integ2->SetCoefficient(metric_coeff2);
1011  }
1012  else { tmop_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
1014  if (fdscheme) { tmop_integ2->EnableFiniteDifferences(x); }
1015  tmop_integ2->SetExactActionFlag(exactaction);
1016
1020  if (normalization) { combo->ParEnableNormalization(x0); }
1021  if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
1022
1024  }
1025  else
1026  {
1028  }
1029
1030  if (pa) { a.Setup(); }
1031
1032  // Compute the minimum det(J) of the starting mesh.
1033  min_detJ = infinity();
1034  const int NE = pmesh->GetNE();
1035  for (int i = 0; i < NE; i++)
1036  {
1037  const IntegrationRule &ir =
1039  ElementTransformation *transf = pmesh->GetElementTransformation(i);
1040  for (int j = 0; j < ir.GetNPoints(); j++)
1041  {
1042  transf->SetIntPoint(&ir.IntPoint(j));
1043  min_detJ = min(min_detJ, transf->Jacobian().Det());
1044  }
1045  }
1046  double minJ0;
1047  MPI_Allreduce(&min_detJ, &minJ0, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
1048  min_detJ = minJ0;
1049  if (myid == 0)
1050  { cout << "Minimum det(J) of the original mesh is " << min_detJ << endl; }
1051
1052  if (min_detJ < 0.0 && barrier_type == 0
1053  && metric_id != 22 && metric_id != 211 && metric_id != 252
1054  && metric_id != 311 && metric_id != 313 && metric_id != 352)
1055  {
1056  MFEM_ABORT("The input mesh is inverted! Try an untangling metric.");
1057  }
1058  if (min_detJ < 0.0)
1059  {
1060  MFEM_VERIFY(target_t == TargetConstructor::IDEAL_SHAPE_UNIT_SIZE,
1061  "Untangling is supported only for ideal targets.");
1062
1063  const DenseMatrix &Wideal =
1065  min_detJ /= Wideal.Det();
1066
1067  double h0min = h0.Min(), h0min_all;
1068  MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
1069  // Slightly below minJ0 to avoid div by 0.
1070  min_detJ -= 0.01 * h0min_all;
1071  }
1072
1073  // For HR tests, the energy is normalized by the number of elements.
1074  const double init_energy = a.GetParGridFunctionEnergy(x) /
1075  (hradaptivity ? pmesh->GetGlobalNE() : 1);
1076  double init_metric_energy = init_energy;
1077  if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
1078  {
1079  lim_coeff.constant = 0.0;
1081  surf_fit_coeff.constant = 0.0;
1082  init_metric_energy = a.GetParGridFunctionEnergy(x) /
1083  (hradaptivity ? pmesh->GetGlobalNE() : 1);
1084  lim_coeff.constant = lim_const;
1086  surf_fit_coeff.constant = surface_fit_const;
1087  }
1088
1089  // Visualize the starting mesh and metric values.
1090  // Note that for combinations of metrics, this only shows the first metric.
1091  if (visualization)
1092  {
1093  char title[] = "Initial metric values";
1094  vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
1095  }
1096
1097  // 14. Fix all boundary nodes, or fix only a given component depending on the
1098  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
1099  // fixed x/y/z components of the node. Attribute dim+1 corresponds to
1100  // an entirely fixed node.
1101  if (move_bnd == false)
1102  {
1103  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
1104  ess_bdr = 1;
1105  a.SetEssentialBC(ess_bdr);
1106  }
1107  else
1108  {
1109  int n = 0;
1110  for (int i = 0; i < pmesh->GetNBE(); i++)
1111  {
1112  const int nd = pfespace->GetBE(i)->GetDof();
1113  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
1114  MFEM_VERIFY(!(dim == 2 && attr == 3),
1115  "Boundary attribute 3 must be used only for 3D meshes. "
1116  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
1118  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
1119  if (attr == 4) { n += nd * dim; }
1120  }
1121  Array<int> ess_vdofs(n);
1122  n = 0;
1123  for (int i = 0; i < pmesh->GetNBE(); i++)
1124  {
1125  const int nd = pfespace->GetBE(i)->GetDof();
1126  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
1127  pfespace->GetBdrElementVDofs(i, vdofs);
1128  if (attr == 1) // Fix x components.
1129  {
1130  for (int j = 0; j < nd; j++)
1131  { ess_vdofs[n++] = vdofs[j]; }
1132  }
1133  else if (attr == 2) // Fix y components.
1134  {
1135  for (int j = 0; j < nd; j++)
1136  { ess_vdofs[n++] = vdofs[j+nd]; }
1137  }
1138  else if (attr == 3) // Fix z components.
1139  {
1140  for (int j = 0; j < nd; j++)
1141  { ess_vdofs[n++] = vdofs[j+2*nd]; }
1142  }
1143  else if (attr == 4) // Fix all components.
1144  {
1145  for (int j = 0; j < vdofs.Size(); j++)
1146  { ess_vdofs[n++] = vdofs[j]; }
1147  }
1148  }
1149  a.SetEssentialVDofs(ess_vdofs);
1150  }
1151
1152  // 15. As we use the Newton method to solve the resulting nonlinear system,
1153  // here we setup the linear solver for the system's Jacobian.
1154  Solver *S = NULL, *S_prec = NULL;
1155  const double linsol_rtol = 1e-12;
1156  if (lin_solver == 0)
1157  {
1158  S = new DSmoother(1, 1.0, max_lin_iter);
1159  }
1160  else if (lin_solver == 1)
1161  {
1162  CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
1163  cg->SetMaxIter(max_lin_iter);
1164  cg->SetRelTol(linsol_rtol);
1165  cg->SetAbsTol(0.0);
1166  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
1167  S = cg;
1168  }
1169  else
1170  {
1171  MINRESSolver *minres = new MINRESSolver(MPI_COMM_WORLD);
1172  minres->SetMaxIter(max_lin_iter);
1173  minres->SetRelTol(linsol_rtol);
1174  minres->SetAbsTol(0.0);
1175  if (verbosity_level > 2) { minres->SetPrintLevel(1); }
1176  else { minres->SetPrintLevel(verbosity_level == 2 ? 3 : -1); }
1177  if (lin_solver == 3 || lin_solver == 4)
1178  {
1179  if (pa)
1180  {
1181  MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
1182  auto js = new OperatorJacobiSmoother;
1183  js->SetPositiveDiagonal(true);
1184  S_prec = js;
1185  }
1186  else
1187  {
1188  auto hs = new HypreSmoother;
1189  hs->SetType((lin_solver == 3) ? HypreSmoother::Jacobi
1190  /* */ : HypreSmoother::l1Jacobi, 1);
1191  hs->SetPositiveDiagonal(true);
1192  S_prec = hs;
1193  }
1194  minres->SetPreconditioner(*S_prec);
1195  }
1196  S = minres;
1197  }
1198
1199  // Perform the nonlinear optimization.
1200  const IntegrationRule &ir =
1202  TMOPNewtonSolver solver(pfespace->GetComm(), ir, solver_type);
1204  if (surface_fit_threshold > 0)
1205  {
1206  solver.SetTerminationWithMaxSurfaceFittingError(surface_fit_threshold);
1207  }
1208  // Provide all integration rules in case of a mixed mesh.
1210  if (solver_type == 0)
1211  {
1212  // Specify linear solver when we use a Newton-based solver.
1213  solver.SetPreconditioner(*S);
1214  }
1215  // For untangling, the solver will update the min det(T) values.
1216  solver.SetMinDetPtr(&min_detJ);
1217  solver.SetMaxIter(solver_iter);
1218  solver.SetRelTol(solver_rtol);
1219  solver.SetAbsTol(0.0);
1220  if (solver_art_type > 0)
1221  {
1223  }
1224  solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
1225
1227  // If hr-adaptivity is disabled, r-adaptivity is done once using the
1228  // TMOPNewtonSolver.
1229  // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
1231  // The solver terminates if an h-adaptivity iteration does not modify
1232  // any element in the mesh.
1233  TMOPHRSolver hr_solver(*pmesh, a, solver,
1235  mesh_poly_deg, h_metric_id,
1236  n_hr_iter, n_h_iter);
1239  {
1242  }
1243  hr_solver.Mult();
1244
1245  // 16. Save the optimized mesh to a file. This output can be viewed later
1246  // using GLVis: "glvis -m optimized -np num_mpi_tasks".
1247  {
1248  ostringstream mesh_name;
1249  mesh_name << "optimized.mesh";
1250  ofstream mesh_ofs(mesh_name.str().c_str());
1251  mesh_ofs.precision(8);
1252  pmesh->PrintAsOne(mesh_ofs);
1253  }
1254
1255  // Report the final energy of the functional.
1256  const double fin_energy = a.GetParGridFunctionEnergy(x) /
1257  (hradaptivity ? pmesh->GetGlobalNE() : 1);
1258  double fin_metric_energy = fin_energy;
1259  if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
1260  {
1261  lim_coeff.constant = 0.0;
1263  surf_fit_coeff.constant = 0.0;
1264  fin_metric_energy = a.GetParGridFunctionEnergy(x) /
1265  (hradaptivity ? pmesh->GetGlobalNE() : 1);
1266  lim_coeff.constant = lim_const;
1268  surf_fit_coeff.constant = surface_fit_const;
1269  }
1270  if (myid == 0)
1271  {
1272  std::cout << std::scientific << std::setprecision(4);
1273  cout << "Initial strain energy: " << init_energy
1274  << " = metrics: " << init_metric_energy
1275  << " + extra terms: " << init_energy - init_metric_energy << endl;
1276  cout << " Final strain energy: " << fin_energy
1277  << " = metrics: " << fin_metric_energy
1278  << " + extra terms: " << fin_energy - fin_metric_energy << endl;
1279  cout << "The strain energy decreased by: "
1280  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
1281  }
1282
1283  // Visualize the final mesh and metric values.
1284  if (visualization)
1285  {
1286  char title[] = "Final metric values";
1287  vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
1288  }
1289
1290  if (adapt_lim_const > 0.0 && visualization)
1291  {
1292  socketstream vis0;
1293  common::VisualizeField(vis0, "localhost", 19916, adapt_lim_gf0, "Xi 0",
1294  600, 600, 300, 300);
1295  }
1296
1297  // Visualize fitting surfaces and report fitting errors.
1298  if (surface_fit_const > 0.0)
1299  {
1300  if (visualization)
1301  {
1302  socketstream vis2, vis3;
1303  common::VisualizeField(vis2, "localhost", 19916, mat,
1304  "Materials", 600, 900, 300, 300);
1305  common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
1306  "Surface dof", 900, 900, 300, 300);
1307  }
1308  double err_avg, err_max;
1309  tmop_integ->GetSurfaceFittingErrors(err_avg, err_max);
1310  if (myid == 0)
1311  {
1312  std::cout << "Avg fitting error: " << err_avg << std::endl
1313  << "Max fitting error: " << err_max << std::endl;
1314  }
1315  }
1316
1317  // Visualize the mesh displacement.
1318  if (visualization)
1319  {
1320  x0 -= x;
1321  socketstream sock;
1322  if (myid == 0)
1323  {
1324  sock.open("localhost", 19916);
1325  sock << "solution\n";
1326  }
1327  pmesh->PrintAsOne(sock);
1328  x0.SaveAsOne(sock);
1329  if (myid == 0)
1330  {
1331  sock << "window_title 'Displacements'\n"
1332  << "window_geometry "
1333  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
1334  << "keys jRmclA" << endl;
1335  }
1336  }
1337
1338  delete S;
1339  delete S_prec;
1340  delete target_c2;
1341  delete metric2;
1342  delete metric_coeff1;
1345  delete target_c;
1348  delete h_metric;
1349  delete metric;
1350  delete untangler_metric;
1351  delete pfespace;
1352  delete fec;
1353  delete pmesh;
1354
1355  return 0;
1356 }
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
Definition: tmop_tools.cpp:899
void ComputeUntangleMetricQuantiles(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:4238
void discrete_aspr_3d(const Vector &x, Vector &v)
virtual void SetParDiscreteTargetSize(const ParGridFunction &tspec_)
Definition: tmop.cpp:1785
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:80
Definition: solvers.hpp:493
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Definition: ex25.cpp:148
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:956
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:4090
double discrete_size_2d(const Vector &x)
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:771
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:339
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:150
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:893
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1108
int Dimension() const
Definition: mesh.hpp:1047
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:899
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:652
Definition: tmop.hpp:1563
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
double surface_level_set(const Vector &x)
double Det() const
Definition: densemat.cpp:488
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
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 & invex).
Definition: tmop.hpp:631
void GetSurfaceFittingErrors(double &err_avg, double &err_max)
Definition: tmop.cpp:2882
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:1320
Parallel non-linear operator on the true dofs.
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
Restriction of the node positions to certain regions.
Definition: tmop.cpp:2811
3D non-barrier metric without a type.
Definition: tmop.hpp:734
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop.hpp:1896
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:1050
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:813
virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
Project a discontinuous vector coefficient as a grid function on a continuous finite element space...
Definition: pgridfunc.cpp:546
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
Definition: tmop.hpp:2051
STL namespace.
MINRES method.
Definition: solvers.hpp:603
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:786
virtual int OrderJ() const =0
Return the order of the elements of the Jacobian of the transformation.
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:338
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:792
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:939
Geometry Geometries
Definition: fe.cpp:49
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:4278
2D Shifted barrier form of shape metric (mu_2).
Definition: tmop.hpp:387
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:460
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:673
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
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:250
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:968
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:616
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:2786
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:1908
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9535
3D barrier metric without a type.
Definition: tmop.hpp:752
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
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:1314
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:915
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:752
int main(int argc, char *argv[])
virtual void SetParDiscreteTargetAspectRatio(const ParGridFunction &tspec_)
Definition: tmop.cpp:1805
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
Parallel smoothers in hypre.
Definition: hypre.hpp:971
static void Init()
Singleton creation with Mpi::Init();.
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:4393
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:558
Definition: tmop_amr.hpp:253
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
Definition: tmop_tools.hpp:225
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:1032
virtual DofTransformation * GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: pfespace.cpp:471
A general vector function coefficient.
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:834
void SetAbsTol(double atol)
Definition: solvers.hpp:200
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
Definition: tmop_amr.hpp:252
void SetRelTol(double rtol)
Definition: solvers.hpp:199
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:535
double weight_fun(const Vector &x)
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 PrintAsOne(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4983
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
Adds a new TMOP_Integrator to the combination.
Definition: tmop.hpp:2081
double material_indicator_2d(const Vector &x)
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:356
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
double a
Definition: lissajous.cpp:41
double discrete_size_3d(const Vector &x)
Class for integration point with weight.
Definition: intrules.hpp:25
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
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+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:1014
3D Shape (S) metric, untangling version of 303.
Definition: tmop.hpp:713
int dim
Definition: ex24.cpp:53
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:623
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:46
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
A general function coefficient.
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1547
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:251
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:874
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition: pmesh.cpp:2076
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:1238
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:2846
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1085
virtual void SetParDiscreteTargetOrientation(const ParGridFunction &tspec_)
Definition: tmop.cpp:1815
Base class for solvers.
Definition: operator.hpp:682
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:5345
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1234
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:1639
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:853
int material_id(int el_id, const GridFunction &g)
Class for parallel meshes.
Definition: pmesh.hpp:32
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:3102
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:3435
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:543
Definition: vector.hpp:468
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:3936
double GetElementVolume(int i)
Definition: mesh.cpp:112
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:612
2D non-barrier metric without a type.
Definition: tmop.hpp:219
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1638
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:372