MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesh-optimizer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
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:
35 // Adapted analytic shape:
36 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -bnd -qt 1 -qo 8
37 // Adapted analytic size+orientation:
38 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
39 // Adapted analytic shape+orientation:
40 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
41 //
42 // Adapted discrete size:
43 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor
44 // Adapted discrete size 3D with PA:
45 // mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa
46 // Adapted discrete size 3D with PA on device (requires CUDA).
47 // * mesh-optimizer -m cube.mesh -o 3 -rs 3 -mid 321 -tid 5 -ls 3 -nor -lc 0.1 -pa -d cuda
48 // Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
49 // mesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
50 // Adapted discrete size+aspect_ratio:
51 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
52 // mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
53 // Adapted discrete size+orientation (requires GSLIB):
54 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -ae 1 -nor
55 // Adapted discrete aspect-ratio+orientation (requires GSLIB):
56 // * mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 8 -ni 10 -bnd -qt 1 -qo 8 -fd -ae 1
57 // Adapted discrete aspect ratio (3D):
58 // mesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -bnd -qt 1 -qo 8
59 //
60 // Adaptive limiting:
61 // mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
62 // Adaptive limiting through the L-BFGS solver:
63 // mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
64 // Adaptive limiting through FD (requires GSLIB):
65 // * mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
66 //
67 // Blade shape:
68 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
69 // Blade shape with FD-based solver:
70 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
71 // Blade limited shape:
72 // mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
73 // ICF shape and equal size:
74 // mesh-optimizer -o 3 -mid 9 -tid 2 -ni 25 -ls 3 -art 2 -qo 5
75 // ICF shape and initial size:
76 // mesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
77 // ICF shape:
78 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
79 // ICF limited shape:
80 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8 -lc 10
81 // ICF combo shape + size (rings, slow convergence):
82 // mesh-optimizer -o 3 -mid 1 -tid 1 -ni 1000 -bnd -qt 1 -qo 8 -cmb 1
83 // Mixed tet / cube / hex mesh with limiting:
84 // mesh-optimizer -m ../../data/fichera-mixed-p2.mesh -o 4 -rs 1 -mid 301 -tid 1 -fix-bnd -qo 6 -nor -lc 0.25
85 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
86 // * mesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -mid 303 -tid 1 -ni 20 -li 500 -fix-bnd
87 // 2D non-conforming shape and equal size:
88 // mesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -bnd -qt 1 -qo 8
89 //
90 // 2D untangling:
91 // mesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
92 // 3D untangling (the mesh is in the mfem/data GitHub repository):
93 // * 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
94 //
95 
96 #include "mfem.hpp"
97 #include "../common/mfem-common.hpp"
98 #include <fstream>
99 #include <iostream>
100 #include "mesh-optimizer.hpp"
101 
102 using namespace mfem;
103 using namespace std;
104 
105 int main(int argc, char *argv[])
106 {
107  // 0. Set the method's default parameters.
108  const char *mesh_file = "icf.mesh";
109  int mesh_poly_deg = 1;
110  int rs_levels = 0;
111  double jitter = 0.0;
112  int metric_id = 1;
113  int target_id = 1;
114  double lim_const = 0.0;
115  double adapt_lim_const = 0.0;
116  int quad_type = 1;
117  int quad_order = 8;
118  int solver_type = 0;
119  int solver_iter = 20;
120  double solver_rtol = 1e-10;
121  int solver_art_type = 0;
122  int lin_solver = 2;
123  int max_lin_iter = 100;
124  bool move_bnd = true;
125  int combomet = 0;
126  bool normalization = false;
127  bool visualization = true;
128  int verbosity_level = 0;
129  bool fdscheme = false;
130  int adapt_eval = 0;
131  bool exactaction = false;
132  const char *devopt = "cpu";
133  bool pa = false;
134 
135  // 1. Parse command-line options.
136  OptionsParser args(argc, argv);
137  args.AddOption(&mesh_file, "-m", "--mesh",
138  "Mesh file to use.");
139  args.AddOption(&mesh_poly_deg, "-o", "--order",
140  "Polynomial degree of mesh finite element space.");
141  args.AddOption(&rs_levels, "-rs", "--refine-serial",
142  "Number of times to refine the mesh uniformly in serial.");
143  args.AddOption(&jitter, "-ji", "--jitter",
144  "Random perturbation scaling factor.");
145  args.AddOption(&metric_id, "-mid", "--metric-id",
146  "Mesh optimization metric:\n\t"
147  "T-metrics\n\t"
148  "1 : |T|^2 -- 2D shape\n\t"
149  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
150  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
151  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
152  "14 : |T-I|^2 -- 2D shape+size+orientation\n\t"
153  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
154  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
155  "55 : (tau-1)^2 -- 2D size\n\t"
156  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
157  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
158  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
159  "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t"
160  "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t"
161  "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t"
162  // "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
163  // "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
164  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
165  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
166  "303: (|T|^2)/3*tau^(2/3)-1 -- 3D shape\n\t"
167  //"311: (tau-1)^2-tau+sqrt(tau^2+eps)-- 3D untangling\n\t"
168  "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t"
169  "315: (tau-1)^2 -- 3D size\n\t"
170  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D size\n\t"
171  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
172  // "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling\n\t"
173  "A-metrics\n\t"
174  "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t"
175  "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t"
176  "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t"
177  "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t"
178  );
179  args.AddOption(&target_id, "-tid", "--target-id",
180  "Target (ideal element) type:\n\t"
181  "1: Ideal shape, unit size\n\t"
182  "2: Ideal shape, equal size\n\t"
183  "3: Ideal shape, initial size\n\t"
184  "4: Given full analytic Jacobian (in physical space)\n\t"
185  "5: Ideal shape, given size (in physical space)");
186  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
187  args.AddOption(&adapt_lim_const, "-alc", "--adapt-limit-const",
188  "Adaptive limiting coefficient constant.");
189  args.AddOption(&quad_type, "-qt", "--quad-type",
190  "Quadrature rule type:\n\t"
191  "1: Gauss-Lobatto\n\t"
192  "2: Gauss-Legendre\n\t"
193  "3: Closed uniform points");
194  args.AddOption(&quad_order, "-qo", "--quad_order",
195  "Order of the quadrature rule.");
196  args.AddOption(&solver_type, "-st", "--solver-type",
197  " Type of solver: (default) 0: Newton, 1: LBFGS");
198  args.AddOption(&solver_iter, "-ni", "--newton-iters",
199  "Maximum number of Newton iterations.");
200  args.AddOption(&solver_rtol, "-rtol", "--newton-rel-tolerance",
201  "Relative tolerance for the Newton solver.");
202  args.AddOption(&solver_art_type, "-art", "--adaptive-rel-tol",
203  "Type of adaptive relative linear solver tolerance:\n\t"
204  "0: None (default)\n\t"
205  "1: Eisenstat-Walker type 1\n\t"
206  "2: Eisenstat-Walker type 2");
207  args.AddOption(&lin_solver, "-ls", "--lin-solver",
208  "Linear solver:\n\t"
209  "0: l1-Jacobi\n\t"
210  "1: CG\n\t"
211  "2: MINRES\n\t"
212  "3: MINRES + Jacobi preconditioner\n\t"
213  "4: MINRES + l1-Jacobi preconditioner");
214  args.AddOption(&max_lin_iter, "-li", "--lin-iter",
215  "Maximum number of iterations in the linear solve.");
216  args.AddOption(&move_bnd, "-bnd", "--move-boundary", "-fix-bnd",
217  "--fix-boundary",
218  "Enable motion along horizontal and vertical boundaries.");
219  args.AddOption(&combomet, "-cmb", "--combo-type",
220  "Combination of metrics options:\n\t"
221  "0: Use single metric\n\t"
222  "1: Shape + space-dependent size given analytically\n\t"
223  "2: Shape + adapted size given discretely; shared target");
224  args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
225  "--no-normalization",
226  "Make all terms in the optimization functional unitless.");
227  args.AddOption(&fdscheme, "-fd", "--fd_approximation",
228  "-no-fd", "--no-fd-approx",
229  "Enable finite difference based derivative computations.");
230  args.AddOption(&exactaction, "-ex", "--exact_action",
231  "-no-ex", "--no-exact-action",
232  "Enable exact action of TMOP_Integrator.");
233  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
234  "--no-visualization",
235  "Enable or disable GLVis visualization.");
236  args.AddOption(&verbosity_level, "-vl", "--verbosity-level",
237  "Set the verbosity level - 0, 1, or 2.");
238  args.AddOption(&adapt_eval, "-ae", "--adaptivity-evaluator",
239  "0 - Advection based (DEFAULT), 1 - GSLIB.");
240  args.AddOption(&devopt, "-d", "--device",
241  "Device configuration string, see Device::Configure().");
242  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
243  "--no-partial-assembly", "Enable Partial Assembly.");
244  args.Parse();
245  if (!args.Good())
246  {
247  args.PrintUsage(cout);
248  return 1;
249  }
250  args.PrintOptions(cout);
251 
252  Device device(devopt);
253  device.Print();
254 
255  // 2. Initialize and refine the starting mesh.
256  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
257  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
258  const int dim = mesh->Dimension();
259  cout << "Mesh curvature: ";
260  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
261  else { cout << "(NONE)"; }
262  cout << endl;
263 
264  // 3. Define a finite element space on the mesh. Here we use vector finite
265  // elements which are tensor products of quadratic finite elements. The
266  // number of components in the vector finite element space is specified by
267  // the last parameter of the FiniteElementSpace constructor.
269  if (mesh_poly_deg <= 0)
270  {
271  fec = new QuadraticPosFECollection;
272  mesh_poly_deg = 2;
273  }
274  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
275  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim);
276 
277  // 4. Make the mesh curved based on the above finite element space. This
278  // means that we define the mesh elements through a fespace-based
279  // transformation of the reference element.
280  mesh->SetNodalFESpace(fespace);
281 
282  // 5. Set up an empty right-hand side vector b, which is equivalent to b=0.
283  Vector b(0);
284 
285  // 6. Get the mesh nodes (vertices and other degrees of freedom in the finite
286  // element space) as a finite element grid function in fespace. Note that
287  // changing x automatically changes the shapes of the mesh elements.
288  GridFunction x(fespace);
289  mesh->SetNodalGridFunction(&x);
290 
291  // 7. Define a vector representing the minimal local mesh size in the mesh
292  // nodes. We index the nodes using the scalar version of the degrees of
293  // freedom in fespace. Note: this is partition-dependent.
294  //
295  // In addition, compute average mesh size and total volume.
296  Vector h0(fespace->GetNDofs());
297  h0 = infinity();
298  double volume = 0.0;
299  Array<int> dofs;
300  for (int i = 0; i < mesh->GetNE(); i++)
301  {
302  // Get the local scalar element degrees of freedom in dofs.
303  fespace->GetElementDofs(i, dofs);
304  // Adjust the value of h0 in dofs based on the local mesh size.
305  const double hi = mesh->GetElementSize(i);
306  for (int j = 0; j < dofs.Size(); j++)
307  {
308  h0(dofs[j]) = min(h0(dofs[j]), hi);
309  }
310  volume += mesh->GetElementVolume(i);
311  }
312  const double small_phys_size = pow(volume, 1.0 / dim) / 100.0;
313 
314  // 8. Add a random perturbation to the nodes in the interior of the domain.
315  // We define a random grid function of fespace and make sure that it is
316  // zero on the boundary and its values are locally of the order of h0.
317  // The latter is based on the DofToVDof() method which maps the scalar to
318  // the vector degrees of freedom in fespace.
319  GridFunction rdm(fespace);
320  rdm.Randomize();
321  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
322  rdm *= jitter;
323  rdm.HostReadWrite();
324  // Scale the random values to be of order of the local mesh size.
325  for (int i = 0; i < fespace->GetNDofs(); i++)
326  {
327  for (int d = 0; d < dim; d++)
328  {
329  rdm(fespace->DofToVDof(i,d)) *= h0(i);
330  }
331  }
332  Array<int> vdofs;
333  for (int i = 0; i < fespace->GetNBE(); i++)
334  {
335  // Get the vector degrees of freedom in the boundary element.
336  fespace->GetBdrElementVDofs(i, vdofs);
337  // Set the boundary values to zero.
338  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
339  }
340  x -= rdm;
341  x.SetTrueVector();
342  x.SetFromTrueVector();
343 
344  // 9. Save the starting (prior to the optimization) mesh to a file. This
345  // output can be viewed later using GLVis: "glvis -m perturbed.mesh".
346  {
347  ofstream mesh_ofs("perturbed.mesh");
348  mesh->Print(mesh_ofs);
349  }
350 
351  // 10. Store the starting (prior to the optimization) positions.
352  GridFunction x0(fespace);
353  x0 = x;
354 
355  // 11. Form the integrator that uses the chosen metric and target.
356  double tauval = -0.1;
357  TMOP_QualityMetric *metric = NULL;
358  switch (metric_id)
359  {
360  // T-metrics
361  case 1: metric = new TMOP_Metric_001; break;
362  case 2: metric = new TMOP_Metric_002; break;
363  case 7: metric = new TMOP_Metric_007; break;
364  case 9: metric = new TMOP_Metric_009; break;
365  case 14: metric = new TMOP_Metric_014; break;
366  case 22: metric = new TMOP_Metric_022(tauval); break;
367  case 50: metric = new TMOP_Metric_050; break;
368  case 55: metric = new TMOP_Metric_055; break;
369  case 56: metric = new TMOP_Metric_056; break;
370  case 58: metric = new TMOP_Metric_058; break;
371  case 77: metric = new TMOP_Metric_077; break;
372  case 80: metric = new TMOP_Metric_080(0.5); break;
373  case 85: metric = new TMOP_Metric_085; break;
374  case 98: metric = new TMOP_Metric_098; break;
375  // case 211: metric = new TMOP_Metric_211; break;
376  // case 252: metric = new TMOP_Metric_252(tauval); break;
377  case 301: metric = new TMOP_Metric_301; break;
378  case 302: metric = new TMOP_Metric_302; break;
379  case 303: metric = new TMOP_Metric_303; break;
380  // case 311: metric = new TMOP_Metric_311; break;
381  case 313: metric = new TMOP_Metric_313(tauval); break;
382  case 315: metric = new TMOP_Metric_315; break;
383  case 316: metric = new TMOP_Metric_316; break;
384  case 321: metric = new TMOP_Metric_321; break;
385  // case 352: metric = new TMOP_Metric_352(tauval); break;
386  // A-metrics
387  case 11: metric = new TMOP_AMetric_011; break;
388  case 36: metric = new TMOP_AMetric_036; break;
389  case 107: metric = new TMOP_AMetric_107a; break;
390  case 126: metric = new TMOP_AMetric_126(0.9); break;
391  default:
392  cout << "Unknown metric_id: " << metric_id << endl;
393  return 3;
394  }
396  TargetConstructor *target_c = NULL;
397  HessianCoefficient *adapt_coeff = NULL;
398  H1_FECollection ind_fec(mesh_poly_deg, dim);
399  FiniteElementSpace ind_fes(mesh, &ind_fec);
400  FiniteElementSpace ind_fesv(mesh, &ind_fec, dim);
401  GridFunction size(&ind_fes), aspr(&ind_fes), disc(&ind_fes), ori(&ind_fes);
402  GridFunction aspr3d(&ind_fesv);
403 
404  const AssemblyLevel al =
406 
407  switch (target_id)
408  {
409  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
410  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
411  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
412  case 4: // Analytic
413  {
415  AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
416  adapt_coeff = new HessianCoefficient(dim, metric_id);
417  tc->SetAnalyticTargetSpec(NULL, NULL, adapt_coeff);
418  target_c = tc;
419  break;
420  }
421  case 5: // Discrete size 2D or 3D
422  {
424  DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
425  if (adapt_eval == 0)
426  {
427  tc->SetAdaptivityEvaluator(new AdvectorCG(al));
428  }
429  else
430  {
431 #ifdef MFEM_USE_GSLIB
433 #else
434  MFEM_ABORT("MFEM is not built with GSLIB.");
435 #endif
436  }
437  if (dim == 2)
438  {
440  size.ProjectCoefficient(ind_coeff);
441  }
442  else if (dim == 3)
443  {
445  size.ProjectCoefficient(ind_coeff);
446  }
447  tc->SetSerialDiscreteTargetSize(size);
448  target_c = tc;
449  break;
450  }
451  case 6: // Discrete size + aspect ratio - 2D
452  {
453  GridFunction d_x(&ind_fes), d_y(&ind_fes);
454 
456  DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
458  disc.ProjectCoefficient(ind_coeff);
459  if (adapt_eval == 0)
460  {
461  tc->SetAdaptivityEvaluator(new AdvectorCG(al));
462  }
463  else
464  {
465 #ifdef MFEM_USE_GSLIB
467 #else
468  MFEM_ABORT("MFEM is not built with GSLIB.");
469 #endif
470  }
471 
472  // Diffuse the interface
473  DiffuseField(disc,2);
474 
475  // Get partials with respect to x and y of the grid function
476  disc.GetDerivative(1,0,d_x);
477  disc.GetDerivative(1,1,d_y);
478 
479  // Compute the squared magnitude of the gradient
480  for (int i = 0; i < size.Size(); i++)
481  {
482  size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
483  }
484  const double max = size.Max();
485 
486  for (int i = 0; i < d_x.Size(); i++)
487  {
488  d_x(i) = std::abs(d_x(i));
489  d_y(i) = std::abs(d_y(i));
490  }
491  const double eps = 0.01;
492  const double aspr_ratio = 20.0;
493  const double size_ratio = 40.0;
494 
495  for (int i = 0; i < size.Size(); i++)
496  {
497  size(i) = (size(i)/max);
498  aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
499  aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
500  if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
501  if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
502  }
503  Vector vals;
504  const int NE = mesh->GetNE();
505  double volume = 0.0, volume_ind = 0.0;
506 
507  for (int i = 0; i < NE; i++)
508  {
510  const IntegrationRule &ir =
511  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr->OrderJ());
512  size.GetValues(i, ir, vals);
513  for (int j = 0; j < ir.GetNPoints(); j++)
514  {
515  const IntegrationPoint &ip = ir.IntPoint(j);
516  Tr->SetIntPoint(&ip);
517  volume += ip.weight * Tr->Weight();
518  volume_ind += vals(j) * ip.weight * Tr->Weight();
519  }
520  }
521 
522  const double avg_zone_size = volume / NE;
523 
524  const double small_avg_ratio = (volume_ind + (volume - volume_ind) /
525  size_ratio) /
526  volume;
527 
528  const double small_zone_size = small_avg_ratio * avg_zone_size;
529  const double big_zone_size = size_ratio * small_zone_size;
530 
531  for (int i = 0; i < size.Size(); i++)
532  {
533  const double val = size(i);
534  const double a = (big_zone_size - small_zone_size) / small_zone_size;
535  size(i) = big_zone_size / (1.0+a*val);
536  }
537 
538  DiffuseField(size, 2);
539  DiffuseField(aspr, 2);
540 
541  tc->SetSerialDiscreteTargetSize(size);
543  target_c = tc;
544  break;
545  }
546  case 7: // Discrete aspect ratio 3D
547  {
549  DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
550  if (adapt_eval == 0)
551  {
552  tc->SetAdaptivityEvaluator(new AdvectorCG(al));
553  }
554  else
555  {
556 #ifdef MFEM_USE_GSLIB
558 #else
559  MFEM_ABORT("MFEM is not built with GSLIB.");
560 #endif
561  }
563  aspr3d.ProjectCoefficient(fd_aspr3d);
564 
566  target_c = tc;
567  break;
568  }
569  case 8: // shape/size + orientation 2D
570  {
572  DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
573  if (adapt_eval == 0)
574  {
575  tc->SetAdaptivityEvaluator(new AdvectorCG(al));
576  }
577  else
578  {
579 #ifdef MFEM_USE_GSLIB
581 #else
582  MFEM_ABORT("MFEM is not built with GSLIB.");
583 #endif
584  }
585 
586  if (metric_id == 14 || metric_id == 36)
587  {
588  ConstantCoefficient ind_coeff(0.1*0.1);
589  size.ProjectCoefficient(ind_coeff);
590  tc->SetSerialDiscreteTargetSize(size);
591  }
592 
593  if (metric_id == 85)
594  {
596  aspr.ProjectCoefficient(aspr_coeff);
597  DiffuseField(aspr,2);
599  }
600 
602  ori.ProjectCoefficient(ori_coeff);
604  target_c = tc;
605  break;
606  }
607  default: cout << "Unknown target_id: " << target_id << endl; return 3;
608  }
609  if (target_c == NULL)
610  {
611  target_c = new TargetConstructor(target_t);
612  }
613  target_c->SetNodes(x0);
614  TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c);
615 
616  // Finite differences for computations of derivatives.
617  if (fdscheme)
618  {
619  MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
620  he_nlf_integ->EnableFiniteDifferences(x);
621  }
622  he_nlf_integ->SetExactActionFlag(exactaction);
623 
624  // Setup the quadrature rules for the TMOP integrator.
625  IntegrationRules *irules = NULL;
626  switch (quad_type)
627  {
628  case 1: irules = &IntRulesLo; break;
629  case 2: irules = &IntRules; break;
630  case 3: irules = &IntRulesCU; break;
631  default: cout << "Unknown quad_type: " << quad_type << endl; return 3;
632  }
633  he_nlf_integ->SetIntegrationRules(*irules, quad_order);
634  if (dim == 2)
635  {
636  cout << "Triangle quadrature points: "
637  << irules->Get(Geometry::TRIANGLE, quad_order).GetNPoints()
638  << "\nQuadrilateral quadrature points: "
639  << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
640  }
641  if (dim == 3)
642  {
643  cout << "Tetrahedron quadrature points: "
644  << irules->Get(Geometry::TETRAHEDRON, quad_order).GetNPoints()
645  << "\nHexahedron quadrature points: "
646  << irules->Get(Geometry::CUBE, quad_order).GetNPoints()
647  << "\nPrism quadrature points: "
648  << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
649  }
650 
651  if (normalization) { he_nlf_integ->EnableNormalization(x0); }
652 
653  // Limit the node movement.
654  // The limiting distances can be given by a general function of space.
655  FiniteElementSpace dist_fespace(mesh, fec); // scalar space
656  GridFunction dist(&dist_fespace);
657  dist = 1.0;
658  // The small_phys_size is relevant only with proper normalization.
659  if (normalization) { dist = small_phys_size; }
660  ConstantCoefficient lim_coeff(lim_const);
661  if (lim_const != 0.0) { he_nlf_integ->EnableLimiting(x0, dist, lim_coeff); }
662 
663  // Adaptive limiting.
664  GridFunction zeta_0(&ind_fes);
665  ConstantCoefficient coef_zeta(adapt_lim_const);
666  AdaptivityEvaluator *adapt_evaluator = NULL;
667  if (adapt_lim_const > 0.0)
668  {
669  MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
670 
672  zeta_0.ProjectCoefficient(alim_coeff);
673 
674  if (adapt_eval == 0) { adapt_evaluator = new AdvectorCG(al); }
675  else if (adapt_eval == 1)
676  {
677 #ifdef MFEM_USE_GSLIB
678  adapt_evaluator = new InterpolatorFP;
679 #else
680  MFEM_ABORT("MFEM is not built with GSLIB support!");
681 #endif
682  }
683  else { MFEM_ABORT("Bad interpolation option."); }
684 
685  he_nlf_integ->EnableAdaptiveLimiting(zeta_0, coef_zeta, *adapt_evaluator);
686  if (visualization)
687  {
688  socketstream vis1;
689  common::VisualizeField(vis1, "localhost", 19916, zeta_0, "Zeta 0",
690  300, 600, 300, 300);
691  }
692  }
693 
694  // 12. Setup the final NonlinearForm (which defines the integral of interest,
695  // its first and second derivatives). Here we can use a combination of
696  // metrics, i.e., optimize the sum of two integrals, where both are
697  // scaled by used-defined space-dependent weights. Note that there are no
698  // command-line options for the weights and the type of the second
699  // metric; one should update those in the code.
700  NonlinearForm a(fespace);
702  ConstantCoefficient *coeff1 = NULL;
703  TMOP_QualityMetric *metric2 = NULL;
704  TargetConstructor *target_c2 = NULL;
706 
707  // Explicit combination of metrics.
708  if (combomet > 0)
709  {
710  // First metric.
711  coeff1 = new ConstantCoefficient(1.0);
712  he_nlf_integ->SetCoefficient(*coeff1);
713 
714  // Second metric.
715  if (dim == 2) { metric2 = new TMOP_Metric_077; }
716  else { metric2 = new TMOP_Metric_315; }
717  TMOP_Integrator *he_nlf_integ2 = NULL;
718  if (combomet == 1)
719  {
720  target_c2 = new TargetConstructor(
722  target_c2->SetVolumeScale(0.01);
723  target_c2->SetNodes(x0);
724  he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
725  he_nlf_integ2->SetCoefficient(coeff2);
726  }
727  else { he_nlf_integ2 = new TMOP_Integrator(metric2, target_c); }
728  he_nlf_integ2->SetIntegrationRules(*irules, quad_order);
729  if (fdscheme) { he_nlf_integ2->EnableFiniteDifferences(x); }
730  he_nlf_integ2->SetExactActionFlag(exactaction);
731 
733  combo->AddTMOPIntegrator(he_nlf_integ);
734  combo->AddTMOPIntegrator(he_nlf_integ2);
735  if (normalization) { combo->EnableNormalization(x0); }
736  if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
737 
738  a.AddDomainIntegrator(combo);
739  }
740  else
741  {
742  a.AddDomainIntegrator(he_nlf_integ);
743  }
744 
745  if (pa) { a.Setup(); }
746 
747  // Compute the minimum det(J) of the starting mesh.
748  tauval = infinity();
749  const int NE = mesh->GetNE();
750  for (int i = 0; i < NE; i++)
751  {
752  const IntegrationRule &ir =
753  irules->Get(fespace->GetFE(i)->GetGeomType(), quad_order);
755  for (int j = 0; j < ir.GetNPoints(); j++)
756  {
757  transf->SetIntPoint(&ir.IntPoint(j));
758  tauval = min(tauval, transf->Jacobian().Det());
759  }
760  }
761  cout << "Minimum det(J) of the original mesh is " << tauval << endl;
762 
763  if (tauval < 0.0 && metric_id != 22 && metric_id != 211 && metric_id != 252
764  && metric_id != 311 && metric_id != 313 && metric_id != 352)
765  {
766  MFEM_ABORT("The input mesh is inverted! Try an untangling metric.");
767  }
768  if (tauval < 0.0)
769  {
770  MFEM_VERIFY(target_t == TargetConstructor::IDEAL_SHAPE_UNIT_SIZE,
771  "Untangling is supported only for ideal targets.");
772 
773  const DenseMatrix &Wideal =
775  tauval /= Wideal.Det();
776 
777  // Slightly below minJ0 to avoid div by 0.
778  tauval -= 0.01 * h0.Min();
779  }
780 
781  const double init_energy = a.GetGridFunctionEnergy(x);
782 
783  // Visualize the starting mesh and metric values.
784  // Note that for combinations of metrics, this only shows the first metric.
785  if (visualization)
786  {
787  char title[] = "Initial metric values";
788  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
789  }
790 
791  // 13. Fix all boundary nodes, or fix only a given component depending on the
792  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
793  // fixed x/y/z components of the node. Attribute 4 corresponds to an
794  // entirely fixed node. Other boundary attributes do not affect the node
795  // movement boundary conditions.
796  if (move_bnd == false)
797  {
798  Array<int> ess_bdr(mesh->bdr_attributes.Max());
799  ess_bdr = 1;
800  a.SetEssentialBC(ess_bdr);
801  }
802  else
803  {
804  int n = 0;
805  for (int i = 0; i < mesh->GetNBE(); i++)
806  {
807  const int nd = fespace->GetBE(i)->GetDof();
808  const int attr = mesh->GetBdrElement(i)->GetAttribute();
809  MFEM_VERIFY(!(dim == 2 && attr == 3),
810  "Boundary attribute 3 must be used only for 3D meshes. "
811  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
812  "components, rest for free nodes), or use -fix-bnd.");
813  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
814  if (attr == 4) { n += nd * dim; }
815  }
816  Array<int> ess_vdofs(n), vdofs;
817  n = 0;
818  for (int i = 0; i < mesh->GetNBE(); i++)
819  {
820  const int nd = fespace->GetBE(i)->GetDof();
821  const int attr = mesh->GetBdrElement(i)->GetAttribute();
822  fespace->GetBdrElementVDofs(i, vdofs);
823  if (attr == 1) // Fix x components.
824  {
825  for (int j = 0; j < nd; j++)
826  { ess_vdofs[n++] = vdofs[j]; }
827  }
828  else if (attr == 2) // Fix y components.
829  {
830  for (int j = 0; j < nd; j++)
831  { ess_vdofs[n++] = vdofs[j+nd]; }
832  }
833  else if (attr == 3) // Fix z components.
834  {
835  for (int j = 0; j < nd; j++)
836  { ess_vdofs[n++] = vdofs[j+2*nd]; }
837  }
838  else if (attr == 4) // Fix all components.
839  {
840  for (int j = 0; j < vdofs.Size(); j++)
841  { ess_vdofs[n++] = vdofs[j]; }
842  }
843  }
844  a.SetEssentialVDofs(ess_vdofs);
845  }
846 
847  // 14. As we use the Newton method to solve the resulting nonlinear system,
848  // here we setup the linear solver for the system's Jacobian.
849  Solver *S = NULL, *S_prec = NULL;
850  const double linsol_rtol = 1e-12;
851  if (lin_solver == 0)
852  {
853  S = new DSmoother(1, 1.0, max_lin_iter);
854  }
855  else if (lin_solver == 1)
856  {
857  CGSolver *cg = new CGSolver;
858  cg->SetMaxIter(max_lin_iter);
859  cg->SetRelTol(linsol_rtol);
860  cg->SetAbsTol(0.0);
861  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
862  S = cg;
863  }
864  else
865  {
866  MINRESSolver *minres = new MINRESSolver;
867  minres->SetMaxIter(max_lin_iter);
868  minres->SetRelTol(linsol_rtol);
869  minres->SetAbsTol(0.0);
870  if (verbosity_level > 2) { minres->SetPrintLevel(1); }
871  minres->SetPrintLevel(verbosity_level == 2 ? 3 : -1);
872  if (lin_solver == 3 || lin_solver == 4)
873  {
874  if (pa)
875  {
876  MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
877  S_prec = new OperatorJacobiSmoother;
878  }
879  else
880  {
881  S_prec = new DSmoother((lin_solver == 3) ? 0 : 1, 1.0, 1);
882  }
883  minres->SetPreconditioner(*S_prec);
884  }
885  S = minres;
886  }
887 
888  // Perform the nonlinear optimization.
889  const IntegrationRule &ir =
890  irules->Get(fespace->GetFE(0)->GetGeomType(), quad_order);
891  TMOPNewtonSolver solver(ir, solver_type);
892  // Provide all integration rules in case of a mixed mesh.
893  solver.SetIntegrationRules(*irules, quad_order);
894  if (solver_type == 0)
895  {
896  // Specify linear solver when we use a Newton-based solver.
897  solver.SetPreconditioner(*S);
898  }
899  // For untangling, the solver will update the min det(T) values.
900  if (tauval < 0.0) { solver.SetMinDetPtr(&tauval); }
901  solver.SetMaxIter(solver_iter);
902  solver.SetRelTol(solver_rtol);
903  solver.SetAbsTol(0.0);
904  if (solver_art_type > 0)
905  {
906  solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
907  }
908  solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
909  solver.SetOperator(a);
910  solver.Mult(b, x.GetTrueVector());
911  x.SetFromTrueVector();
912  if (solver.GetConverged() == false)
913  {
914  cout << "Nonlinear solver: rtol = " << solver_rtol << " not achieved.\n";
915  }
916 
917  // 15. Save the optimized mesh to a file. This output can be viewed later
918  // using GLVis: "glvis -m optimized.mesh".
919  {
920  ofstream mesh_ofs("optimized.mesh");
921  mesh_ofs.precision(14);
922  mesh->Print(mesh_ofs);
923  }
924 
925  // 16. Compute the amount of energy decrease.
926  const double fin_energy = a.GetGridFunctionEnergy(x);
927  double metric_part = fin_energy;
928  if (lim_const > 0.0 || adapt_lim_const > 0.0)
929  {
930  lim_coeff.constant = 0.0;
931  coef_zeta.constant = 0.0;
932  metric_part = a.GetGridFunctionEnergy(x);
933  lim_coeff.constant = lim_const;
934  coef_zeta.constant = adapt_lim_const;
935  }
936  cout << "Initial strain energy: " << init_energy
937  << " = metrics: " << init_energy
938  << " + limiting term: " << 0.0 << endl;
939  cout << " Final strain energy: " << fin_energy
940  << " = metrics: " << metric_part
941  << " + limiting term: " << fin_energy - metric_part << endl;
942  cout << "The strain energy decreased by: " << setprecision(12)
943  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
944 
945  // 17. Visualize the final mesh and metric values.
946  if (visualization)
947  {
948  char title[] = "Final metric values";
949  vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
950  }
951 
952  if (adapt_lim_const > 0.0 && visualization)
953  {
954  socketstream vis0;
955  common::VisualizeField(vis0, "localhost", 19916, zeta_0, "Xi 0",
956  600, 600, 300, 300);
957  }
958 
959  // 18. Visualize the mesh displacement.
960  if (visualization)
961  {
962  x0 -= x;
963  osockstream sock(19916, "localhost");
964  sock << "solution\n";
965  mesh->Print(sock);
966  x0.Save(sock);
967  sock.send();
968  sock << "window_title 'Displacements'\n"
969  << "window_geometry "
970  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
971  << "keys jRmclA" << endl;
972  }
973 
974  // 19. Free the used memory.
975  delete S;
976  delete S_prec;
977  delete target_c2;
978  delete metric2;
979  delete coeff1;
980  delete adapt_evaluator;
981  delete target_c;
982  delete adapt_coeff;
983  delete metric;
984  delete fespace;
985  delete fec;
986  delete mesh;
987 
988  return 0;
989 }
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:134
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1387
Conjugate gradient method.
Definition: solvers.hpp:316
Definition: ex25.cpp:148
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:537
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:3070
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:233
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:920
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:172
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:143
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
Definition: tmop.cpp:1454
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
Definition: tmop.cpp:2938
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:849
3D barrier Shape (S) metric.
Definition: tmop.hpp:484
void SetAdaptivityEvaluator(AdaptivityEvaluator *ae)
Definition: tmop.hpp:1158
virtual void Setup()
Setup the NonlinearForm: based on the current AssemblyLevel and the current mesh, optionally...
double Det() const
Definition: densemat.cpp:451
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:358
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
Container class for integration rules.
Definition: intrules.hpp:311
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
3D barrier Shape (S) metric.
Definition: tmop.hpp:466
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:950
2D barrier shape (S) metric (polyconvex).
Definition: tmop.hpp:186
void EnableAdaptiveLimiting(const GridFunction &z0, Coefficient &coeff, AdaptivityEvaluator &ae)
Restriction of the node positions to certain regions.
Definition: tmop.cpp:2239
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:1418
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
Definition: tmop.cpp:3171
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:734
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:733
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
Definition: tmop.hpp:1523
virtual void GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2298
MINRES method.
Definition: solvers.hpp:426
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:763
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:178
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:204
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:971
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:3098
2D Shifted barrier form of shape metric (mu_2).
Definition: tmop.hpp:253
virtual void SetSerialDiscreteTargetAspectRatio(const GridFunction &tspec_)
Definition: tmop.cpp:1445
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:97
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:323
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3484
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.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:13507
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:137
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:439
void SetPrintLevel(int print_lvl)
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:2214
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:576
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Definition: tmop.hpp:1430
2D non-barrier size (V) metric (not polyconvex).
Definition: tmop.hpp:288
3D barrier Size (V) metric.
Definition: tmop.hpp:560
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:130
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:111
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:944
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
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:579
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
double discrete_ori_2d(const Vector &x)
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:396
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4843
int GetConverged() const
Definition: solvers.hpp:104
int Dimension() const
Definition: mesh.hpp:911
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:457
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:716
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:127
A general vector function coefficient.
void SetAbsTol(double atol)
Definition: solvers.hpp:99
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:204
void SetRelTol(double rtol)
Definition: solvers.hpp:98
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:212
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp: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:1745
double adapt_lim_fun(const Vector &x)
void AddTMOPIntegrator(TMOP_Integrator *ti)
Adds a new TMOP_Integrator to the combination.
Definition: tmop.hpp:1541
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:347
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:698
3D Shape (S) metric, untangling version of 303.
Definition: tmop.hpp:521
int dim
Definition: ex24.cpp:53
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
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:2388
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
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:2278
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
A general function coefficient.
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: tmop_tools.hpp:187
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
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:872
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Specify essential boundary conditions.
Base class for solvers.
Definition: operator.hpp:648
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:2686
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:268
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:4871
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:868
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
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
Definition: tmop.cpp:1428
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:377
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:381
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:446
int main()
double GetElementVolume(int i)
Definition: mesh.cpp:112
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:946
3D barrier Shape (S) metric.
Definition: tmop.hpp:450
2D non-barrier metric without a type.
Definition: tmop.hpp:108
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1642
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:1198
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:238