MFEM  v4.6.0 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 94 -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 NC mesh;
52 // mpirun -np 4 pmesh-optimizer -m amr-quad-q2.mesh -o 2 -rs 2 -mid 94 -tid 5 -ni 50 -qo 4 -nor
53 // Adapted discrete size 3D with PA:
54 // mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa
55 // Adapted discrete size 3D with PA on device (requires CUDA):
56 // * 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
57 // Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
58 // 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
60 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
61 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
63 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -nor
64 // Adapted discrete aspect ratio (3D):
65 // mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -bnd -qt 1 -qo 8
66 //
68 // mpirun -np 4 pmesh-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 // 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
71 // Adaptive limiting through FD (requires GSLIB):
72 // * 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
73 //
75 // 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
76 // (requires CUDA):
77 // * 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
78 // Blade shape with FD-based solver:
79 // mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
81 // mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
82 // ICF shape and equal size:
83 // mpirun -np 4 pmesh-optimizer -o 3 -mid 80 -bec -tid 2 -ni 25 -ls 3 -art 2 -qo 5
84 // ICF shape and initial size:
85 // mpirun -np 4 pmesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
86 // ICF shape:
87 // mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
88 // ICF limited shape:
89 // mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8 -lc 10
90 // ICF combo shape + size (rings, slow convergence):
91 // mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 1000 -bnd -qt 1 -qo 8 -cmb 1
92 // Mixed tet / cube / hex mesh with limiting:
93 // 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
94 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
95 // * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -mid 303 -tid 1 -ni 20 -li 500 -fix-bnd
96 // 2D non-conforming shape and equal size:
97 // 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
98 //
99 // 2D untangling:
100 // mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
101 // 2D untangling with shifted barrier metric:
102 // mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
103 // 3D untangling (the mesh is in the mfem/data GitHub repository):
104 // * 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
105 // Shape optimization for a Kershaw transformed mesh using partial assembly:
106 // Mesh for Kershaw transformation must be a Cartesian mesh with nx % 6 = ny % 2 = nz % 2 = 0.
107 // Kershaw transformation can be imposed using the transformation ('t') feature in the mesh-explorer miniapp.
108 // * 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
109
110 #include "mfem.hpp"
111 #include "../common/mfem-common.hpp"
112 #include <iostream>
113 #include <fstream>
114 #include "mesh-optimizer.hpp"
115
116 using namespace mfem;
117 using namespace std;
118
119 int main (int argc, char *argv[])
120 {
121  // 0. Initialize MPI and HYPRE.
122  Mpi::Init(argc, argv);
123  int myid = Mpi::WorldRank();
124  Hypre::Init();
125
126  // 1. Set the method's default parameters.
127  const char *mesh_file = "icf.mesh";
128  int mesh_poly_deg = 1;
129  int rs_levels = 0;
130  int rp_levels = 0;
131  double jitter = 0.0;
132  int metric_id = 1;
133  int target_id = 1;
134  double lim_const = 0.0;
138  int solver_type = 0;
139  int solver_iter = 20;
140  double solver_rtol = 1e-10;
141  int solver_art_type = 0;
142  int lin_solver = 2;
143  int max_lin_iter = 100;
144  bool move_bnd = true;
145  int combomet = 0;
146  bool bal_expl_combo = false;
148  int h_metric_id = -1;
149  bool normalization = false;
150  bool visualization = true;
151  int verbosity_level = 0;
152  bool fdscheme = false;
154  bool exactaction = false;
155  bool integ_over_targ = true;
156  const char *devopt = "cpu";
157  bool pa = false;
158  int n_hr_iter = 5;
159  int n_h_iter = 1;
160  int mesh_node_ordering = 0;
161  int barrier_type = 0;
162  int worst_case_type = 0;
163
164  // 2. Parse command-line options.
165  OptionsParser args(argc, argv);
167  "Mesh file to use.");
169  "Polynomial degree of mesh finite element space.");
171  "Number of times to refine the mesh uniformly in serial.");
173  "Number of times to refine the mesh uniformly in parallel.");
175  "Random perturbation scaling factor.");
177  "Mesh optimization metric:\n\t"
178  "T-metrics\n\t"
179  "1 : |T|^2 -- 2D no type\n\t"
180  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
181  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
182  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
183  "14 : |T-I|^2 -- 2D shape+size+orientation\n\t"
184  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
185  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
186  "55 : (tau-1)^2 -- 2D size\n\t"
187  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
188  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
189  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
190  "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t"
191  "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t"
192  "90 : balanced combo mu_50 & mu_77 -- 2D shape+size\n\t"
193  "94 : balanced combo mu_2 & mu_56 -- 2D shape+size\n\t"
194  "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t"
195  // "211: (tau-1)^2-tau+sqrt(tau^2+eps) -- 2D untangling\n\t"
196  // "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
197  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
198  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
199  "303: (|T|^2)/3/tau^(2/3)-1 -- 3D shape\n\t"
200  "304: (|T|^3)/3^{3/2}/tau-1 -- 3D shape\n\t"
201  // "311: (tau-1)^2-tau+sqrt(tau^2+eps)-- 3D untangling\n\t"
202  "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t"
203  "315: (tau-1)^2 -- 3D no type\n\t"
204  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D no type\n\t"
205  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
206  "322: |T-adjT^-t|^2 -- 3D shape+size\n\t"
207  "323: |J|^3-3sqrt(3)ln(det(J))-3sqrt(3) -- 3D shape+size\n\t"
208  "328: balanced combo mu_301 & mu_316 -- 3D shape+size\n\t"
209  "332: (1-gamma) mu_302 + gamma mu_315 -- 3D shape+size\n\t"
210  "333: (1-gamma) mu_302 + gamma mu_316 -- 3D shape+size\n\t"
211  "334: (1-gamma) mu_303 + gamma mu_316 -- 3D shape+size\n\t"
212  "328: balanced combo mu_302 & mu_318 -- 3D shape+size\n\t"
213  "347: (1-gamma) mu_304 + gamma mu_316 -- 3D shape+size\n\t"
214  // "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling\n\t"
215  "360: (|T|^3)/3^{3/2}-tau -- 3D shape\n\t"
216  "A-metrics\n\t"
217  "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t"
218  "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t"
219  "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t"
220  "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t"
221  );
223  "Target (ideal element) type:\n\t"
224  "1: Ideal shape, unit size\n\t"
225  "2: Ideal shape, equal size\n\t"
226  "3: Ideal shape, initial size\n\t"
227  "4: Given full analytic Jacobian (in physical space)\n\t"
228  "5: Ideal shape, given size (in physical space)");
229  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
234  "1: Gauss-Lobatto\n\t"
235  "2: Gauss-Legendre\n\t"
236  "3: Closed uniform points");
238  "Order of the quadrature rule.");
240  " Type of solver: (default) 0: Newton, 1: LBFGS");
242  "Maximum number of Newton iterations.");
244  "Relative tolerance for the Newton solver.");
246  "Type of adaptive relative linear solver tolerance:\n\t"
247  "0: None (default)\n\t"
248  "1: Eisenstat-Walker type 1\n\t"
249  "2: Eisenstat-Walker type 2");
251  "Linear solver:\n\t"
252  "0: l1-Jacobi\n\t"
253  "1: CG\n\t"
254  "2: MINRES\n\t"
255  "3: MINRES + Jacobi preconditioner\n\t"
256  "4: MINRES + l1-Jacobi preconditioner");
258  "Maximum number of iterations in the linear solve.");
260  "--fix-boundary",
261  "Enable motion along horizontal and vertical boundaries.");
263  "Combination of metrics options:\n\t"
264  "0: Use single metric\n\t"
265  "1: Shape + space-dependent size given analytically\n\t"
266  "2: Shape + adapted size given discretely; shared target");
268  "-no-bec", "--balance-explicit-combo",
269  "Automatic balancing of explicit combo metrics.");
274  "Same options as metric_id. Used to determine refinement"
275  " type for each element if h-adaptivity is enabled.");
277  "--no-normalization",
278  "Make all terms in the optimization functional unitless.");
280  "-no-fd", "--no-fd-approx",
281  "Enable finite difference based derivative computations.");
283  "-no-ex", "--no-exact-action",
284  "Enable exact action of TMOP_Integrator.");
286  "-ir", "--integrate-reference",
287  "Integrate over target (-it) or reference (-ir) element.");
289  "--no-visualization",
290  "Enable or disable GLVis visualization.");
292  "Verbosity level for the involved iterative solvers:\n\t"
293  "0: no output\n\t"
294  "1: Newton iterations\n\t"
295  "2: Newton iterations + linear solver summaries\n\t"
296  "3: newton iterations + linear solver iterations");
298  "0 - Advection based (DEFAULT), 1 - GSLIB.");
300  "Device configuration string, see Device::Configure().");
302  "--no-partial-assembly", "Enable Partial Assembly.");
307  "iteration.");
309  "Ordering of mesh nodes."
310  "0 (default): byNodes, 1: byVDIM");
312  "0 - None,"
313  "1 - Shifted Barrier,"
314  "2 - Pseudo Barrier.");
316  "0 - None,"
317  "1 - Beta,"
318  "2 - PMean.");
319
320  args.Parse();
321  if (!args.Good())
322  {
323  if (myid == 0) { args.PrintUsage(cout); }
324  return 1;
325  }
326  if (myid == 0) { args.PrintOptions(cout); }
327  if (h_metric_id < 0) { h_metric_id = metric_id; }
328
330  {
331  MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
332  " supported on cpus.");
333  }
334  Device device(devopt);
335  if (myid == 0) { device.Print();}
336
337  // 3. Initialize and refine the starting mesh.
338  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
339  for (int lev = 0; lev < rs_levels; lev++)
340  {
341  mesh->UniformRefinement();
342  }
343  const int dim = mesh->Dimension();
344
345  if (hradaptivity) { mesh->EnsureNCMesh(); }
346  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
347  delete mesh;
348  for (int lev = 0; lev < rp_levels; lev++)
349  {
350  pmesh->UniformRefinement();
351  }
352
353  // 4. Define a finite element space on the mesh. Here we use vector finite
354  // elements which are tensor products of quadratic finite elements. The
355  // number of components in the vector finite element space is specified by
356  // the last parameter of the FiniteElementSpace constructor.
358  if (mesh_poly_deg <= 0)
359  {
361  mesh_poly_deg = 2;
362  }
363  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
364  ParFiniteElementSpace *pfespace = new ParFiniteElementSpace(pmesh, fec, dim,
365  mesh_node_ordering);
366
367  // 5. Make the mesh curved based on the above finite element space. This
368  // means that we define the mesh elements through a fespace-based
369  // transformation of the reference element.
370  pmesh->SetNodalFESpace(pfespace);
371
372  // 7. Get the mesh nodes (vertices and other degrees of freedom in the finite
373  // element space) as a finite element grid function in fespace. Note that
374  // changing x automatically changes the shapes of the mesh elements.
375  ParGridFunction x(pfespace);
376  pmesh->SetNodalGridFunction(&x);
377
378  // 8. Define a vector representing the minimal local mesh size in the mesh
379  // nodes. We index the nodes using the scalar version of the degrees of
380  // freedom in pfespace. Note: this is partition-dependent.
381  //
382  // In addition, compute average mesh size and total volume.
383  Vector h0(pfespace->GetNDofs());
384  h0 = infinity();
385  double vol_loc = 0.0;
386  Array<int> dofs;
387  for (int i = 0; i < pmesh->GetNE(); i++)
388  {
389  // Get the local scalar element degrees of freedom in dofs.
390  pfespace->GetElementDofs(i, dofs);
391  // Adjust the value of h0 in dofs based on the local mesh size.
392  const double hi = pmesh->GetElementSize(i);
393  for (int j = 0; j < dofs.Size(); j++)
394  {
395  h0(dofs[j]) = min(h0(dofs[j]), hi);
396  }
397  vol_loc += pmesh->GetElementVolume(i);
398  }
399  double vol_glb;
400  MPI_Allreduce(&vol_loc, &vol_glb, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
401  const double small_phys_size = pow(vol_glb, 1.0 / dim) / 100.0;
402
403  // 9. Add a random perturbation to the nodes in the interior of the domain.
404  // We define a random grid function of fespace and make sure that it is
405  // zero on the boundary and its values are locally of the order of h0.
406  // The latter is based on the DofToVDof() method which maps the scalar to
407  // the vector degrees of freedom in pfespace.
408  ParGridFunction rdm(pfespace);
409  rdm.Randomize();
410  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
411  rdm *= jitter;
413  // Scale the random values to be of order of the local mesh size.
414  for (int i = 0; i < pfespace->GetNDofs(); i++)
415  {
416  for (int d = 0; d < dim; d++)
417  {
418  rdm(pfespace->DofToVDof(i,d)) *= h0(i);
419  }
420  }
421  Array<int> vdofs;
422  for (int i = 0; i < pfespace->GetNBE(); i++)
423  {
424  // Get the vector degrees of freedom in the boundary element.
425  pfespace->GetBdrElementVDofs(i, vdofs);
426  // Set the boundary values to zero.
427  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
428  }
429  x -= rdm;
430  // Set the perturbation of all nodes from the true nodes.
431  x.SetTrueVector();
432  x.SetFromTrueVector();
433
434  // 10. Save the starting (prior to the optimization) mesh to a file. This
435  // output can be viewed later using GLVis: "glvis -m perturbed -np
437  {
438  ostringstream mesh_name;
439  mesh_name << "perturbed.mesh";
440  ofstream mesh_ofs(mesh_name.str().c_str());
441  mesh_ofs.precision(8);
442  pmesh->PrintAsOne(mesh_ofs);
443  }
444
445  // 11. Store the starting (prior to the optimization) positions.
446  ParGridFunction x0(pfespace);
447  x0 = x;
448
449  // 12. Form the integrator that uses the chosen metric and target.
450  double min_detJ = -0.1;
451  TMOP_QualityMetric *metric = NULL;
452  switch (metric_id)
453  {
454  // T-metrics
455  case 1: metric = new TMOP_Metric_001; break;
456  case 2: metric = new TMOP_Metric_002; break;
457  case 4: metric = new TMOP_Metric_004; break;
458  case 7: metric = new TMOP_Metric_007; break;
459  case 9: metric = new TMOP_Metric_009; break;
460  case 14: metric = new TMOP_Metric_014; break;
461  case 22: metric = new TMOP_Metric_022(min_detJ); break;
462  case 50: metric = new TMOP_Metric_050; break;
463  case 55: metric = new TMOP_Metric_055; break;
464  case 56: metric = new TMOP_Metric_056; break;
465  case 58: metric = new TMOP_Metric_058; break;
466  case 66: metric = new TMOP_Metric_066(0.5); break;
467  case 77: metric = new TMOP_Metric_077; break;
468  case 80: metric = new TMOP_Metric_080(0.5); break;
469  case 85: metric = new TMOP_Metric_085; break;
470  case 90: metric = new TMOP_Metric_090; break;
471  case 94: metric = new TMOP_Metric_094; break;
472  case 98: metric = new TMOP_Metric_098; break;
473  // case 211: metric = new TMOP_Metric_211; break;
474  // case 252: metric = new TMOP_Metric_252(min_detJ); break;
475  case 301: metric = new TMOP_Metric_301; break;
476  case 302: metric = new TMOP_Metric_302; break;
477  case 303: metric = new TMOP_Metric_303; break;
478  case 304: metric = new TMOP_Metric_304; break;
479  // case 311: metric = new TMOP_Metric_311; break;
480  case 313: metric = new TMOP_Metric_313(min_detJ); break;
481  case 315: metric = new TMOP_Metric_315; break;
482  case 316: metric = new TMOP_Metric_316; break;
483  case 321: metric = new TMOP_Metric_321; break;
484  case 322: metric = new TMOP_Metric_322; break;
485  case 323: metric = new TMOP_Metric_323; break;
486  case 328: metric = new TMOP_Metric_328; break;
487  case 332: metric = new TMOP_Metric_332(0.5); break;
488  case 333: metric = new TMOP_Metric_333(0.5); break;
489  case 334: metric = new TMOP_Metric_334(0.5); break;
490  case 338: metric = new TMOP_Metric_338; break;
491  case 347: metric = new TMOP_Metric_347(0.5); break;
492  // case 352: metric = new TMOP_Metric_352(min_detJ); break;
493  case 360: metric = new TMOP_Metric_360; break;
494  // A-metrics
495  case 11: metric = new TMOP_AMetric_011; break;
496  case 36: metric = new TMOP_AMetric_036; break;
497  case 107: metric = new TMOP_AMetric_107a; break;
498  case 126: metric = new TMOP_AMetric_126(0.9); break;
499  default:
500  if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
501  return 3;
502  }
503  TMOP_QualityMetric *h_metric = NULL;
505  {
506  switch (h_metric_id)
507  {
508  case 1: h_metric = new TMOP_Metric_001; break;
509  case 2: h_metric = new TMOP_Metric_002; break;
510  case 7: h_metric = new TMOP_Metric_007; break;
511  case 9: h_metric = new TMOP_Metric_009; break;
512  case 55: h_metric = new TMOP_Metric_055; break;
513  case 56: h_metric = new TMOP_Metric_056; break;
514  case 58: h_metric = new TMOP_Metric_058; break;
515  case 77: h_metric = new TMOP_Metric_077; break;
516  case 315: h_metric = new TMOP_Metric_315; break;
517  case 316: h_metric = new TMOP_Metric_316; break;
518  case 321: h_metric = new TMOP_Metric_321; break;
519  default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
520  endl;
521  return 3;
522  }
523  }
524
526  switch (barrier_type)
527  {
529  break;
531  break;
533  break;
534  default: cout << "barrier_type not supported: " << barrier_type << endl;
535  return 3;
536  }
537
539  switch (worst_case_type)
540  {
542  break;
544  break;
546  break;
547  default: cout << "worst_case_type not supported: " << worst_case_type << endl;
548  return 3;
549  }
550
551  TMOP_QualityMetric *untangler_metric = NULL;
552  if (barrier_type > 0 || worst_case_type > 0)
553  {
554  if (barrier_type > 0)
555  {
556  MFEM_VERIFY(metric_id == 4 || metric_id == 14 || metric_id == 66,
557  "Metric not supported for shifted/pseudo barriers.");
558  }
559  untangler_metric = new TMOP_WorstCaseUntangleOptimizer_Metric(*metric,
560  2,
561  1.5,
562  0.001,//0.01 for pseudo barrier
563  0.001,
564  btype,
565  wctype);
566  }
567
568  if (metric_id < 300 || h_metric_id < 300)
569  {
570  MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
571  }
572  if (metric_id >= 300 || h_metric_id >= 300)
573  {
574  MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
575  }
576
578  TargetConstructor *target_c = NULL;
581  H1_FECollection ind_fec(mesh_poly_deg, dim);
582  ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
583  ParFiniteElementSpace ind_fesv(pmesh, &ind_fec, dim);
584  ParGridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
585  ParGridFunction aspr3d(&ind_fesv);
586
587  const AssemblyLevel al =
589
590  switch (target_id)
591  {
592  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
593  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
594  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
595  case 4:
596  {
599  adapt_coeff = new HessianCoefficient(dim, metric_id);
601  target_c = tc;
602  break;
603  }
604  case 5: // Discrete size 2D or 3D
605  {
609  {
611  }
612  else
613  {
614 #ifdef MFEM_USE_GSLIB
616 #else
617  MFEM_ABORT("MFEM is not built with GSLIB.");
618 #endif
619  }
620  ConstructSizeGF(size);
621  tc->SetParDiscreteTargetSize(size);
622  target_c = tc;
623  break;
624  }
625  case 6: // material indicator 2D
626  {
627  ParGridFunction d_x(&ind_fes), d_y(&ind_fes), disc(&ind_fes);
628
632  disc.ProjectCoefficient(mat_coeff);
634  {
636  }
637  else
638  {
639 #ifdef MFEM_USE_GSLIB
641 #else
642  MFEM_ABORT("MFEM is not built with GSLIB.");
643 #endif
644  }
645  // Diffuse the interface
646  DiffuseField(disc,2);
647
648  // Get partials with respect to x and y of the grid function
649  disc.GetDerivative(1,0,d_x);
650  disc.GetDerivative(1,1,d_y);
651
652  // Compute the squared magnitude of the gradient
653  for (int i = 0; i < size.Size(); i++)
654  {
655  size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
656  }
657  const double max = size.Max();
658  double max_all;
659  MPI_Allreduce(&max, &max_all, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
660
661  for (int i = 0; i < d_x.Size(); i++)
662  {
663  d_x(i) = std::abs(d_x(i));
664  d_y(i) = std::abs(d_y(i));
665  }
666  const double eps = 0.01;
667  const double aspr_ratio = 20.0;
668  const double size_ratio = 40.0;
669
670  for (int i = 0; i < size.Size(); i++)
671  {
672  size(i) = (size(i)/max_all);
673  aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
674  aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
675  if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
676  if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
677  }
678  Vector vals;
679  const int NE = pmesh->GetNE();
680  double volume = 0.0, volume_ind = 0.0;
681
682  for (int i = 0; i < NE; i++)
683  {
685  const IntegrationRule &ir =
686  IntRules.Get(pmesh->GetElementBaseGeometry(i), Tr->OrderJ());
687  size.GetValues(i, ir, vals);
688  for (int j = 0; j < ir.GetNPoints(); j++)
689  {
690  const IntegrationPoint &ip = ir.IntPoint(j);
691  Tr->SetIntPoint(&ip);
692  volume += ip.weight * Tr->Weight();
693  volume_ind += vals(j) * ip.weight * Tr->Weight();
694  }
695  }
696  double volume_all, volume_ind_all;
697  MPI_Allreduce(&volume, &volume_all, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
698  MPI_Allreduce(&volume_ind, &volume_ind_all, 1, MPI_DOUBLE, MPI_SUM,
699  MPI_COMM_WORLD);
700  const int NE_ALL = pmesh->GetGlobalNE();
701
702  const double avg_zone_size = volume_all / NE_ALL;
703
704  const double small_avg_ratio =
705  (volume_ind_all + (volume_all - volume_ind_all) / size_ratio)
706  / volume_all;
707
708  const double small_zone_size = small_avg_ratio * avg_zone_size;
709  const double big_zone_size = size_ratio * small_zone_size;
710
711  for (int i = 0; i < size.Size(); i++)
712  {
713  const double val = size(i);
714  const double a = (big_zone_size - small_zone_size) / small_zone_size;
715  size(i) = big_zone_size / (1.0+a*val);
716  }
717
718  DiffuseField(size, 2);
719  DiffuseField(aspr, 2);
720
721  tc->SetParDiscreteTargetSize(size);
723  target_c = tc;
724  break;
725  }
726  case 7: // Discrete aspect ratio 3D
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  }
743  aspr3d.ProjectCoefficient(fd_aspr3d);
745  target_c = tc;
746  break;
747  }
748  case 8: // shape/size + orientation 2D
749  {
753  {
755  }
756  else
757  {
758 #ifdef MFEM_USE_GSLIB
760 #else
761  MFEM_ABORT("MFEM is not built with GSLIB.");
762 #endif
763  }
764
765  ConstantCoefficient size_coeff(0.1*0.1);
766  size.ProjectCoefficient(size_coeff);
767  tc->SetParDiscreteTargetSize(size);
768
770  ori.ProjectCoefficient(ori_coeff);
772  target_c = tc;
773  break;
774  }
775  // Targets used for hr-adaptivity tests.
776  case 9: // size target in an annular region.
777  case 10: // size+aspect-ratio in an annular region.
778  case 11: // size+aspect-ratio target for a rotate sine wave
779  {
782  hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
784  target_c = tc;
785  break;
786  }
787  default:
788  if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
789  return 3;
790  }
791  if (target_c == NULL)
792  {
793  target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
794  }
795  target_c->SetNodes(x0);
796
797  // Automatically balanced gamma in composite metrics.
798  auto metric_combo = dynamic_cast<TMOP_Combo_QualityMetric *>(metric);
799  if (metric_combo && bal_expl_combo)
800  {
801  Vector bal_weights;
802  metric_combo->ComputeBalancedWeights(x, *target_c, bal_weights);
803  metric_combo->SetWeights(bal_weights);
804  }
805
806  TMOP_QualityMetric *metric_to_use = barrier_type > 0 || worst_case_type > 0
807  ? untangler_metric
808  : metric;
809  auto tmop_integ = new TMOP_Integrator(metric_to_use, target_c, h_metric);
810  tmop_integ->IntegrateOverTarget(integ_over_targ);
811  if (barrier_type > 0 || worst_case_type > 0)
812  {
813  tmop_integ->ComputeUntangleMetricQuantiles(x, *pfespace);
814  }
815
816  // Finite differences for computations of derivatives.
817  if (fdscheme)
818  {
819  MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
820  tmop_integ->EnableFiniteDifferences(x);
821  }
822  tmop_integ->SetExactActionFlag(exactaction);
823
824  // Setup the quadrature rules for the TMOP integrator.
825  IntegrationRules *irules = NULL;
827  {
828  case 1: irules = &IntRulesLo; break;
829  case 2: irules = &IntRules; break;
830  case 3: irules = &IntRulesCU; break;
831  default:
832  if (myid == 0) { cout << "Unknown quad_type: " << quad_type << endl; }
833  return 3;
834  }
836  if (myid == 0 && dim == 2)
837  {
838  cout << "Triangle quadrature points: "
841  << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
842  }
843  if (myid == 0 && dim == 3)
844  {
845  cout << "Tetrahedron quadrature points: "
847  << "\nHexahedron quadrature points: "
849  << "\nPrism quadrature points: "
850  << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
851  }
852
853  // Limit the node movement.
854  // The limiting distances can be given by a general function of space.
855  ParFiniteElementSpace dist_pfespace(pmesh, fec); // scalar space
856  ParGridFunction dist(&dist_pfespace);
857  dist = 1.0;
858  // The small_phys_size is relevant only with proper normalization.
859  if (normalization) { dist = small_phys_size; }
860  ConstantCoefficient lim_coeff(lim_const);
861  if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
862
868  {
869  MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
870
873
875  else if (adapt_eval == 1)
876  {
877 #ifdef MFEM_USE_GSLIB
879 #else
880  MFEM_ABORT("MFEM is not built with GSLIB support!");
881 #endif
882  }
883  else { MFEM_ABORT("Bad interpolation option."); }
884
887  if (visualization)
888  {
889  socketstream vis1;
890  common::VisualizeField(vis1, "localhost", 19916, adapt_lim_gf0, "Zeta 0",
891  300, 600, 300, 300);
892  }
893  }
894
895  // Has to be after the enabling of the limiting / alignment, as it computes
896  // normalization factors for these terms as well.
897  if (normalization) { tmop_integ->ParEnableNormalization(x0); }
898
899  // 13. Setup the final NonlinearForm (which defines the integral of interest,
900  // its first and second derivatives). Here we can use a combination of
901  // metrics, i.e., optimize the sum of two integrals, where both are
902  // scaled by used-defined space-dependent weights. Note that there are
903  // no command-line options for the weights and the type of the second
904  // metric; one should update those in the code.
905  ParNonlinearForm a(pfespace);
906  if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
907  ConstantCoefficient *metric_coeff1 = NULL;
908  TMOP_QualityMetric *metric2 = NULL;
909  TargetConstructor *target_c2 = NULL;
910  FunctionCoefficient metric_coeff2(weight_fun);
911
912  // Explicit combination of metrics.
913  if (combomet > 0)
914  {
915  // First metric.
916  metric_coeff1 = new ConstantCoefficient(1.0);
917  tmop_integ->SetCoefficient(*metric_coeff1);
918
919  // Second metric.
920  if (dim == 2) { metric2 = new TMOP_Metric_077; }
921  else { metric2 = new TMOP_Metric_315; }
922  TMOP_Integrator *tmop_integ2 = NULL;
923  if (combomet == 1)
924  {
925  target_c2 = new TargetConstructor(
927  target_c2->SetVolumeScale(0.01);
928  target_c2->SetNodes(x0);
929  tmop_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
930  tmop_integ2->SetCoefficient(metric_coeff2);
931  }
932  else { tmop_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
933  tmop_integ2->IntegrateOverTarget(integ_over_targ);
935  if (fdscheme) { tmop_integ2->EnableFiniteDifferences(x); }
936  tmop_integ2->SetExactActionFlag(exactaction);
937
941  if (normalization) { combo->ParEnableNormalization(x0); }
942  if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
943
945  }
946  else
947  {
949  }
950
951  if (pa) { a.Setup(); }
952
953  // Compute the minimum det(J) of the starting mesh.
954  min_detJ = infinity();
955  const int NE = pmesh->GetNE();
956  for (int i = 0; i < NE; i++)
957  {
958  const IntegrationRule &ir =
961  for (int j = 0; j < ir.GetNPoints(); j++)
962  {
963  transf->SetIntPoint(&ir.IntPoint(j));
964  min_detJ = min(min_detJ, transf->Jacobian().Det());
965  }
966  }
967  double minJ0;
968  MPI_Allreduce(&min_detJ, &minJ0, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
969  min_detJ = minJ0;
970  if (myid == 0)
971  { cout << "Minimum det(J) of the original mesh is " << min_detJ << endl; }
972
973  if (min_detJ < 0.0 && barrier_type == 0
974  && metric_id != 22 && metric_id != 211 && metric_id != 252
975  && metric_id != 311 && metric_id != 313 && metric_id != 352)
976  {
977  MFEM_ABORT("The input mesh is inverted! Try an untangling metric.");
978  }
979  if (min_detJ < 0.0)
980  {
981  MFEM_VERIFY(target_t == TargetConstructor::IDEAL_SHAPE_UNIT_SIZE,
982  "Untangling is supported only for ideal targets.");
983
984  const DenseMatrix &Wideal =
986  min_detJ /= Wideal.Det();
987
988  double h0min = h0.Min(), h0min_all;
989  MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
990  // Slightly below minJ0 to avoid div by 0.
991  min_detJ -= 0.01 * h0min_all;
992  }
993
994  // For HR tests, the energy is normalized by the number of elements.
995  const double init_energy = a.GetParGridFunctionEnergy(x) /
996  (hradaptivity ? pmesh->GetGlobalNE() : 1);
997  double init_metric_energy = init_energy;
998  if (lim_const > 0.0 || adapt_lim_const > 0.0)
999  {
1000  lim_coeff.constant = 0.0;
1002  init_metric_energy = a.GetParGridFunctionEnergy(x) /
1003  (hradaptivity ? pmesh->GetGlobalNE() : 1);
1004  lim_coeff.constant = lim_const;
1006  }
1007
1008  // Visualize the starting mesh and metric values.
1009  // Note that for combinations of metrics, this only shows the first metric.
1010  if (visualization)
1011  {
1012  char title[] = "Initial metric values";
1013  vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
1014  }
1015
1016  // 14. Fix all boundary nodes, or fix only a given component depending on the
1017  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
1018  // fixed x/y/z components of the node. Attribute dim+1 corresponds to
1019  // an entirely fixed node.
1020  if (move_bnd == false)
1021  {
1022  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
1023  ess_bdr = 1;
1024  a.SetEssentialBC(ess_bdr);
1025  }
1026  else
1027  {
1028  int n = 0;
1029  for (int i = 0; i < pmesh->GetNBE(); i++)
1030  {
1031  const int nd = pfespace->GetBE(i)->GetDof();
1032  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
1033  MFEM_VERIFY(!(dim == 2 && attr == 3),
1034  "Boundary attribute 3 must be used only for 3D meshes. "
1035  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
1037  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
1038  if (attr == 4) { n += nd * dim; }
1039  }
1040  Array<int> ess_vdofs(n);
1041  n = 0;
1042  for (int i = 0; i < pmesh->GetNBE(); i++)
1043  {
1044  const int nd = pfespace->GetBE(i)->GetDof();
1045  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
1046  pfespace->GetBdrElementVDofs(i, vdofs);
1047  if (attr == 1) // Fix x components.
1048  {
1049  for (int j = 0; j < nd; j++)
1050  { ess_vdofs[n++] = vdofs[j]; }
1051  }
1052  else if (attr == 2) // Fix y components.
1053  {
1054  for (int j = 0; j < nd; j++)
1055  { ess_vdofs[n++] = vdofs[j+nd]; }
1056  }
1057  else if (attr == 3) // Fix z components.
1058  {
1059  for (int j = 0; j < nd; j++)
1060  { ess_vdofs[n++] = vdofs[j+2*nd]; }
1061  }
1062  else if (attr == 4) // Fix all components.
1063  {
1064  for (int j = 0; j < vdofs.Size(); j++)
1065  { ess_vdofs[n++] = vdofs[j]; }
1066  }
1067  }
1068  a.SetEssentialVDofs(ess_vdofs);
1069  }
1070
1071  // As we use the inexact Newton method to solve the resulting nonlinear
1072  // system, here we setup the linear solver for the system's Jacobian.
1073  Solver *S = NULL, *S_prec = NULL;
1074  const double linsol_rtol = 1e-12;
1075  // Level of output.
1076  IterativeSolver::PrintLevel linsolver_print;
1077  if (verbosity_level == 2)
1078  { linsolver_print.Errors().Warnings().FirstAndLast(); }
1079  if (verbosity_level > 2)
1080  { linsolver_print.Errors().Warnings().Iterations(); }
1081  if (lin_solver == 0)
1082  {
1083  S = new DSmoother(1, 1.0, max_lin_iter);
1084  }
1085  else if (lin_solver == 1)
1086  {
1087  CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
1088  cg->SetMaxIter(max_lin_iter);
1089  cg->SetRelTol(linsol_rtol);
1090  cg->SetAbsTol(0.0);
1091  cg->SetPrintLevel(linsolver_print);
1092  S = cg;
1093  }
1094  else
1095  {
1096  MINRESSolver *minres = new MINRESSolver(MPI_COMM_WORLD);
1097  minres->SetMaxIter(max_lin_iter);
1098  minres->SetRelTol(linsol_rtol);
1099  minres->SetAbsTol(0.0);
1100  minres->SetPrintLevel(linsolver_print);
1101  if (lin_solver == 3 || lin_solver == 4)
1102  {
1103  if (pa)
1104  {
1105  MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
1106  auto js = new OperatorJacobiSmoother;
1107  js->SetPositiveDiagonal(true);
1108  S_prec = js;
1109  }
1110  else
1111  {
1112  auto hs = new HypreSmoother;
1113  hs->SetType((lin_solver == 3) ? HypreSmoother::Jacobi
1114  /* */ : HypreSmoother::l1Jacobi, 1);
1115  hs->SetPositiveDiagonal(true);
1116  S_prec = hs;
1117  }
1118  minres->SetPreconditioner(*S_prec);
1119  }
1120  S = minres;
1121  }
1122
1123  //
1124  // Perform the nonlinear optimization.
1125  //
1126  const IntegrationRule &ir =
1128  TMOPNewtonSolver solver(pfespace->GetComm(), ir, solver_type);
1129  // Provide all integration rules in case of a mixed mesh.
1131  // Specify linear solver when we use a Newton-based solver.
1132  if (solver_type == 0) { solver.SetPreconditioner(*S); }
1133  // For untangling, the solver will update the min det(T) values.
1134  solver.SetMinDetPtr(&min_detJ);
1135  solver.SetMaxIter(solver_iter);
1136  solver.SetRelTol(solver_rtol);
1137  solver.SetAbsTol(0.0);
1138  if (solver_art_type > 0)
1139  {
1141  }
1142  // Level of output.
1143  IterativeSolver::PrintLevel newton_print;
1144  if (verbosity_level > 0)
1145  { newton_print.Errors().Warnings().Iterations(); }
1146  solver.SetPrintLevel(newton_print);
1148  // If hr-adaptivity is disabled, r-adaptivity is done once using the
1149  // TMOPNewtonSolver.
1150  // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
1152  // The solver terminates if an h-adaptivity iteration does not modify
1153  // any element in the mesh.
1154  TMOPHRSolver hr_solver(*pmesh, a, solver,
1156  mesh_poly_deg, h_metric_id,
1157  n_hr_iter, n_h_iter);
1160  {
1163  }
1164  hr_solver.Mult();
1165
1166  // 16. Save the optimized mesh to a file. This output can be viewed later
1167  // using GLVis: "glvis -m optimized -np num_mpi_tasks".
1168  {
1169  ostringstream mesh_name;
1170  mesh_name << "optimized.mesh";
1171  ofstream mesh_ofs(mesh_name.str().c_str());
1172  mesh_ofs.precision(8);
1173  pmesh->PrintAsOne(mesh_ofs);
1174  }
1175
1176  // Report the final energy of the functional.
1177  const double fin_energy = a.GetParGridFunctionEnergy(x) /
1178  (hradaptivity ? pmesh->GetGlobalNE() : 1);
1179  double fin_metric_energy = fin_energy;
1180  if (lim_const > 0.0 || adapt_lim_const > 0.0)
1181  {
1182  lim_coeff.constant = 0.0;
1184  fin_metric_energy = a.GetParGridFunctionEnergy(x) /
1185  (hradaptivity ? pmesh->GetGlobalNE() : 1);
1186  lim_coeff.constant = lim_const;
1188  }
1189  if (myid == 0)
1190  {
1191  std::cout << std::scientific << std::setprecision(4);
1192  cout << "Initial strain energy: " << init_energy
1193  << " = metrics: " << init_metric_energy
1194  << " + extra terms: " << init_energy - init_metric_energy << endl;
1195  cout << " Final strain energy: " << fin_energy
1196  << " = metrics: " << fin_metric_energy
1197  << " + extra terms: " << fin_energy - fin_metric_energy << endl;
1198  cout << "The strain energy decreased by: "
1199  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
1200  }
1201
1202  // Visualize the final mesh and metric values.
1203  if (visualization)
1204  {
1205  char title[] = "Final metric values";
1206  vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
1207  }
1208
1209  if (adapt_lim_const > 0.0 && visualization)
1210  {
1211  socketstream vis0;
1212  common::VisualizeField(vis0, "localhost", 19916, adapt_lim_gf0, "Xi 0",
1213  600, 600, 300, 300);
1214  }
1215
1216  // Visualize the mesh displacement.
1217  if (visualization)
1218  {
1219  x0 -= x;
1220  socketstream sock;
1221  if (myid == 0)
1222  {
1223  sock.open("localhost", 19916);
1224  sock << "solution\n";
1225  }
1226  pmesh->PrintAsOne(sock);
1227  x0.SaveAsOne(sock);
1228  if (myid == 0)
1229  {
1230  sock << "window_title 'Displacements'\n"
1231  << "window_geometry "
1232  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
1233  << "keys jRmclA" << endl;
1234  }
1235  }
1236
1237  delete S;
1238  delete S_prec;
1239  delete target_c2;
1240  delete metric2;
1241  delete metric_coeff1;
1243  delete target_c;
1246  delete h_metric;
1247  delete metric;
1248  delete untangler_metric;
1249  delete pfespace;
1250  delete fec;
1251  delete pmesh;
1252
1253  return 0;
1254 }
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
Definition: tmop_tools.cpp:879
void discrete_aspr_3d(const Vector &x, Vector &v)
virtual void SetParDiscreteTargetSize(const ParGridFunction &tspec_)
Definition: tmop.cpp:1847
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:253
Definition: ex25.cpp:149
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:96
3D non-barrier Shape (S) metric.
Definition: tmop.hpp:1051
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:4497
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:843
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:980
void SetPositiveDiagonal(bool pos_diag=true)
Replace diagonal entries with their absolute values.
Definition: solvers.hpp:339
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:207
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:967
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:331
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:935
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:703
Definition: tmop.hpp:1658
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
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
void ConstructSizeGF(GridFunction &size)
double Det() const
Definition: densemat.cpp:487
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:412
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:682
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:1415
Parallel non-linear operator on the true dofs.
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:706
3D Size (V) metric.
Definition: tmop.hpp:785
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop.hpp:2004
3D compound barrier Shape+Size (VS) metric (polyconvex, balanced).
Definition: tmop.hpp:989
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:1145
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:885
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
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:561
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
Definition: tmop.hpp:2224
STL namespace.
MINRES method.
Definition: solvers.hpp:603
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:817
void IntegrateOverTarget(bool integ_over_target_)
Definition: tmop.hpp:2013
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:340
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:864
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
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:4685
2D Shifted barrier form of shape metric (mu_2).
Definition: tmop.hpp:389
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:466
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:724
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
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:256
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
2D compound barrier Shape+Size (VS) metric (balanced).
Definition: tmop.hpp:588
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:1115
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:83
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:2028
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9888
3D Size (V) metric.
Definition: tmop.hpp:803
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:10232
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:1409
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:1010
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:764
int main(int argc, char *argv[])
virtual void SetParDiscreteTargetAspectRatio(const ParGridFunction &tspec_)
Definition: tmop.cpp:1867
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:4800
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:609
Definition: tmop_amr.hpp:253
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:1127
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:470
A general vector function coefficient.
3D compound barrier Shape+Size (VS) metric (polyconvex, balanced).
Definition: tmop.hpp:906
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:534
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:4969
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
Adds a new TMOP_Integrator to the combination.
Definition: tmop.hpp:2254
double material_indicator_2d(const Vector &x)
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:358
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
double a
Definition: lissajous.cpp:41
Class for integration point with weight.
Definition: intrules.hpp:31
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:355
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:22
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i&#39;th boundary element. The returned indices are offsets int...
Definition: fespace.cpp:297
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:1109
3D Shape (S) metric, untangling version of 303.
Definition: tmop.hpp:764
int dim
Definition: ex24.cpp:53
2D compound barrier Shape+Size (VS) metric (balanced).
Definition: tmop.hpp:567
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:745
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:44
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:269
A general function coefficient.
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition: fespace.cpp:251
Vector data type.
Definition: vector.hpp:58
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:948
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition: pmesh.cpp:2076
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:1333
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
virtual void SetParDiscreteTargetOrientation(const ParGridFunction &tspec_)
Definition: tmop.cpp:1877
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:5624
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1329
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:1701
3D compound barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:927
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:3166
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:3429
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:552
Definition: vector.hpp:473
double GetElementVolume(int i)
Definition: mesh.cpp:119
Settings for the output behavior of the IterativeSolver.
Definition: solvers.hpp:78
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:663
2D non-barrier metric without a type.
Definition: tmop.hpp:221
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1733
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:374