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