MFEM  v4.4.0 Finite element discretization library
mesh-optimizer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // --------------------------------------------------
13 // Mesh Optimizer Miniapp: Optimize high-order meshes
14 // --------------------------------------------------
15 //
16 // This miniapp performs mesh optimization using the Target-Matrix Optimization
17 // Paradigm (TMOP) by P.Knupp et al., and a global variational minimization
18 // approach. It minimizes the quantity sum_T int_T mu(J(x)), where T are the
19 // target (ideal) elements, J is the Jacobian of the transformation from the
20 // target to the physical element, and mu is the mesh quality metric. This
21 // metric can measure shape, size or alignment of the region around each
22 // quadrature point. The combination of targets & quality metrics is used to
23 // optimize the physical node positions, i.e., they must be as close as possible
24 // to the shape / size / alignment of their targets. This code also demonstrates
25 // a possible use of nonlinear operators (the class TMOP_QualityMetric, defining
26 // mu(J), and the class TMOP_Integrator, defining int mu(J)), as well as their
27 // coupling to Newton methods for solving minimization problems. Note that the
28 // utilized Newton methods are oriented towards avoiding invalid meshes with
29 // negative Jacobian determinants. Each Newton step requires the inversion of a
30 // Jacobian matrix, which is done through an inner linear solver.
31 //
32 // Compile with: make mesh-optimizer
33 //
34 // Sample runs:
36 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -bnd -qt 1 -qo 8
38 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
40 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
41 //
43 // mesh-optimizer -m square01.mesh -o 2 -tid 9 -ni 50 -li 20 -hmid 55 -mid 7 -hr
44 // mesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
45 // mesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
46 //
48 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
49 // Adapted discrete size 3D with PA:
50 // mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa
51 // Adapted discrete size 3D with PA on device (requires CUDA).
52 // * mesh-optimizer -m cube.mesh -o 3 -rs 3 -mid 321 -tid 5 -ls 3 -nor -lc 0.1 -pa -d cuda
53 // Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
54 // mesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
56 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
57 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
58 // Adapted discrete size+orientation (requires GSLIB):
59 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -ae 1 -nor
60 // Adapted discrete aspect-ratio+orientation (requires GSLIB):
61 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 8 -ni 10 -bnd -qt 1 -qo 8 -fd -ae 1
62 // Adapted discrete aspect ratio (3D):
63 // mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -bnd -qt 1 -qo 8
64 //
66 // mesh-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 // mesh-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 // * mesh-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 // mesh-optimizer -m square01.mesh -o 3 -rs 1 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 5e4 -rtol 1e-5 -nor
74 // mesh-optimizer -m square01-tri.mesh -o 3 -rs 0 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 1e4 -rtol 1e-5 -nor
75 // Surface fitting with weight adaptation and termination based on fitting error
76 // mesh-optimizer -m square01.mesh -o 2 -rs 1 -mid 2 -tid 1 -ni 100 -vl 2 -sfc 10 -rtol 1e-20 -st 0 -sfa -sft 1e-5
77 //
79 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
80 // Blade shape with FD-based solver:
81 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
83 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
84 // ICF shape and equal size:
85 // mesh-optimizer -o 3 -mid 9 -tid 2 -ni 25 -ls 3 -art 2 -qo 5
86 // ICF shape and initial size:
87 // mesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
88 // ICF shape:
89 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
90 // ICF limited shape:
91 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8 -lc 10
92 // ICF combo shape + size (rings, slow convergence):
93 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 1000 -bnd -qt 1 -qo 8 -cmb 1
94 // Mixed tet / cube / hex mesh with limiting:
95 // mesh-optimizer -m ../../data/fichera-mixed-p2.mesh -o 4 -rs 1 -mid 301 -tid 1 -fix-bnd -qo 6 -nor -lc 0.25
96 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
97 // * mesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -mid 303 -tid 1 -ni 20 -li 500 -fix-bnd
98 // 2D non-conforming shape and equal size:
99 // mesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -bnd -qt 1 -qo 8
100 //
101 // 2D untangling:
102 // mesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
103 // 3D untangling (the mesh is in the mfem/data GitHub repository):
104 // * mesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
105
106 #include "mfem.hpp"
107 #include "../common/mfem-common.hpp"
108 #include <fstream>
109 #include <iostream>
110 #include "mesh-optimizer.hpp"
111
112 using namespace mfem;
113 using namespace std;
114
115 int main(int argc, char *argv[])
116 {
117  // 0. Set the method's default parameters.
118  const char *mesh_file = "icf.mesh";
119  int mesh_poly_deg = 1;
120  int rs_levels = 0;
121  double jitter = 0.0;
122  int metric_id = 1;
123  int target_id = 1;
124  double lim_const = 0.0;
126  double surface_fit_const = 0.0;
129  int solver_type = 0;
130  int solver_iter = 20;
131  double solver_rtol = 1e-10;
132  int solver_art_type = 0;
133  int lin_solver = 2;
134  int max_lin_iter = 100;
135  bool move_bnd = true;
136  int combomet = 0;
138  int h_metric_id = -1;
139  bool normalization = false;
140  bool visualization = true;
141  int verbosity_level = 0;
142  bool fdscheme = false;
144  bool exactaction = false;
145  const char *devopt = "cpu";
146  bool pa = false;
147  int n_hr_iter = 5;
148  int n_h_iter = 1;
150  double surface_fit_threshold = -10;
151
152  // 1. Parse command-line options.
153  OptionsParser args(argc, argv);
155  "Mesh file to use.");
157  "Polynomial degree of mesh finite element space.");
159  "Number of times to refine the mesh uniformly in serial.");
161  "Random perturbation scaling factor.");
163  "Mesh optimization metric:\n\t"
164  "T-metrics\n\t"
165  "1 : |T|^2 -- 2D shape\n\t"
166  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
167  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
168  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
169  "14 : |T-I|^2 -- 2D shape+size+orientation\n\t"
170  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
171  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
172  "55 : (tau-1)^2 -- 2D size\n\t"
173  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
174  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
175  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
176  "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t"
177  "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t"
178  "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t"
179  // "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
180  // "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
181  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
182  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
183  "303: (|T|^2)/3*tau^(2/3)-1 -- 3D shape\n\t"
184  //"311: (tau-1)^2-tau+sqrt(tau^2+eps)-- 3D untangling\n\t"
185  "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t"
186  "315: (tau-1)^2 -- 3D size\n\t"
187  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D size\n\t"
188  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
189  // "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling\n\t"
190  "A-metrics\n\t"
191  "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t"
192  "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t"
193  "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t"
194  "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t"
195  );
197  "Target (ideal element) type:\n\t"
198  "1: Ideal shape, unit size\n\t"
199  "2: Ideal shape, equal size\n\t"
200  "3: Ideal shape, initial size\n\t"
201  "4: Given full analytic Jacobian (in physical space)\n\t"
202  "5: Ideal shape, given size (in physical space)");
203  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
207  "Surface preservation constant.");
210  "1: Gauss-Lobatto\n\t"
211  "2: Gauss-Legendre\n\t"
212  "3: Closed uniform points");
214  "Order of the quadrature rule.");
216  " Type of solver: (default) 0: Newton, 1: LBFGS");
218  "Maximum number of Newton iterations.");
220  "Relative tolerance for the Newton solver.");
222  "Type of adaptive relative linear solver tolerance:\n\t"
223  "0: None (default)\n\t"
224  "1: Eisenstat-Walker type 1\n\t"
225  "2: Eisenstat-Walker type 2");
227  "Linear solver:\n\t"
228  "0: l1-Jacobi\n\t"
229  "1: CG\n\t"
230  "2: MINRES\n\t"
231  "3: MINRES + Jacobi preconditioner\n\t"
232  "4: MINRES + l1-Jacobi preconditioner");
234  "Maximum number of iterations in the linear solve.");
236  "--fix-boundary",
237  "Enable motion along horizontal and vertical boundaries.");
239  "Combination of metrics options:\n\t"
240  "0: Use single metric\n\t"
241  "1: Shape + space-dependent size given analytically\n\t"
242  "2: Shape + adapted size given discretely; shared target");
247  "Same options as metric_id. Used to determine refinement"
248  " type for each element if h-adaptivity is enabled.");
250  "--no-normalization",
251  "Make all terms in the optimization functional unitless.");
253  "-no-fd", "--no-fd-approx",
254  "Enable finite difference based derivative computations.");
256  "-no-ex", "--no-exact-action",
257  "Enable exact action of TMOP_Integrator.");
259  "--no-visualization",
260  "Enable or disable GLVis visualization.");
262  "Set the verbosity level - 0, 1, or 2.");
264  "0 - Advection based (DEFAULT), 1 - GSLIB.");
266  "Device configuration string, see Device::Configure().");
268  "--no-partial-assembly", "Enable Partial Assembly.");
273  "iteration.");
276  "Enable or disable adaptive surface fitting.");
278  "Set threshold for surface fitting. TMOP solver will"
279  "terminate when max surface fitting error is below this limit");
280  args.Parse();
281  if (!args.Good())
282  {
283  args.PrintUsage(cout);
284  return 1;
285  }
286  args.PrintOptions(cout);
287
288  if (h_metric_id < 0) { h_metric_id = metric_id; }
289
291  {
292  MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
293  " supported on cpus.");
294  }
295  Device device(devopt);
296  device.Print();
297
298  // 2. Initialize and refine the starting mesh.
299  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
300  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
301  const int dim = mesh->Dimension();
302
303  if (hradaptivity) { mesh->EnsureNCMesh(); }
304
305  // 3. Define a finite element space on the mesh-> Here we use vector finite
306  // elements which are tensor products of quadratic finite elements. The
307  // number of components in the vector finite element space is specified by
308  // the last parameter of the FiniteElementSpace constructor.
310  if (mesh_poly_deg <= 0)
311  {
313  mesh_poly_deg = 2;
314  }
315  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
316  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim);
317
318  // 4. Make the mesh curved based on the above finite element space. This
319  // means that we define the mesh elements through a fespace-based
320  // transformation of the reference element.
321  mesh->SetNodalFESpace(fespace);
322
323  // 5. Set up an empty right-hand side vector b, which is equivalent to b=0.
324  Vector b(0);
325
326  // 6. Get the mesh nodes (vertices and other degrees of freedom in the finite
327  // element space) as a finite element grid function in fespace. Note that
328  // changing x automatically changes the shapes of the mesh elements.
329  GridFunction x(fespace);
330  mesh->SetNodalGridFunction(&x);
331
332  // 7. Define a vector representing the minimal local mesh size in the mesh
333  // nodes. We index the nodes using the scalar version of the degrees of
334  // freedom in fespace. Note: this is partition-dependent.
335  //
336  // In addition, compute average mesh size and total volume.
337  Vector h0(fespace->GetNDofs());
338  h0 = infinity();
339  double mesh_volume = 0.0;
340  Array<int> dofs;
341  for (int i = 0; i < mesh->GetNE(); i++)
342  {
343  // Get the local scalar element degrees of freedom in dofs.
344  fespace->GetElementDofs(i, dofs);
345  // Adjust the value of h0 in dofs based on the local mesh size.
346  const double hi = mesh->GetElementSize(i);
347  for (int j = 0; j < dofs.Size(); j++)
348  {
349  h0(dofs[j]) = min(h0(dofs[j]), hi);
350  }
351  mesh_volume += mesh->GetElementVolume(i);
352  }
353  const double small_phys_size = pow(mesh_volume, 1.0 / dim) / 100.0;
354
355  // 8. Add a random perturbation to the nodes in the interior of the domain.
356  // We define a random grid function of fespace and make sure that it is
357  // zero on the boundary and its values are locally of the order of h0.
358  // The latter is based on the DofToVDof() method which maps the scalar to
359  // the vector degrees of freedom in fespace.
360  GridFunction rdm(fespace);
361  rdm.Randomize();
362  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
363  rdm *= jitter;
365  // Scale the random values to be of order of the local mesh size.
366  for (int i = 0; i < fespace->GetNDofs(); i++)
367  {
368  for (int d = 0; d < dim; d++)
369  {
370  rdm(fespace->DofToVDof(i,d)) *= h0(i);
371  }
372  }
373  Array<int> vdofs;
374  for (int i = 0; i < fespace->GetNBE(); i++)
375  {
376  // Get the vector degrees of freedom in the boundary element.
377  fespace->GetBdrElementVDofs(i, vdofs);
378  // Set the boundary values to zero.
379  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
380  }
381  x -= rdm;
382  x.SetTrueVector();
383  x.SetFromTrueVector();
384
385  // 9. Save the starting (prior to the optimization) mesh to a file. This
386  // output can be viewed later using GLVis: "glvis -m perturbed.mesh".
387  {
388  ofstream mesh_ofs("perturbed.mesh");
389  mesh->Print(mesh_ofs);
390  }
391
392  // 10. Store the starting (prior to the optimization) positions.
393  GridFunction x0(fespace);
394  x0 = x;
395
396  // 11. Form the integrator that uses the chosen metric and target.
397  double tauval = -0.1;
398  TMOP_QualityMetric *metric = NULL;
399  switch (metric_id)
400  {
401  // T-metrics
402  case 1: metric = new TMOP_Metric_001; break;
403  case 2: metric = new TMOP_Metric_002; break;
404  case 7: metric = new TMOP_Metric_007; break;
405  case 9: metric = new TMOP_Metric_009; break;
406  case 14: metric = new TMOP_Metric_014; break;
407  case 22: metric = new TMOP_Metric_022(tauval); break;
408  case 50: metric = new TMOP_Metric_050; break;
409  case 55: metric = new TMOP_Metric_055; break;
410  case 56: metric = new TMOP_Metric_056; break;
411  case 58: metric = new TMOP_Metric_058; break;
412  case 77: metric = new TMOP_Metric_077; break;
413  case 80: metric = new TMOP_Metric_080(0.5); break;
414  case 85: metric = new TMOP_Metric_085; break;
415  case 98: metric = new TMOP_Metric_098; break;
416  // case 211: metric = new TMOP_Metric_211; break;
417  // case 252: metric = new TMOP_Metric_252(tauval); break;
418  case 301: metric = new TMOP_Metric_301; break;
419  case 302: metric = new TMOP_Metric_302; break;
420  case 303: metric = new TMOP_Metric_303; break;
421  // case 311: metric = new TMOP_Metric_311; break;
422  case 313: metric = new TMOP_Metric_313(tauval); break;
423  case 315: metric = new TMOP_Metric_315; break;
424  case 316: metric = new TMOP_Metric_316; break;
425  case 321: metric = new TMOP_Metric_321; break;
426  case 328: metric = new TMOP_Metric_328(0.5); break;
427  case 332: metric = new TMOP_Metric_332(0.5); break;
428  case 333: metric = new TMOP_Metric_333(0.5); break;
429  case 334: metric = new TMOP_Metric_334(0.5); break;
430  // case 352: metric = new TMOP_Metric_352(tauval); break;
431  // A-metrics
432  case 11: metric = new TMOP_AMetric_011; break;
433  case 36: metric = new TMOP_AMetric_036; break;
434  case 107: metric = new TMOP_AMetric_107a; break;
435  case 126: metric = new TMOP_AMetric_126(0.9); break;
436  default:
437  cout << "Unknown metric_id: " << metric_id << endl;
438  return 3;
439  }
440  TMOP_QualityMetric *h_metric = NULL;
442  {
443  switch (h_metric_id)
444  {
445  case 1: h_metric = new TMOP_Metric_001; break;
446  case 2: h_metric = new TMOP_Metric_002; break;
447  case 7: h_metric = new TMOP_Metric_007; break;
448  case 9: h_metric = new TMOP_Metric_009; break;
449  case 55: h_metric = new TMOP_Metric_055; break;
450  case 56: h_metric = new TMOP_Metric_056; break;
451  case 58: h_metric = new TMOP_Metric_058; break;
452  case 77: h_metric = new TMOP_Metric_077; break;
453  case 315: h_metric = new TMOP_Metric_315; break;
454  case 316: h_metric = new TMOP_Metric_316; break;
455  case 321: h_metric = new TMOP_Metric_321; break;
456  default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
457  endl;
458  return 3;
459  }
460  }
461
462  if (metric_id < 300 || h_metric_id < 300)
463  {
464  MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
465  }
466  if (metric_id >= 300 || h_metric_id >= 300)
467  {
468  MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
469  }
470
472  TargetConstructor *target_c = NULL;
475  H1_FECollection ind_fec(mesh_poly_deg, dim);
476  FiniteElementSpace ind_fes(mesh, &ind_fec);
477  FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
478  GridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
479  GridFunction aspr3d(&ind_fesv);
480
481  const AssemblyLevel al =
483
484  switch (target_id)
485  {
486  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
487  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
488  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
489  case 4: // Analytic
490  {
493  adapt_coeff = new HessianCoefficient(dim, metric_id);
495  target_c = tc;
496  break;
497  }
498  case 5: // Discrete size 2D or 3D
499  {
503  {
505  }
506  else
507  {
508 #ifdef MFEM_USE_GSLIB
510 #else
511  MFEM_ABORT("MFEM is not built with GSLIB.");
512 #endif
513  }
514  if (dim == 2)
515  {
517  size.ProjectCoefficient(size_coeff);
518  }
519  else if (dim == 3)
520  {
522  size.ProjectCoefficient(size_coeff);
523  }
524  tc->SetSerialDiscreteTargetSize(size);
525  target_c = tc;
526  break;
527  }
528  case 6: // Discrete size + aspect ratio - 2D
529  {
530  GridFunction d_x(&ind_fes), d_y(&ind_fes), disc(&ind_fes);
531
535  disc.ProjectCoefficient(mat_coeff);
537  {
539  }
540  else
541  {
542 #ifdef MFEM_USE_GSLIB
544 #else
545  MFEM_ABORT("MFEM is not built with GSLIB.");
546 #endif
547  }
548
549  // Diffuse the interface
550  DiffuseField(disc,2);
551
552  // Get partials with respect to x and y of the grid function
553  disc.GetDerivative(1,0,d_x);
554  disc.GetDerivative(1,1,d_y);
555
556  // Compute the squared magnitude of the gradient
557  for (int i = 0; i < size.Size(); i++)
558  {
559  size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
560  }
561  const double max = size.Max();
562
563  for (int i = 0; i < d_x.Size(); i++)
564  {
565  d_x(i) = std::abs(d_x(i));
566  d_y(i) = std::abs(d_y(i));
567  }
568  const double eps = 0.01;
569  const double aspr_ratio = 20.0;
570  const double size_ratio = 40.0;
571
572  for (int i = 0; i < size.Size(); i++)
573  {
574  size(i) = (size(i)/max);
575  aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
576  aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
577  if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
578  if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
579  }
580  Vector vals;
581  const int NE = mesh->GetNE();
582  double volume = 0.0, volume_ind = 0.0;
583
584  for (int i = 0; i < NE; i++)
585  {
587  const IntegrationRule &ir =
588  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr->OrderJ());
589  size.GetValues(i, ir, vals);
590  for (int j = 0; j < ir.GetNPoints(); j++)
591  {
592  const IntegrationPoint &ip = ir.IntPoint(j);
593  Tr->SetIntPoint(&ip);
594  volume += ip.weight * Tr->Weight();
595  volume_ind += vals(j) * ip.weight * Tr->Weight();
596  }
597  }
598
599  const double avg_zone_size = volume / NE;
600
601  const double small_avg_ratio = (volume_ind + (volume - volume_ind) /
602  size_ratio) /
603  volume;
604
605  const double small_zone_size = small_avg_ratio * avg_zone_size;
606  const double big_zone_size = size_ratio * small_zone_size;
607
608  for (int i = 0; i < size.Size(); i++)
609  {
610  const double val = size(i);
611  const double a = (big_zone_size - small_zone_size) / small_zone_size;
612  size(i) = big_zone_size / (1.0+a*val);
613  }
614
615  DiffuseField(size, 2);
616  DiffuseField(aspr, 2);
617
618  tc->SetSerialDiscreteTargetSize(size);
620  target_c = tc;
621  break;
622  }
623  case 7: // Discrete aspect ratio 3D
624  {
628  {
630  }
631  else
632  {
633 #ifdef MFEM_USE_GSLIB
635 #else
636  MFEM_ABORT("MFEM is not built with GSLIB.");
637 #endif
638  }
640  aspr3d.ProjectCoefficient(fd_aspr3d);
641
643  target_c = tc;
644  break;
645  }
646  case 8: // shape/size + orientation 2D
647  {
651  {
653  }
654  else
655  {
656 #ifdef MFEM_USE_GSLIB
658 #else
659  MFEM_ABORT("MFEM is not built with GSLIB.");
660 #endif
661  }
662
663  if (metric_id == 14 || metric_id == 36)
664  {
665  ConstantCoefficient size_coeff(0.1*0.1);
666  size.ProjectCoefficient(size_coeff);
667  tc->SetSerialDiscreteTargetSize(size);
668  }
669
670  if (metric_id == 85)
671  {
673  aspr.ProjectCoefficient(aspr_coeff);
674  DiffuseField(aspr,2);
676  }
677
679  ori.ProjectCoefficient(ori_coeff);
681  target_c = tc;
682  break;
683  }
684  // Targets used for hr-adaptivity tests.
685  case 9: // size target in an annular region.
686  case 10: // size+aspect-ratio in an annular region.
687  case 11: // size+aspect-ratio target for a rotate sine wave
688  {
691  hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
693  target_c = tc;
694  break;
695  }
696  default: cout << "Unknown target_id: " << target_id << endl; return 3;
697  }
698  if (target_c == NULL)
699  {
700  target_c = new TargetConstructor(target_t);
701  }
702  target_c->SetNodes(x0);
703  TMOP_Integrator *tmop_integ = new TMOP_Integrator(metric, target_c,
704  h_metric);
705
706  // Finite differences for computations of derivatives.
707  if (fdscheme)
708  {
709  MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
710  tmop_integ->EnableFiniteDifferences(x);
711  }
712  tmop_integ->SetExactActionFlag(exactaction);
713
714  // Setup the quadrature rules for the TMOP integrator.
715  IntegrationRules *irules = NULL;
717  {
718  case 1: irules = &IntRulesLo; break;
719  case 2: irules = &IntRules; break;
720  case 3: irules = &IntRulesCU; break;
721  default: cout << "Unknown quad_type: " << quad_type << endl; return 3;
722  }
724  if (dim == 2)
725  {
726  cout << "Triangle quadrature points: "
729  << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
730  }
731  if (dim == 3)
732  {
733  cout << "Tetrahedron quadrature points: "
735  << "\nHexahedron quadrature points: "
737  << "\nPrism quadrature points: "
738  << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
739  }
740
741  // Limit the node movement.
742  // The limiting distances can be given by a general function of space.
743  FiniteElementSpace dist_fespace(mesh, fec); // scalar space
744  GridFunction dist(&dist_fespace);
745  dist = 1.0;
746  // The small_phys_size is relevant only with proper normalization.
747  if (normalization) { dist = small_phys_size; }
748  ConstantCoefficient lim_coeff(lim_const);
749  if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
750
756  {
757  MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
758
761
763  else if (adapt_eval == 1)
764  {
765 #ifdef MFEM_USE_GSLIB
767 #else
768  MFEM_ABORT("MFEM is not built with GSLIB support!");
769 #endif
770  }
771  else { MFEM_ABORT("Bad interpolation option."); }
772
775  if (visualization)
776  {
777  socketstream vis1;
778  common::VisualizeField(vis1, "localhost", 19916, adapt_lim_gf0, "Zeta 0",
779  300, 600, 300, 300);
780  }
781  }
782
783  // Surface fitting.
784  L2_FECollection mat_coll(0, dim);
785  H1_FECollection surf_fit_fec(mesh_poly_deg, dim);
786  FiniteElementSpace surf_fit_fes(mesh, &surf_fit_fec);
787  FiniteElementSpace mat_fes(mesh, &mat_coll);
788  GridFunction mat(&mat_fes);
789  GridFunction surf_fit_mat_gf(&surf_fit_fes);
790  GridFunction surf_fit_gf0(&surf_fit_fes);
791  Array<bool> surf_fit_marker(surf_fit_gf0.Size());
792  ConstantCoefficient surf_fit_coeff(surface_fit_const);
794  if (surface_fit_const > 0.0)
795  {
797  "Surface fitting with HR is not implemented yet.");
798  MFEM_VERIFY(pa == false,
799  "Surface fitting with PA is not implemented yet.");
800
802  surf_fit_gf0.ProjectCoefficient(ls_coeff);
803
804  for (int i = 0; i < mesh->GetNE(); i++)
805  {
806  mat(i) = material_id(i, surf_fit_gf0);
807  mesh->SetAttribute(i, mat(i) + 1);
808  }
809
810  GridFunctionCoefficient mat_coeff(&mat);
811  surf_fit_mat_gf.ProjectDiscCoefficient(mat_coeff, GridFunction::ARITHMETIC);
812  for (int j = 0; j < surf_fit_marker.Size(); j++)
813  {
814  if (surf_fit_mat_gf(j) > 0.1 && surf_fit_mat_gf(j) < 0.9)
815  {
816  surf_fit_marker[j] = true;
817  surf_fit_mat_gf(j) = 1.0;
818  }
819  else
820  {
821  surf_fit_marker[j] = false;
822  surf_fit_mat_gf(j) = 0.0;
823  }
824  }
825
827  else if (adapt_eval == 1)
828  {
829 #ifdef MFEM_USE_GSLIB
831 #else
832  MFEM_ABORT("MFEM is not built with GSLIB support!");
833 #endif
834  }
835  else { MFEM_ABORT("Bad interpolation option."); }
836
837  tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker,
839  if (visualization)
840  {
841  socketstream vis1, vis2, vis3;
842  common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0, "Level Set 0",
843  300, 600, 300, 300);
844  common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
845  600, 600, 300, 300);
846  common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
847  "Dofs to Move",
848  900, 600, 300, 300);
849  }
850  }
851
852  // Has to be after the enabling of the limiting / alignment, as it computes
853  // normalization factors for these terms as well.
854  if (normalization) { tmop_integ->EnableNormalization(x0); }
855
856  // 12. Setup the final NonlinearForm (which defines the integral of interest,
857  // its first and second derivatives). Here we can use a combination of
858  // metrics, i.e., optimize the sum of two integrals, where both are
859  // scaled by used-defined space-dependent weights. Note that there are no
860  // command-line options for the weights and the type of the second
861  // metric; one should update those in the code.
862  NonlinearForm a(fespace);
864  ConstantCoefficient *metric_coeff1 = NULL;
865  TMOP_QualityMetric *metric2 = NULL;
866  TargetConstructor *target_c2 = NULL;
867  FunctionCoefficient metric_coeff2(weight_fun);
868
869  // Explicit combination of metrics.
870  if (combomet > 0)
871  {
872  // First metric.
873  metric_coeff1 = new ConstantCoefficient(1.0);
874  tmop_integ->SetCoefficient(*metric_coeff1);
875
876  // Second metric.
877  if (dim == 2) { metric2 = new TMOP_Metric_077; }
878  else { metric2 = new TMOP_Metric_315; }
879  TMOP_Integrator *tmop_integ2 = NULL;
880  if (combomet == 1)
881  {
882  target_c2 = new TargetConstructor(
884  target_c2->SetVolumeScale(0.01);
885  target_c2->SetNodes(x0);
886  tmop_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
887  tmop_integ2->SetCoefficient(metric_coeff2);
888  }
889  else { tmop_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
891  if (fdscheme) { tmop_integ2->EnableFiniteDifferences(x); }
892  tmop_integ2->SetExactActionFlag(exactaction);
893
897  if (normalization) { combo->EnableNormalization(x0); }
898  if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
899
901  }
902  else
903  {
905  }
906
907  if (pa) { a.Setup(); }
908
909  // Compute the minimum det(J) of the starting mesh.
910  tauval = infinity();
911  const int NE = mesh->GetNE();
912  for (int i = 0; i < NE; i++)
913  {
914  const IntegrationRule &ir =
917  for (int j = 0; j < ir.GetNPoints(); j++)
918  {
919  transf->SetIntPoint(&ir.IntPoint(j));
920  tauval = min(tauval, transf->Jacobian().Det());
921  }
922  }
923  cout << "Minimum det(J) of the original mesh is " << tauval << endl;
924
925  if (tauval < 0.0 && metric_id != 22 && metric_id != 211 && metric_id != 252
926  && metric_id != 311 && metric_id != 313 && metric_id != 352)
927  {
928  MFEM_ABORT("The input mesh is inverted! Try an untangling metric.");
929  }
930  if (tauval < 0.0)
931  {
932  MFEM_VERIFY(target_t == TargetConstructor::IDEAL_SHAPE_UNIT_SIZE,
933  "Untangling is supported only for ideal targets.");
934
935  const DenseMatrix &Wideal =
937  tauval /= Wideal.Det();
938
939  // Slightly below minJ0 to avoid div by 0.
940  tauval -= 0.01 * h0.Min();
941  }
942
943  // For HR tests, the energy is normalized by the number of elements.
944  const double init_energy = a.GetGridFunctionEnergy(x) /
945  (hradaptivity ? mesh->GetNE() : 1);
946  double init_metric_energy = init_energy;
947  if (lim_const > 0.0 || adapt_lim_const > 0.0 || surface_fit_const > 0.0)
948  {
949  lim_coeff.constant = 0.0;
951  surf_fit_coeff.constant = 0.0;
952  init_metric_energy = a.GetGridFunctionEnergy(x) /
953  (hradaptivity ? mesh->GetNE() : 1);
954  lim_coeff.constant = lim_const;
956  surf_fit_coeff.constant = surface_fit_const;
957  }
958
959  // Visualize the starting mesh and metric values.
960  // Note that for combinations of metrics, this only shows the first metric.
961  if (visualization)
962  {
963  char title[] = "Initial metric values";
964  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
965  }
966
967  // 13. Fix all boundary nodes, or fix only a given component depending on the
968  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
969  // fixed x/y/z components of the node. Attribute 4 corresponds to an
970  // entirely fixed node. Other boundary attributes do not affect the node
971  // movement boundary conditions.
972  if (move_bnd == false)
973  {
974  Array<int> ess_bdr(mesh->bdr_attributes.Max());
975  ess_bdr = 1;
976  a.SetEssentialBC(ess_bdr);
977  }
978  else
979  {
980  int n = 0;
981  for (int i = 0; i < mesh->GetNBE(); i++)
982  {
983  const int nd = fespace->GetBE(i)->GetDof();
984  const int attr = mesh->GetBdrElement(i)->GetAttribute();
985  MFEM_VERIFY(!(dim == 2 && attr == 3),
986  "Boundary attribute 3 must be used only for 3D meshes. "
987  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
989  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
990  if (attr == 4) { n += nd * dim; }
991  }
992  Array<int> ess_vdofs(n);
993  n = 0;
994  for (int i = 0; i < mesh->GetNBE(); i++)
995  {
996  const int nd = fespace->GetBE(i)->GetDof();
997  const int attr = mesh->GetBdrElement(i)->GetAttribute();
998  fespace->GetBdrElementVDofs(i, vdofs);
999  if (attr == 1) // Fix x components.
1000  {
1001  for (int j = 0; j < nd; j++)
1002  { ess_vdofs[n++] = vdofs[j]; }
1003  }
1004  else if (attr == 2) // Fix y components.
1005  {
1006  for (int j = 0; j < nd; j++)
1007  { ess_vdofs[n++] = vdofs[j+nd]; }
1008  }
1009  else if (attr == 3) // Fix z components.
1010  {
1011  for (int j = 0; j < nd; j++)
1012  { ess_vdofs[n++] = vdofs[j+2*nd]; }
1013  }
1014  else if (attr == 4) // Fix all components.
1015  {
1016  for (int j = 0; j < vdofs.Size(); j++)
1017  { ess_vdofs[n++] = vdofs[j]; }
1018  }
1019  }
1020  a.SetEssentialVDofs(ess_vdofs);
1021  }
1022
1023  // 14. As we use the Newton method to solve the resulting nonlinear system,
1024  // here we setup the linear solver for the system's Jacobian.
1025  Solver *S = NULL, *S_prec = NULL;
1026  const double linsol_rtol = 1e-12;
1027  if (lin_solver == 0)
1028  {
1029  S = new DSmoother(1, 1.0, max_lin_iter);
1030  }
1031  else if (lin_solver == 1)
1032  {
1033  CGSolver *cg = new CGSolver;
1034  cg->SetMaxIter(max_lin_iter);
1035  cg->SetRelTol(linsol_rtol);
1036  cg->SetAbsTol(0.0);
1037  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
1038  S = cg;
1039  }
1040  else
1041  {
1042  MINRESSolver *minres = new MINRESSolver;
1043  minres->SetMaxIter(max_lin_iter);
1044  minres->SetRelTol(linsol_rtol);
1045  minres->SetAbsTol(0.0);
1046  if (verbosity_level > 2) { minres->SetPrintLevel(1); }
1047  minres->SetPrintLevel(verbosity_level == 2 ? 3 : -1);
1048  if (lin_solver == 3 || lin_solver == 4)
1049  {
1050  if (pa)
1051  {
1052  MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
1053  auto js = new OperatorJacobiSmoother;
1054  js->SetPositiveDiagonal(true);
1055  S_prec = js;
1056  }
1057  else
1058  {
1059  auto ds = new DSmoother((lin_solver == 3) ? 0 : 1, 1.0, 1);
1060  ds->SetPositiveDiagonal(true);
1061  S_prec = ds;
1062  }
1063  minres->SetPreconditioner(*S_prec);
1064  }
1065  S = minres;
1066  }
1067
1068  // Perform the nonlinear optimization.
1069  const IntegrationRule &ir =
1071  TMOPNewtonSolver solver(ir, solver_type);
1073  if (surface_fit_threshold > 0)
1074  {
1075  solver.SetTerminationWithMaxSurfaceFittingError(surface_fit_threshold);
1076  }
1077  // Provide all integration rules in case of a mixed mesh.
1079  if (solver_type == 0)
1080  {
1081  // Specify linear solver when we use a Newton-based solver.
1082  solver.SetPreconditioner(*S);
1083  }
1084  // For untangling, the solver will update the min det(T) values.
1085  if (tauval < 0.0) { solver.SetMinDetPtr(&tauval); }
1086  solver.SetMaxIter(solver_iter);
1087  solver.SetRelTol(solver_rtol);
1088  solver.SetAbsTol(0.0);
1089  if (solver_art_type > 0)
1090  {
1092  }
1093  solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
1094
1096  // If hr-adaptivity is disabled, r-adaptivity is done once using the
1097  // TMOPNewtonSolver.
1098  // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
1100  // The solver terminates if an h-adaptivity iteration does not modify
1101  // any element in the mesh.
1102  TMOPHRSolver hr_solver(*mesh, a, solver,
1104  mesh_poly_deg, h_metric_id,
1105  n_hr_iter, n_h_iter);
1108  {
1111  }
1112  hr_solver.Mult();
1113
1114  // 15. Save the optimized mesh to a file. This output can be viewed later
1115  // using GLVis: "glvis -m optimized.mesh".
1116  {
1117  ofstream mesh_ofs("optimized.mesh");
1118  mesh_ofs.precision(14);
1119  mesh->Print(mesh_ofs);
1120  }
1121
1122  const double fin_energy = a.GetGridFunctionEnergy(x) /
1123  (hradaptivity ? mesh->GetNE() : 1);
1124  double fin_metric_energy = fin_energy;
1125  if (lim_const > 0.0 || adapt_lim_const > 0.0)
1126  {
1127  lim_coeff.constant = 0.0;
1129  surf_fit_coeff.constant = 0.0;
1130  fin_metric_energy = a.GetGridFunctionEnergy(x) /
1131  (hradaptivity ? mesh->GetNE() : 1);
1132  lim_coeff.constant = lim_const;
1134  surf_fit_coeff.constant = surface_fit_const;
1135  }
1136  std::cout << std::scientific << std::setprecision(4);
1137  cout << "Initial strain energy: " << init_energy
1138  << " = metrics: " << init_metric_energy
1139  << " + extra terms: " << init_energy - init_metric_energy << endl;
1140  cout << " Final strain energy: " << fin_energy
1141  << " = metrics: " << fin_metric_energy
1142  << " + extra terms: " << fin_energy - fin_metric_energy << endl;
1143  cout << "The strain energy decreased by: "
1144  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
1145
1146  // 16. Visualize the final mesh and metric values.
1147  if (visualization)
1148  {
1149  char title[] = "Final metric values";
1150  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
1151  }
1152
1153  if (adapt_lim_const > 0.0 && visualization)
1154  {
1155  socketstream vis0;
1156  common::VisualizeField(vis0, "localhost", 19916, adapt_lim_gf0, "Xi 0",
1157  600, 600, 300, 300);
1158  }
1159
1160  if (surface_fit_const > 0.0)
1161  {
1162  if (visualization)
1163  {
1164  socketstream vis2, vis3;
1165  common::VisualizeField(vis2, "localhost", 19916, mat, "Materials",
1166  600, 900, 300, 300);
1167  common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf, "Surface dof",
1168  900, 900, 300, 300);
1169  }
1170  double err_avg, err_max;
1171  tmop_integ->GetSurfaceFittingErrors(err_avg, err_max);
1172  std::cout << "Avg fitting error: " << err_avg << std::endl
1173  << "Max fitting error: " << err_max << std::endl;
1174  }
1175
1176  // 17. Visualize the mesh displacement.
1177  if (visualization)
1178  {
1179  osockstream sock(19916, "localhost");
1180  sock << "solution\n";
1181  mesh->Print(sock);
1182  x0 -= x;
1183  x0.Save(sock);
1184  sock.send();
1185  sock << "window_title 'Displacements'\n"
1186  << "window_geometry "
1187  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
1188  << "keys jRmclA" << endl;
1189  }
1190
1191  delete S;
1192  delete S_prec;
1193  delete target_c2;
1194  delete metric2;
1195  delete metric_coeff1;
1198  delete target_c;
1201  delete h_metric;
1202  delete metric;
1203  delete fespace;
1204  delete fec;
1205  delete mesh;
1206
1207  return 0;
1208 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
void discrete_aspr_3d(const Vector &x, Vector &v)
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Definition: solvers.hpp:465
Definition: ex25.cpp:148
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:560
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
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:3642
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:251
double discrete_size_2d(const Vector &x)
3D barrier Shape+Size (VS) metric.
Definition: tmop.hpp:578
Data type for scaled Jacobi-type smoother of sparse matrix.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
void SetPositiveDiagonal(bool pos_diag=true)
Replace diagonal entries with their absolute values.
Definition: solvers.hpp:311
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:192
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:150
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:661
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
Definition: tmop.cpp:1474
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
Definition: tmop.cpp:3480
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:931
3D barrier Shape (S) metric.
Definition: tmop.hpp:484
Definition: tmop.hpp:1233
virtual void Setup()
Setup the NonlinearForm: based on the current AssemblyLevel and the current mesh, optionally...
double Det() const
Definition: densemat.cpp:436
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2569
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:358
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
double surface_level_set(const Vector &x)
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Container class for integration rules.
Definition: intrules.hpp:311
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
3D barrier Shape (S) metric.
Definition: tmop.hpp:466
void GetSurfaceFittingErrors(double &err_avg, double &err_max)
Definition: tmop.cpp:2436
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:994
2D barrier shape (S) metric (polyconvex).
Definition: tmop.hpp:186
Restriction of the node positions to certain regions.
Definition: tmop.cpp:2365
3D non-barrier Size (V) metric.
Definition: tmop.hpp:542
2D barrier size (V) metric (polyconvex).
Definition: tmop.hpp:340
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop.hpp:1556
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
Definition: tmop.cpp:3769
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:774
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
Definition: tmop_tools.cpp:891
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
Definition: tmop.hpp:1711
MINRES method.
Definition: solvers.hpp:575
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:772
virtual int OrderJ() const =0
Return the order of the elements of the Jacobian of the transformation.
void SetMinDetPtr(double *md_ptr)
Definition: tmop_tools.hpp:198
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:204
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1062
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:3670
2D Shifted barrier form of shape metric (mu_2).
Definition: tmop.hpp:253
virtual void SetSerialDiscreteTargetAspectRatio(const GridFunction &tspec_)
Definition: tmop.cpp:1466
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:297
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:323
void SetTerminationWithMaxSurfaceFittingError(double max_error)
Definition: tmop_tools.hpp:219
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3619
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
Geometry Geometries
Definition: fe.cpp:49
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void GetDerivative(int comp, int der_comp, GridFunction &der)
Compute a certain derivative of a function&#39;s component. Derivatives of the function are computed at t...
Definition: gridfunc.cpp:1421
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:588
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:76
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Limiting of the mesh displacements (general version).
Definition: tmop.cpp:2340
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:599
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:1568
2D non-barrier size (V) metric (not polyconvex).
Definition: tmop.hpp:288
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:9394
3D barrier Size (V) metric.
Definition: tmop.hpp:560
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
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:988
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:716
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
double discrete_ori_2d(const Vector &x)
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:396
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5273
Definition: tmop_amr.hpp:253
int Dimension() const
Definition: mesh.hpp:999
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level. The default is AssemblyLevel::LEGACY.
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
Definition: tmop_tools.hpp:215
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:756
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2680
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
A general vector function coefficient.
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:596
void SetAbsTol(double atol)
Definition: solvers.hpp:199
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
Definition: tmop_amr.hpp:252
void SetRelTol(double rtol)
Definition: solvers.hpp:198
2D barrier (not a shape) metric (polyconvex).
Definition: tmop.hpp:272
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: tmop_tools.hpp:249
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
double weight_fun(const Vector &x)
2D barrier size (V) metric (polyconvex).
Definition: tmop.hpp:305
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetAdaptiveLinRtol(const int type=2, const double rtol0=0.5, const double rtol_max=0.9, const double alpha=0.5 *(1.0+sqrt(5.0)), const double gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
Definition: solvers.cpp:1912
Adds a new TMOP_Integrator to the combination.
Definition: tmop.hpp:1735
double material_indicator_2d(const Vector &x)
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:222
double a
Definition: lissajous.cpp:41
double discrete_size_3d(const Vector &x)
double discrete_aspr_2d(const Vector &x)
Class for integration point with weight.
Definition: intrules.hpp:25
double GetGridFunctionEnergy(const Vector &x) const
Compute the enery corresponding to the state x.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1646
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:738
3D Shape (S) metric, untangling version of 303.
Definition: tmop.hpp:521
int dim
Definition: ex24.cpp:53
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2395
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
A general function coefficient.
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1458
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:94
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:640
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:912
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:2400
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Specify essential boundary conditions.
Base class for solvers.
Definition: operator.hpp:651
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 SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:5301
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:908
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:1227
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:617
int material_id(int el_id, const GridFunction &g)
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
Definition: tmop.cpp:1450
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:381
Definition: vector.hpp:458
int main()
double GetElementVolume(int i)
Definition: mesh.cpp:112
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1037
3D barrier Shape (S) metric.
Definition: tmop.hpp:450
2D non-barrier metric without a type.
Definition: tmop.hpp:108
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1303
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:238