MFEM  v4.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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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 // Blade shape:
36 // mesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
37 // Blade limited shape:
38 // mesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8 -lc 5000
39 // ICF shape and equal size:
40 // mesh-optimizer -o 3 -rs 0 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
41 // ICF shape and initial size:
42 // mesh-optimizer -o 3 -rs 0 -mid 9 -tid 3 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8
43 // ICF shape:
44 // mesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8
45 // ICF limited shape:
46 // mesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -lc 10
47 // ICF combo shape + size (rings, slow convergence):
48 // mesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 1000 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb
49 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
50 // * mesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -rs 0 -mid 303 -tid 1 -ni 20 -ls 2 -li 500 -fix-bnd
51 
52 #include "mfem.hpp"
53 #include <fstream>
54 #include <iostream>
55 
56 using namespace mfem;
57 using namespace std;
58 
59 double weight_fun(const Vector &x);
60 
61 // Metric values are visualized by creating an L2 finite element functions and
62 // computing the metric values at the nodes.
63 void vis_metric(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc,
64  Mesh &mesh, char *title, int position)
65 {
67  FiniteElementSpace fes(&mesh, &fec, 1);
68  GridFunction metric(&fes);
69  InterpolateTMOP_QualityMetric(qm, tc, mesh, metric);
70  osockstream sock(19916, "localhost");
71  sock << "solution\n";
72  mesh.Print(sock);
73  metric.Save(sock);
74  sock.send();
75  sock << "window_title '"<< title << "'\n"
76  << "window_geometry "
77  << position << " " << 0 << " " << 600 << " " << 600 << "\n"
78  << "keys jRmclA" << endl;
79 }
80 
81 class RelaxedNewtonSolver : public NewtonSolver
82 {
83 private:
84  // Quadrature points that are checked for negative Jacobians etc.
85  const IntegrationRule &ir;
86  FiniteElementSpace *fes;
87  mutable GridFunction x_gf;
88 
89 public:
90  RelaxedNewtonSolver(const IntegrationRule &irule, FiniteElementSpace *f)
91  : ir(irule), fes(f) { }
92 
93  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
94 };
95 
96 double RelaxedNewtonSolver::ComputeScalingFactor(const Vector &x,
97  const Vector &b) const
98 {
99  const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
100  MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
101  const bool have_b = (b.Size() == Height());
102 
103  const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
104  dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
105  Array<int> xdofs(dof * dim);
106  DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
107  Vector posV(pos.Data(), dof * dim);
108 
109  Vector x_out(x.Size());
110  bool x_out_ok = false;
111  const double energy_in = nlf->GetEnergy(x);
112  double scale = 1.0, energy_out;
113  double norm0 = Norm(r);
114  x_gf.MakeTRef(fes, x_out, 0);
115 
116  // Decreases the scaling of the update until the new mesh is valid.
117  for (int i = 0; i < 12; i++)
118  {
119  add(x, -scale, c, x_out);
120  x_gf.SetFromTrueVector();
121 
122  energy_out = nlf->GetGridFunctionEnergy(x_gf);
123  if (energy_out > 1.2*energy_in || std::isnan(energy_out) != 0)
124  {
125  if (print_level >= 0)
126  { cout << "Scale = " << scale << " Increasing energy." << endl; }
127  scale *= 0.5; continue;
128  }
129 
130  int jac_ok = 1;
131  for (int i = 0; i < NE; i++)
132  {
133  fes->GetElementVDofs(i, xdofs);
134  x_gf.GetSubVector(xdofs, posV);
135  for (int j = 0; j < nsp; j++)
136  {
137  fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
138  MultAtB(pos, dshape, Jpr);
139  if (Jpr.Det() <= 0.0) { jac_ok = 0; goto break2; }
140  }
141  }
142  break2:
143  if (jac_ok == 0)
144  {
145  if (print_level >= 0)
146  { cout << "Scale = " << scale << " Neg det(J) found." << endl; }
147  scale *= 0.5; continue;
148  }
149 
150  oper->Mult(x_out, r);
151  if (have_b) { r -= b; }
152  double norm = Norm(r);
153 
154  if (norm > 1.2*norm0)
155  {
156  if (print_level >= 0)
157  { cout << "Scale = " << scale << " Norm increased." << endl; }
158  scale *= 0.5; continue;
159  }
160  else { x_out_ok = true; break; }
161  }
162 
163  if (print_level >= 0)
164  {
165  cout << "Energy decrease: "
166  << (energy_in - energy_out) / energy_in * 100.0
167  << "% with " << scale << " scaling." << endl;
168  }
169 
170  if (x_out_ok == false) { scale = 0.0; }
171 
172  return scale;
173 }
174 
175 // Allows negative Jacobians. Used in untangling metrics.
176 class DescentNewtonSolver : public NewtonSolver
177 {
178 private:
179  // Quadrature points that are checked for negative Jacobians etc.
180  const IntegrationRule &ir;
181  FiniteElementSpace *fes;
182  mutable GridFunction x_gf;
183 
184 public:
185  DescentNewtonSolver(const IntegrationRule &irule, FiniteElementSpace *f)
186  : ir(irule), fes(f) { }
187 
188  virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const;
189 };
190 
191 double DescentNewtonSolver::ComputeScalingFactor(const Vector &x,
192  const Vector &b) const
193 {
194  const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
195  MFEM_VERIFY(nlf != NULL, "invalid Operator subclass");
196 
197  const int NE = fes->GetMesh()->GetNE(), dim = fes->GetFE(0)->GetDim(),
198  dof = fes->GetFE(0)->GetDof(), nsp = ir.GetNPoints();
199  Array<int> xdofs(dof * dim);
200  DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
201  Vector posV(pos.Data(), dof * dim);
202 
203  x_gf.MakeTRef(fes, x.GetData());
204  x_gf.SetFromTrueVector();
205 
206  double min_detJ = infinity();
207  for (int i = 0; i < NE; i++)
208  {
209  fes->GetElementVDofs(i, xdofs);
210  x_gf.GetSubVector(xdofs, posV);
211  for (int j = 0; j < nsp; j++)
212  {
213  fes->GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
214  MultAtB(pos, dshape, Jpr);
215  min_detJ = min(min_detJ, Jpr.Det());
216  }
217  }
218  cout << "Minimum det(J) = " << min_detJ << endl;
219 
220  Vector x_out(x.Size());
221  bool x_out_ok = false;
222  const double energy_in = nlf->GetGridFunctionEnergy(x_gf);
223  double scale = 1.0, energy_out;
224 
225  for (int i = 0; i < 7; i++)
226  {
227  add(x, -scale, c, x_out);
228 
229  energy_out = nlf->GetEnergy(x_out);
230  if (energy_out > energy_in || std::isnan(energy_out) != 0)
231  {
232  scale *= 0.5;
233  }
234  else { x_out_ok = true; break; }
235  }
236 
237  cout << "Energy decrease: " << (energy_in - energy_out) / energy_in * 100.0
238  << "% with " << scale << " scaling." << endl;
239 
240  if (x_out_ok == false) { return 0.0; }
241 
242  return scale;
243 }
244 
245 // Additional IntegrationRules that can be used with the --quad-type option.
248 
249 
250 int main (int argc, char *argv[])
251 {
252  // 0. Set the method's default parameters.
253  const char *mesh_file = "icf.mesh";
254  int mesh_poly_deg = 1;
255  int rs_levels = 0;
256  double jitter = 0.0;
257  int metric_id = 1;
258  int target_id = 1;
259  double lim_const = 0.0;
260  int quad_type = 1;
261  int quad_order = 8;
262  int newton_iter = 10;
263  double newton_rtol = 1e-12;
264  int lin_solver = 2;
265  int max_lin_iter = 100;
266  bool move_bnd = true;
267  bool combomet = 0;
268  bool normalization = false;
269  bool visualization = true;
270  int verbosity_level = 0;
271 
272  // 1. Parse command-line options.
273  OptionsParser args(argc, argv);
274  args.AddOption(&mesh_file, "-m", "--mesh",
275  "Mesh file to use.");
276  args.AddOption(&mesh_poly_deg, "-o", "--order",
277  "Polynomial degree of mesh finite element space.");
278  args.AddOption(&rs_levels, "-rs", "--refine-serial",
279  "Number of times to refine the mesh uniformly in serial.");
280  args.AddOption(&jitter, "-ji", "--jitter",
281  "Random perturbation scaling factor.");
282  args.AddOption(&metric_id, "-mid", "--metric-id",
283  "Mesh optimization metric:\n\t"
284  "1 : |T|^2 -- 2D shape\n\t"
285  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
286  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
287  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
288  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
289  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
290  "55 : (tau-1)^2 -- 2D size\n\t"
291  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
292  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
293  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
294  "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
295  "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
296  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
297  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
298  "303: (|T|^2)/3*tau^(2/3)-1 -- 3D shape\n\t"
299  "315: (tau-1)^2 -- 3D size\n\t"
300  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D size\n\t"
301  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
302  "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling");
303  args.AddOption(&target_id, "-tid", "--target-id",
304  "Target (ideal element) type:\n\t"
305  "1: Ideal shape, unit size\n\t"
306  "2: Ideal shape, equal size\n\t"
307  "3: Ideal shape, initial size");
308  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
309  args.AddOption(&quad_type, "-qt", "--quad-type",
310  "Quadrature rule type:\n\t"
311  "1: Gauss-Lobatto\n\t"
312  "2: Gauss-Legendre\n\t"
313  "3: Closed uniform points");
314  args.AddOption(&quad_order, "-qo", "--quad_order",
315  "Order of the quadrature rule.");
316  args.AddOption(&newton_iter, "-ni", "--newton-iters",
317  "Maximum number of Newton iterations.");
318  args.AddOption(&newton_rtol, "-rtol", "--newton-rel-tolerance",
319  "Relative tolerance for the Newton solver.");
320  args.AddOption(&lin_solver, "-ls", "--lin-solver",
321  "Linear solver: 0 - l1-Jacobi, 1 - CG, 2 - MINRES.");
322  args.AddOption(&max_lin_iter, "-li", "--lin-iter",
323  "Maximum number of iterations in the linear solve.");
324  args.AddOption(&move_bnd, "-bnd", "--move-boundary", "-fix-bnd",
325  "--fix-boundary",
326  "Enable motion along horizontal and vertical boundaries.");
327  args.AddOption(&combomet, "-cmb", "--combo-met", "-no-cmb", "--no-combo-met",
328  "Combination of metrics.");
329  args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
330  "--no-normalization",
331  "Make all terms in the optimization functional unitless.");
332  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
333  "--no-visualization",
334  "Enable or disable GLVis visualization.");
335  args.AddOption(&verbosity_level, "-vl", "--verbosity-level",
336  "Set the verbosity level - 0, 1, or 2.");
337  args.Parse();
338  if (!args.Good())
339  {
340  args.PrintUsage(cout);
341  return 1;
342  }
343  args.PrintOptions(cout);
344 
345  // 2. Initialize and refine the starting mesh.
346  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
347  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
348  const int dim = mesh->Dimension();
349  cout << "Mesh curvature: ";
350  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
351  else { cout << "(NONE)"; }
352  cout << endl;
353 
354  // 3. Define a finite element space on the mesh. Here we use vector finite
355  // elements which are tensor products of quadratic finite elements. The
356  // number of components in the vector finite element space is specified by
357  // the last parameter of the FiniteElementSpace constructor.
359  if (mesh_poly_deg <= 0)
360  {
361  fec = new QuadraticPosFECollection;
362  mesh_poly_deg = 2;
363  }
364  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
365  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim);
366 
367  // 4. Make the mesh curved based on the above finite element space. This
368  // means that we define the mesh elements through a fespace-based
369  // transformation of the reference element.
370  mesh->SetNodalFESpace(fespace);
371 
372  // 5. Set up an empty right-hand side vector b, which is equivalent to b=0.
373  Vector b(0);
374 
375  // 6. Get the mesh nodes (vertices and other degrees of freedom in the finite
376  // element space) as a finite element grid function in fespace. Note that
377  // changing x automatically changes the shapes of the mesh elements.
378  GridFunction *x = mesh->GetNodes();
379 
380  // 7. Define a vector representing the minimal local mesh size in the mesh
381  // nodes. We index the nodes using the scalar version of the degrees of
382  // freedom in pfespace. Note: this is partition-dependent.
383  //
384  // In addition, compute average mesh size and total volume.
385  Vector h0(fespace->GetNDofs());
386  h0 = infinity();
387  double volume = 0.0;
388  Array<int> dofs;
389  for (int i = 0; i < mesh->GetNE(); i++)
390  {
391  // Get the local scalar element degrees of freedom in dofs.
392  fespace->GetElementDofs(i, dofs);
393  // Adjust the value of h0 in dofs based on the local mesh size.
394  const double hi = mesh->GetElementSize(i);
395  for (int j = 0; j < dofs.Size(); j++)
396  {
397  h0(dofs[j]) = min(h0(dofs[j]), hi);
398  }
399  volume += mesh->GetElementVolume(i);
400  }
401  const double small_phys_size = pow(volume, 1.0 / dim) / 100.0;
402 
403  // 8. Add a random perturbation to the nodes in the interior of the domain.
404  // We define a random grid function of fespace and make sure that it is
405  // zero on the boundary and its values are locally of the order of h0.
406  // The latter is based on the DofToVDof() method which maps the scalar to
407  // the vector degrees of freedom in fespace.
408  GridFunction rdm(fespace);
409  rdm.Randomize();
410  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
411  rdm *= jitter;
412  // Scale the random values to be of order of the local mesh size.
413  for (int i = 0; i < fespace->GetNDofs(); i++)
414  {
415  for (int d = 0; d < dim; d++)
416  {
417  rdm(fespace->DofToVDof(i,d)) *= h0(i);
418  }
419  }
420  Array<int> vdofs;
421  for (int i = 0; i < fespace->GetNBE(); i++)
422  {
423  // Get the vector degrees of freedom in the boundary element.
424  fespace->GetBdrElementVDofs(i, vdofs);
425  // Set the boundary values to zero.
426  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
427  }
428  *x -= rdm;
429  // Set the perturbation of all nodes from the true nodes.
430  x->SetTrueVector();
431  x->SetFromTrueVector();
432 
433  // 9. Save the starting (prior to the optimization) mesh to a file. This
434  // output can be viewed later using GLVis: "glvis -m perturbed.mesh".
435  {
436  ofstream mesh_ofs("perturbed.mesh");
437  mesh->Print(mesh_ofs);
438  }
439 
440  // 10. Store the starting (prior to the optimization) positions.
441  GridFunction x0(fespace);
442  x0 = *x;
443 
444  // 11. Form the integrator that uses the chosen metric and target.
445  double tauval = -0.1;
446  TMOP_QualityMetric *metric = NULL;
447  switch (metric_id)
448  {
449  case 1: metric = new TMOP_Metric_001; break;
450  case 2: metric = new TMOP_Metric_002; break;
451  case 7: metric = new TMOP_Metric_007; break;
452  case 9: metric = new TMOP_Metric_009; break;
453  case 22: metric = new TMOP_Metric_022(tauval); break;
454  case 50: metric = new TMOP_Metric_050; break;
455  case 55: metric = new TMOP_Metric_055; break;
456  case 56: metric = new TMOP_Metric_056; break;
457  case 58: metric = new TMOP_Metric_058; break;
458  case 77: metric = new TMOP_Metric_077; break;
459  case 211: metric = new TMOP_Metric_211; break;
460  case 252: metric = new TMOP_Metric_252(tauval); break;
461  case 301: metric = new TMOP_Metric_301; break;
462  case 302: metric = new TMOP_Metric_302; break;
463  case 303: metric = new TMOP_Metric_303; break;
464  case 315: metric = new TMOP_Metric_315; break;
465  case 316: metric = new TMOP_Metric_316; break;
466  case 321: metric = new TMOP_Metric_321; break;
467  case 352: metric = new TMOP_Metric_352(tauval); break;
468  default: cout << "Unknown metric_id: " << metric_id << endl; return 3;
469  }
471  switch (target_id)
472  {
473  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
474  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
475  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
476  default: cout << "Unknown target_id: " << target_id << endl;
477  delete metric; return 3;
478  }
479  TargetConstructor *target_c = new TargetConstructor(target_t);
480  target_c->SetNodes(x0);
481  TMOP_Integrator *he_nlf_integ = new TMOP_Integrator(metric, target_c);
482 
483  // 12. Setup the quadrature rule for the non-linear form integrator.
484  const IntegrationRule *ir = NULL;
485  const int geom_type = fespace->GetFE(0)->GetGeomType();
486  switch (quad_type)
487  {
488  case 1: ir = &IntRulesLo.Get(geom_type, quad_order); break;
489  case 2: ir = &IntRules.Get(geom_type, quad_order); break;
490  case 3: ir = &IntRulesCU.Get(geom_type, quad_order); break;
491  default: cout << "Unknown quad_type: " << quad_type << endl;
492  delete he_nlf_integ; return 3;
493  }
494  cout << "Quadrature points per cell: " << ir->GetNPoints() << endl;
495  he_nlf_integ->SetIntegrationRule(*ir);
496 
497  if (normalization) { he_nlf_integ->EnableNormalization(x0); }
498 
499  // 13. Limit the node movement.
500  // The limiting distances can be given by a general function of space.
501  GridFunction dist(fespace);
502  dist = 1.0;
503  // The small_phys_size is relevant only with proper normalization.
504  if (normalization) { dist = small_phys_size; }
505  ConstantCoefficient lim_coeff(lim_const);
506  if (lim_const != 0.0) { he_nlf_integ->EnableLimiting(x0, dist, lim_coeff); }
507 
508  // 14. Setup the final NonlinearForm (which defines the integral of interest,
509  // its first and second derivatives). Here we can use a combination of
510  // metrics, i.e., optimize the sum of two integrals, where both are
511  // scaled by used-defined space-dependent weights. Note that there are no
512  // command-line options for the weights and the type of the second
513  // metric; one should update those in the code.
514  NonlinearForm a(fespace);
515  ConstantCoefficient *coeff1 = NULL;
516  TMOP_QualityMetric *metric2 = NULL;
517  TargetConstructor *target_c2 = NULL;
519 
520  if (combomet == 1)
521  {
522  // TODO normalization of combinations.
523  // We will probably drop this example and replace it with adaptivity.
524  if (normalization) { MFEM_ABORT("Not implemented."); }
525 
526  // Weight of the original metric.
527  coeff1 = new ConstantCoefficient(1.0);
528  he_nlf_integ->SetCoefficient(*coeff1);
529  a.AddDomainIntegrator(he_nlf_integ);
530 
531  metric2 = new TMOP_Metric_077;
532  target_c2 = new TargetConstructor(
534  target_c2->SetVolumeScale(0.01);
535  target_c2->SetNodes(x0);
536  TMOP_Integrator *he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
537  he_nlf_integ2->SetIntegrationRule(*ir);
538 
539  // Weight of metric2.
540  he_nlf_integ2->SetCoefficient(coeff2);
541  a.AddDomainIntegrator(he_nlf_integ2);
542  }
543  else { a.AddDomainIntegrator(he_nlf_integ); }
544 
545  const double init_energy = a.GetGridFunctionEnergy(*x);
546 
547  // 15. Visualize the starting mesh and metric values.
548  if (visualization)
549  {
550  char title[] = "Initial metric values";
551  vis_metric(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
552  }
553 
554  // 16. Fix all boundary nodes, or fix only a given component depending on the
555  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
556  // fixed x/y/z components of the node. Attribute 4 corresponds to an
557  // entirely fixed node. Other boundary attributes do not affect the node
558  // movement boundary conditions.
559  if (move_bnd == false)
560  {
561  Array<int> ess_bdr(mesh->bdr_attributes.Max());
562  ess_bdr = 1;
563  a.SetEssentialBC(ess_bdr);
564  }
565  else
566  {
567  const int nd = fespace->GetBE(0)->GetDof();
568  int n = 0;
569  for (int i = 0; i < mesh->GetNBE(); i++)
570  {
571  const int attr = mesh->GetBdrElement(i)->GetAttribute();
572  MFEM_VERIFY(!(dim == 2 && attr == 3),
573  "Boundary attribute 3 must be used only for 3D meshes. "
574  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
575  "components, rest for free nodes), or use -fix-bnd.");
576  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
577  if (attr == 4) { n += nd * dim; }
578  }
579  Array<int> ess_vdofs(n), vdofs;
580  n = 0;
581  for (int i = 0; i < mesh->GetNBE(); i++)
582  {
583  const int attr = mesh->GetBdrElement(i)->GetAttribute();
584  fespace->GetBdrElementVDofs(i, vdofs);
585  if (attr == 1) // Fix x components.
586  {
587  for (int j = 0; j < nd; j++)
588  { ess_vdofs[n++] = vdofs[j]; }
589  }
590  else if (attr == 2) // Fix y components.
591  {
592  for (int j = 0; j < nd; j++)
593  { ess_vdofs[n++] = vdofs[j+nd]; }
594  }
595  else if (attr == 3) // Fix z components.
596  {
597  for (int j = 0; j < nd; j++)
598  { ess_vdofs[n++] = vdofs[j+2*nd]; }
599  }
600  else if (attr == 4) // Fix all components.
601  {
602  for (int j = 0; j < vdofs.Size(); j++)
603  { ess_vdofs[n++] = vdofs[j]; }
604  }
605  }
606  a.SetEssentialVDofs(ess_vdofs);
607  }
608 
609  // 17. As we use the Newton method to solve the resulting nonlinear system,
610  // here we setup the linear solver for the system's Jacobian.
611  Solver *S = NULL;
612  const double linsol_rtol = 1e-12;
613  if (lin_solver == 0)
614  {
615  S = new DSmoother(1, 1.0, max_lin_iter);
616  }
617  else if (lin_solver == 1)
618  {
619  CGSolver *cg = new CGSolver;
620  cg->SetMaxIter(max_lin_iter);
621  cg->SetRelTol(linsol_rtol);
622  cg->SetAbsTol(0.0);
623  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
624  S = cg;
625  }
626  else
627  {
628  MINRESSolver *minres = new MINRESSolver;
629  minres->SetMaxIter(max_lin_iter);
630  minres->SetRelTol(linsol_rtol);
631  minres->SetAbsTol(0.0);
632  minres->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
633  S = minres;
634  }
635 
636  // 18. Compute the minimum det(J) of the starting mesh.
637  tauval = infinity();
638  const int NE = mesh->GetNE();
639  for (int i = 0; i < NE; i++)
640  {
642  for (int j = 0; j < ir->GetNPoints(); j++)
643  {
644  transf->SetIntPoint(&ir->IntPoint(j));
645  tauval = min(tauval, transf->Jacobian().Det());
646  }
647  }
648  cout << "Minimum det(J) of the original mesh is " << tauval << endl;
649 
650  // 19. Finally, perform the nonlinear optimization.
651  NewtonSolver *newton = NULL;
652  if (tauval > 0.0)
653  {
654  tauval = 0.0;
655  newton = new RelaxedNewtonSolver(*ir, fespace);
656  cout << "The RelaxedNewtonSolver is used (as all det(J)>0)." << endl;
657  }
658  else
659  {
660  if ( (dim == 2 && metric_id != 22 && metric_id != 252) ||
661  (dim == 3 && metric_id != 352) )
662  {
663  cout << "The mesh is inverted. Use an untangling metric." << endl;
664  return 3;
665  }
666  tauval -= 0.01 * h0.Min(); // Slightly below minJ0 to avoid div by 0.
667  newton = new DescentNewtonSolver(*ir, fespace);
668  cout << "The DescentNewtonSolver is used (as some det(J)<0)." << endl;
669  }
670  newton->SetPreconditioner(*S);
671  newton->SetMaxIter(newton_iter);
672  newton->SetRelTol(newton_rtol);
673  newton->SetAbsTol(0.0);
674  newton->SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
675  newton->SetOperator(a);
676  newton->Mult(b, x->GetTrueVector());
677  x->SetFromTrueVector();
678  if (newton->GetConverged() == false)
679  {
680  cout << "NewtonIteration: rtol = " << newton_rtol << " not achieved."
681  << endl;
682  }
683  delete newton;
684 
685  // 20. Save the optimized mesh to a file. This output can be viewed later
686  // using GLVis: "glvis -m optimized.mesh".
687  {
688  ofstream mesh_ofs("optimized.mesh");
689  mesh_ofs.precision(14);
690  mesh->Print(mesh_ofs);
691  }
692 
693  // 21. Compute the amount of energy decrease.
694  const double fin_energy = a.GetGridFunctionEnergy(*x);
695  double metric_part = fin_energy;
696  if (lim_const != 0.0)
697  {
698  lim_coeff.constant = 0.0;
699  metric_part = a.GetGridFunctionEnergy(*x);
700  lim_coeff.constant = lim_const;
701  }
702  cout << "Initial strain energy: " << init_energy
703  << " = metrics: " << init_energy
704  << " + limiting term: " << 0.0 << endl;
705  cout << " Final strain energy: " << fin_energy
706  << " = metrics: " << metric_part
707  << " + limiting term: " << fin_energy - metric_part << endl;
708  cout << "The strain energy decreased by: " << setprecision(12)
709  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
710 
711  // 22. Visualize the final mesh and metric values.
712  if (visualization)
713  {
714  char title[] = "Final metric values";
715  vis_metric(mesh_poly_deg, *metric, *target_c, *mesh, title, 600);
716  }
717 
718  // 23. Visualize the mesh displacement.
719  if (visualization)
720  {
721  x0 -= *x;
722  osockstream sock(19916, "localhost");
723  sock << "solution\n";
724  mesh->Print(sock);
725  x0.Save(sock);
726  sock.send();
727  sock << "window_title 'Displacements'\n"
728  << "window_geometry "
729  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
730  << "keys jRmclA" << endl;
731  }
732 
733  // 24. Free the used memory.
734  delete S;
735  delete target_c2;
736  delete metric2;
737  delete coeff1;
738  delete target_c;
739  delete metric;
740  delete fespace;
741  delete fec;
742  delete mesh;
743 
744  return 0;
745 }
746 
747 // Defined with respect to the icf mesh.
748 double weight_fun(const Vector &x)
749 {
750  const double r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
751  const double den = 0.002;
752  double l2 = 0.2 + 0.5*std::tanh((r-0.16)/den) - 0.5*std::tanh((r-0.17)/den)
753  + 0.5*std::tanh((r-0.23)/den) - 0.5*std::tanh((r-0.24)/den);
754  return l2;
755 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:237
int Size() const
Logical size of the array.
Definition: array.hpp:118
Shifted barrier form of metric 56 (area, ideal barrier metric), 2D.
Definition: tmop.hpp:325
Shifted barrier form of 3D metric 16 (volume, ideal barrier metric), 3D.
Definition: tmop.hpp:443
virtual void Print(std::ostream &out=mfem::out) const
Definition: mesh.hpp:1156
Conjugate gradient method.
Definition: solvers.hpp:111
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:344
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:85
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:143
Shape &amp; volume, ideal barrier metric, 3D.
Definition: tmop.hpp:427
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:877
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:139
Subclass constant coefficient.
Definition: coefficient.hpp:67
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
Definition: tmop.cpp:1190
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:679
void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric, const TargetConstructor &tc, const Mesh &mesh, GridFunction &metric_gf)
Interpolates the metric&#39;s values at the nodes of metric_gf.
Definition: tmop.cpp:1259
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:377
double Det() const
Definition: densemat.cpp:472
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:53
aka closed Newton-Cotes
Definition: intrules.hpp:287
Container class for integration rules.
Definition: intrules.hpp:299
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:361
int Size() const
Returns the size of the vector.
Definition: vector.hpp:150
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:589
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:152
Volume metric, 3D.
Definition: tmop.hpp:393
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:289
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:676
MINRES method.
Definition: solvers.hpp:221
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:690
int main(int argc, char *argv[])
Definition: ex1.cpp:58
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:159
double weight_fun(const Vector &x)
Shape &amp; area, ideal barrier metric, 2D.
Definition: tmop.hpp:168
void vis_metric(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
Shifted barrier form of metric 2 (shape, ideal barrier metric), 2D.
Definition: tmop.hpp:200
virtual double GetEnergy(const Vector &x) const
Compute the enery corresponding to the state x.
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:238
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:271
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2627
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:308
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:24
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:240
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:133
int dim
Definition: ex3.cpp:48
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:67
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Adds a limiting term to the integrator (general version).
Definition: tmop.cpp:867
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:383
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Definition: tmop.hpp:663
Area metric, 2D.
Definition: tmop.hpp:235
Volume, ideal barrier metric, 3D.
Definition: tmop.hpp:409
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7591
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:68
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:586
void SetMaxIter(int max_it)
Definition: solvers.hpp:63
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:1439
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:259
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:383
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3621
int GetConverged() const
Definition: solvers.hpp:67
int Dimension() const
Definition: mesh.hpp:713
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: solvers.cpp:1250
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:121
double GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Definition: mesh.cpp:74
void SetAbsTol(double atol)
Definition: solvers.hpp:62
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:179
void SetRelTol(double rtol)
Definition: solvers.hpp:61
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:219
Untangling metric, 2D.
Definition: tmop.hpp:306
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:85
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:311
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:252
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)
Definition: optparser.hpp:76
Shape &amp; area metric, 2D.
Definition: tmop.hpp:184
double GetGridFunctionEnergy(const Vector &x) const
Compute the enery corresponding to the state x.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:336
void SetIntegrationRule(const IntegrationRule &irule)
Prescribe a fixed IntegrationRule to use.
Definition: nonlininteg.hpp:40
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1541
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:42
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
class for C-function coefficient
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:3725
Vector data type.
Definition: vector.hpp:48
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5837
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:88
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:530
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
(DEPRECATED) Specify essential boundary conditions.
Base class for solvers.
Definition: operator.hpp:279
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i&#39;th boundary element.
Definition: fespace.cpp:1781
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:178
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:526
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:368
double GetElementVolume(int i)
Definition: mesh.cpp:101
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:748
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:345
Metric without a type, 2D.
Definition: tmop.hpp:76
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1239
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:134
bool Good() const
Definition: optparser.hpp:122
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:607