MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
pmesh-optimizer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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 - Parallel Version
14 // ---------------------------------------------------------------------
15 //
16 // This miniapp performs mesh optimization using the Target-Matrix Optimization
17 // Paradigm (TMOP) by P.Knupp et al., and a global variational minimization
18 // approach. It minimizes the quantity sum_T int_T mu(J(x)), where T are the
19 // target (ideal) elements, J is the Jacobian of the transformation from the
20 // target to the physical element, and mu is the mesh quality metric. This
21 // metric can measure shape, size or alignment of the region around each
22 // quadrature point. The combination of targets & quality metrics is used to
23 // optimize the physical node positions, i.e., they must be as close as possible
24 // to the shape / size / alignment of their targets. This code also demonstrates
25 // a possible use of nonlinear operators (the class TMOP_QualityMetric, defining
26 // mu(J), and the class TMOP_Integrator, defining int mu(J)), as well as their
27 // coupling to Newton methods for solving minimization problems. Note that the
28 // utilized Newton methods are oriented towards avoiding invalid meshes with
29 // negative Jacobian determinants. Each Newton step requires the inversion of a
30 // Jacobian matrix, which is done through an inner linear solver.
31 //
32 // Compile with: make pmesh-optimizer
33 //
34 // Sample runs:
35 // Adapted analytic Hessian:
36 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
37 // Adapted discrete size:
38 // mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 5 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
39 //
40 // Blade shape:
41 // mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -rs 0 -mid 2 -tid 1 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
42 // Blade limited shape:
43 // mpirun -np 4 pmesh-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
44 // ICF shape and equal size:
45 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
46 // ICF shape and initial size:
47 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 9 -tid 3 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8
48 // ICF shape:
49 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8
50 // ICF limited shape:
51 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 100 -ls 2 -li 100 -bnd -qt 1 -qo 8 -lc 10
52 // ICF combo shape + size (rings, slow convergence):
53 // mpirun -np 4 pmesh-optimizer -o 3 -rs 0 -mid 1 -tid 1 -ni 1000 -ls 2 -li 100 -bnd -qt 1 -qo 8 -cmb
54 // 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
55 // * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -rs 0 -mid 303 -tid 1 -ni 20 -ls 2 -li 500 -fix-bnd
56 // 2D non-conforming shape and equal size:
57 // mpirun -np 4 pmesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -ls 2 -li 100 -bnd -qt 1 -qo 8
58 
59 #include "mfem.hpp"
60 #include <iostream>
61 #include <fstream>
62 
63 using namespace mfem;
64 using namespace std;
65 
66 double weight_fun(const Vector &x);
67 
68 double ind_values(const Vector &x)
69 {
70  const int opt = 6;
71  const double small = 0.001, big = 0.01;
72 
73  // Sine wave.
74  if (opt==1)
75  {
76  const double X = x(0), Y = x(1);
77  const double ind = std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) + 1) -
78  std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) - 1);
79 
80  return ind * small + (1.0 - ind) * big;
81  }
82 
83  if (opt==2)
84  {
85  // Circle in the middle.
86  double val = 0.;
87  const double xc = x(0) - 0.5, yc = x(1) - 0.5;
88  const double r = sqrt(xc*xc + yc*yc);
89  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
90  val = 0.5*(std::tanh(sf*(r-r1)) - std::tanh(sf*(r-r2)));
91  if (val > 1.) {val = 1;}
92 
93  return val * small + (1.0 - val) * big;
94  }
95 
96  if (opt == 3)
97  {
98  // cross
99  const double X = x(0), Y = x(1);
100  const double r1 = 0.45, r2 = 0.55;
101  const double sf = 40.0;
102 
103  double val = 0.5 * ( std::tanh(sf*(X-r1)) - std::tanh(sf*(X-r2)) +
104  std::tanh(sf*(Y-r1)) - std::tanh(sf*(Y-r2)) );
105  if (val > 1.) { val = 1.0; }
106 
107  return val * small + (1.0 - val) * big;
108  }
109 
110  if (opt==4)
111  {
112  // Multiple circles
113  double r1,r2,val,rval;
114  double sf = 10;
115  val = 0.;
116  // circle 1
117  r1= 0.25; r2 = 0.25; rval = 0.1;
118  double xc = x(0) - r1, yc = x(1) - r2;
119  double r = sqrt(xc*xc+yc*yc);
120  val = 0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
121  (r-rval)));// std::exp(val1);
122  // circle 2
123  r1= 0.75; r2 = 0.75;
124  xc = x(0) - r1, yc = x(1) - r2;
125  r = sqrt(xc*xc+yc*yc);
126  val += (0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
127  (r-rval))));// std::exp(val1);
128  // circle 3
129  r1= 0.75; r2 = 0.25;
130  xc = x(0) - r1, yc = x(1) - r2;
131  r = sqrt(xc*xc+yc*yc);
132  val += 0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*
133  (r-rval)));// std::exp(val1);
134  // circle 4
135  r1= 0.25; r2 = 0.75;
136  xc = x(0) - r1, yc = x(1) - r2;
137  r = sqrt(xc*xc+yc*yc);
138  val += 0.5*(1+std::tanh(sf*(r+rval))) - 0.5*(1+std::tanh(sf*(r-rval)));
139  if (val > 1.0) {val = 1.;}
140  if (val < 0.0) {val = 0.;}
141 
142  return val * small + (1.0 - val) * big;
143  }
144 
145  if (opt==5)
146  {
147  // cross
148  double val = 0.;
149  double X = x(0)-0.5, Y = x(1)-0.5;
150  double rval = std::sqrt(X*X + Y*Y);
151  double thval = 60.*M_PI/180.;
152  double Xmod,Ymod;
153  Xmod = X*std::cos(thval) + Y*std::sin(thval);
154  Ymod= -X*std::sin(thval) + Y*std::cos(thval);
155  X = Xmod+0.5; Y = Ymod+0.5;
156  double r1 = 0.45; double r2 = 0.55; double sf=30.0;
157  val = ( 0.5*(1+std::tanh(sf*(X-r1))) - 0.5*(1+std::tanh(sf*(X-r2)))
158  + 0.5*(1+std::tanh(sf*(Y-r1))) - 0.5*(1+std::tanh(sf*(Y-r2))) );
159  if (rval > 0.4) {val = 0.;}
160  if (val > 1.0) {val = 1.;}
161  if (val < 0.0) {val = 0.;}
162 
163  return val * small + (1.0 - val) * big;
164  }
165 
166  if (opt==6)
167  {
168  double val = 0.;
169  const double xc = x(0) - 0.0, yc = x(1) - 0.5;
170  const double r = sqrt(xc*xc + yc*yc);
171  double r1 = 0.45; double r2 = 0.55; double sf=30.0;
172  val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
173  if (val > 1.) {val = 1;}
174  if (val < 0.) {val = 0;}
175 
176  return val * small + (1.0 - val) * big;
177  }
178 
179  return 0.0;
180 }
181 
182 class HessianCoefficient : public MatrixCoefficient
183 {
184 private:
185  int type;
186 
187 public:
188  HessianCoefficient(int dim, int type_)
189  : MatrixCoefficient(dim), type(type_) { }
190 
191  virtual void Eval(DenseMatrix &K, ElementTransformation &T,
192  const IntegrationPoint &ip)
193  {
194  Vector pos(3);
195  T.Transform(ip, pos);
196 
197  if (type == 0)
198  {
199  K(0, 0) = 1.0 + 3.0 * std::sin(M_PI*pos(0));
200  K(0, 1) = 0.0;
201  K(1, 0) = 0.0;
202  K(1, 1) = 1.0;
203  }
204  else
205  {
206  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
207  const double r = sqrt(xc*xc + yc*yc);
208  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
209  const double eps = 0.5;
210 
211  const double tan1 = std::tanh(sf*(r-r1)),
212  tan2 = std::tanh(sf*(r-r2));
213 
214  K(0, 0) = eps + 1.0 * (tan1 - tan2);
215  K(0, 1) = 0.0;
216  K(1, 0) = 0.0;
217  K(1, 1) = 1.0;
218  }
219  }
220 };
221 
222 // Additional IntegrationRules that can be used with the --quad-type option.
225 
226 int main (int argc, char *argv[])
227 {
228  // 0. Initialize MPI.
229  int num_procs, myid;
230  MPI_Init(&argc, &argv);
231  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
232  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
233 
234  // 1. Set the method's default parameters.
235  const char *mesh_file = "icf.mesh";
236  int mesh_poly_deg = 1;
237  int rs_levels = 0;
238  int rp_levels = 0;
239  double jitter = 0.0;
240  int metric_id = 1;
241  int target_id = 1;
242  double lim_const = 0.0;
243  int quad_type = 1;
244  int quad_order = 8;
245  int newton_iter = 10;
246  double newton_rtol = 1e-12;
247  int lin_solver = 2;
248  int max_lin_iter = 100;
249  bool move_bnd = true;
250  bool combomet = 0;
251  bool normalization = false;
252  bool visualization = true;
253  int verbosity_level = 0;
254 
255  // 2. Parse command-line options.
256  OptionsParser args(argc, argv);
257  args.AddOption(&mesh_file, "-m", "--mesh",
258  "Mesh file to use.");
259  args.AddOption(&mesh_poly_deg, "-o", "--order",
260  "Polynomial degree of mesh finite element space.");
261  args.AddOption(&rs_levels, "-rs", "--refine-serial",
262  "Number of times to refine the mesh uniformly in serial.");
263  args.AddOption(&rp_levels, "-rp", "--refine-parallel",
264  "Number of times to refine the mesh uniformly in parallel.");
265  args.AddOption(&jitter, "-ji", "--jitter",
266  "Random perturbation scaling factor.");
267  args.AddOption(&metric_id, "-mid", "--metric-id",
268  "Mesh optimization metric:\n\t"
269  "1 : |T|^2 -- 2D shape\n\t"
270  "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
271  "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
272  "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
273  "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
274  "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
275  "55 : (tau-1)^2 -- 2D size\n\t"
276  "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
277  "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
278  "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
279  "211: (tau-1)^2-tau+sqrt(tau^2) -- 2D untangling\n\t"
280  "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
281  "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
282  "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
283  "303: (|T|^2)/3*tau^(2/3)-1 -- 3D shape\n\t"
284  "315: (tau-1)^2 -- 3D size\n\t"
285  "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D size\n\t"
286  "321: |T-T^-t|^2 -- 3D shape+size\n\t"
287  "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling");
288  args.AddOption(&target_id, "-tid", "--target-id",
289  "Target (ideal element) type:\n\t"
290  "1: Ideal shape, unit size\n\t"
291  "2: Ideal shape, equal size\n\t"
292  "3: Ideal shape, initial size\n\t"
293  "4: Given full analytic Jacobian (in physical space)\n\t"
294  "5: Ideal shape, given size (in physical space)");
295  args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
296  args.AddOption(&quad_type, "-qt", "--quad-type",
297  "Quadrature rule type:\n\t"
298  "1: Gauss-Lobatto\n\t"
299  "2: Gauss-Legendre\n\t"
300  "3: Closed uniform points");
301  args.AddOption(&quad_order, "-qo", "--quad_order",
302  "Order of the quadrature rule.");
303  args.AddOption(&newton_iter, "-ni", "--newton-iters",
304  "Maximum number of Newton iterations.");
305  args.AddOption(&newton_rtol, "-rtol", "--newton-rel-tolerance",
306  "Relative tolerance for the Newton solver.");
307  args.AddOption(&lin_solver, "-ls", "--lin-solver",
308  "Linear solver: 0 - l1-Jacobi, 1 - CG, 2 - MINRES.");
309  args.AddOption(&max_lin_iter, "-li", "--lin-iter",
310  "Maximum number of iterations in the linear solve.");
311  args.AddOption(&move_bnd, "-bnd", "--move-boundary", "-fix-bnd",
312  "--fix-boundary",
313  "Enable motion along horizontal and vertical boundaries.");
314  args.AddOption(&combomet, "-cmb", "--combo-met", "-no-cmb", "--no-combo-met",
315  "Combination of metrics.");
316  args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
317  "--no-normalization",
318  "Make all terms in the optimization functional unitless.");
319  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
320  "--no-visualization",
321  "Enable or disable GLVis visualization.");
322  args.AddOption(&verbosity_level, "-vl", "--verbosity-level",
323  "Set the verbosity level - 0, 1, or 2.");
324  args.Parse();
325  if (!args.Good())
326  {
327  if (myid == 0) { args.PrintUsage(cout); }
328  return 1;
329  }
330  if (myid == 0) { args.PrintOptions(cout); }
331 
332  // 3. Initialize and refine the starting mesh.
333  Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
334  for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
335  const int dim = mesh->Dimension();
336  if (myid == 0)
337  {
338  cout << "Mesh curvature: ";
339  if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
340  else { cout << "(NONE)"; }
341  cout << endl;
342  }
343  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
344 
345  delete mesh;
346  for (int lev = 0; lev < rp_levels; lev++) { pmesh->UniformRefinement(); }
347 
348  // 4. Define a finite element space on the mesh. Here we use vector finite
349  // elements which are tensor products of quadratic finite elements. The
350  // number of components in the vector finite element space is specified by
351  // the last parameter of the FiniteElementSpace constructor.
353  if (mesh_poly_deg <= 0)
354  {
355  fec = new QuadraticPosFECollection;
356  mesh_poly_deg = 2;
357  }
358  else { fec = new H1_FECollection(mesh_poly_deg, dim); }
359  ParFiniteElementSpace *pfespace = new ParFiniteElementSpace(pmesh, fec, dim);
360 
361  // 5. Make the mesh curved based on the above finite element space. This
362  // means that we define the mesh elements through a fespace-based
363  // transformation of the reference element.
364  pmesh->SetNodalFESpace(pfespace);
365 
366  // 6. Set up an empty right-hand side vector b, which is equivalent to b=0.
367  Vector b(0);
368 
369  // 7. Get the mesh nodes (vertices and other degrees of freedom in the finite
370  // element space) as a finite element grid function in fespace. Note that
371  // changing x automatically changes the shapes of the mesh elements.
372  ParGridFunction x(pfespace);
373  pmesh->SetNodalGridFunction(&x);
374 
375  // 8. Define a vector representing the minimal local mesh size in the mesh
376  // nodes. We index the nodes using the scalar version of the degrees of
377  // freedom in pfespace. Note: this is partition-dependent.
378  //
379  // In addition, compute average mesh size and total volume.
380  Vector h0(pfespace->GetNDofs());
381  h0 = infinity();
382  double vol_loc = 0.0;
383  Array<int> dofs;
384  for (int i = 0; i < pmesh->GetNE(); i++)
385  {
386  // Get the local scalar element degrees of freedom in dofs.
387  pfespace->GetElementDofs(i, dofs);
388  // Adjust the value of h0 in dofs based on the local mesh size.
389  const double hi = pmesh->GetElementSize(i);
390  for (int j = 0; j < dofs.Size(); j++)
391  {
392  h0(dofs[j]) = min(h0(dofs[j]), hi);
393  }
394  vol_loc += pmesh->GetElementVolume(i);
395  }
396  double volume;
397  MPI_Allreduce(&vol_loc, &volume, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
398  const double small_phys_size = pow(volume, 1.0 / dim) / 100.0;
399 
400  // 9. Add a random perturbation to the nodes in the interior of the domain.
401  // We define a random grid function of fespace and make sure that it is
402  // zero on the boundary and its values are locally of the order of h0.
403  // The latter is based on the DofToVDof() method which maps the scalar to
404  // the vector degrees of freedom in fespace.
405  ParGridFunction rdm(pfespace);
406  rdm.Randomize();
407  rdm -= 0.25; // Shift to random values in [-0.5,0.5].
408  rdm *= jitter;
409  // Scale the random values to be of order of the local mesh size.
410  for (int i = 0; i < pfespace->GetNDofs(); i++)
411  {
412  for (int d = 0; d < dim; d++)
413  {
414  rdm(pfespace->DofToVDof(i,d)) *= h0(i);
415  }
416  }
417  Array<int> vdofs;
418  for (int i = 0; i < pfespace->GetNBE(); i++)
419  {
420  // Get the vector degrees of freedom in the boundary element.
421  pfespace->GetBdrElementVDofs(i, vdofs);
422  // Set the boundary values to zero.
423  for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
424  }
425  x -= rdm;
426  // Set the perturbation of all nodes from the true nodes.
427  x.SetTrueVector();
428  x.SetFromTrueVector();
429 
430  // 10. Save the starting (prior to the optimization) mesh to a file. This
431  // output can be viewed later using GLVis: "glvis -m perturbed -np
432  // num_mpi_tasks".
433  {
434  ostringstream mesh_name;
435  mesh_name << "perturbed." << setfill('0') << setw(6) << myid;
436  ofstream mesh_ofs(mesh_name.str().c_str());
437  mesh_ofs.precision(8);
438  pmesh->Print(mesh_ofs);
439  }
440 
441  // 11. Store the starting (prior to the optimization) positions.
442  ParGridFunction x0(pfespace);
443  x0 = x;
444 
445  // 12. Form the integrator that uses the chosen metric and target.
446  double tauval = -0.1;
447  TMOP_QualityMetric *metric = NULL;
448  switch (metric_id)
449  {
450  case 1: metric = new TMOP_Metric_001; break;
451  case 2: metric = new TMOP_Metric_002; break;
452  case 7: metric = new TMOP_Metric_007; break;
453  case 9: metric = new TMOP_Metric_009; break;
454  case 22: metric = new TMOP_Metric_022(tauval); break;
455  case 50: metric = new TMOP_Metric_050; break;
456  case 55: metric = new TMOP_Metric_055; break;
457  case 56: metric = new TMOP_Metric_056; break;
458  case 58: metric = new TMOP_Metric_058; break;
459  case 77: metric = new TMOP_Metric_077; break;
460  case 211: metric = new TMOP_Metric_211; break;
461  case 252: metric = new TMOP_Metric_252(tauval); break;
462  case 301: metric = new TMOP_Metric_301; break;
463  case 302: metric = new TMOP_Metric_302; break;
464  case 303: metric = new TMOP_Metric_303; break;
465  case 315: metric = new TMOP_Metric_315; break;
466  case 316: metric = new TMOP_Metric_316; break;
467  case 321: metric = new TMOP_Metric_321; break;
468  case 352: metric = new TMOP_Metric_352(tauval); break;
469  default:
470  if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
471  return 3;
472  }
474  TargetConstructor *target_c = NULL;
475  HessianCoefficient *adapt_coeff = NULL;
476  H1_FECollection ind_fec(3, dim);
477  ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
478  ParGridFunction size;
479  switch (target_id)
480  {
481  case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
482  case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
483  case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
484  case 4:
485  {
487  AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
488  adapt_coeff = new HessianCoefficient(dim, 1);
489  tc->SetAnalyticTargetSpec(NULL, NULL, adapt_coeff);
490  target_c = tc;
491  break;
492  }
493  case 5:
494  {
496  DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
497  size.SetSpace(&ind_fes);
498  FunctionCoefficient ind_coeff(ind_values);
499  size.ProjectCoefficient(ind_coeff);
500  tc->SetParDiscreteTargetSpec(size);
501  target_c = tc;
502  break;
503  }
504  default:
505  if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
506  return 3;
507  }
508 
509  if (target_c == NULL)
510  {
511  target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
512  }
513  target_c->SetNodes(x0);
514  TMOP_Integrator *he_nlf_integ= new TMOP_Integrator(metric, target_c);
515 
516  // 13. Setup the quadrature rule for the non-linear form integrator.
517  const IntegrationRule *ir = NULL;
518  const int geom_type = pfespace->GetFE(0)->GetGeomType();
519  switch (quad_type)
520  {
521  case 1: ir = &IntRulesLo.Get(geom_type, quad_order); break;
522  case 2: ir = &IntRules.Get(geom_type, quad_order); break;
523  case 3: ir = &IntRulesCU.Get(geom_type, quad_order); break;
524  default:
525  if (myid == 0) { cout << "Unknown quad_type: " << quad_type << endl; }
526  return 3;
527  }
528  if (myid == 0)
529  { cout << "Quadrature points per cell: " << ir->GetNPoints() << endl; }
530  he_nlf_integ->SetIntegrationRule(*ir);
531 
532  if (normalization) { he_nlf_integ->ParEnableNormalization(x0); }
533 
534  // 14. Limit the node movement.
535  // The limiting distances can be given by a general function of space.
536  ParGridFunction dist(pfespace);
537  dist = 1.0;
538  // The small_phys_size is relevant only with proper normalization.
539  if (normalization) { dist = small_phys_size; }
540  ConstantCoefficient lim_coeff(lim_const);
541  if (lim_const != 0.0) { he_nlf_integ->EnableLimiting(x0, dist, lim_coeff); }
542 
543  // 15. Setup the final NonlinearForm (which defines the integral of interest,
544  // its first and second derivatives). Here we can use a combination of
545  // metrics, i.e., optimize the sum of two integrals, where both are
546  // scaled by used-defined space-dependent weights. Note that there are
547  // no command-line options for the weights and the type of the second
548  // metric; one should update those in the code.
549  ParNonlinearForm a(pfespace);
550  ConstantCoefficient *coeff1 = NULL;
551  TMOP_QualityMetric *metric2 = NULL;
552  TargetConstructor *target_c2 = NULL;
554 
555  if (combomet == 1)
556  {
557  // TODO normalization of combinations.
558  // We will probably drop this example and replace it with adaptivity.
559  if (normalization) { MFEM_ABORT("Not implemented."); }
560 
561  // First metric.
562  coeff1 = new ConstantCoefficient(1.0);
563  he_nlf_integ->SetCoefficient(*coeff1);
564  a.AddDomainIntegrator(he_nlf_integ);
565 
566  // Second metric.
567  metric2 = new TMOP_Metric_077;
568  target_c2 = new TargetConstructor(
570  target_c2->SetVolumeScale(0.01);
571  target_c2->SetNodes(x0);
572  TMOP_Integrator *he_nlf_integ2;
573  he_nlf_integ2 = new TMOP_Integrator(metric2, target_c2);
574  he_nlf_integ2->SetIntegrationRule(*ir);
575 
576  // Weight of metric2.
577  he_nlf_integ2->SetCoefficient(coeff2);
578  a.AddDomainIntegrator(he_nlf_integ2);
579  }
580  else { a.AddDomainIntegrator(he_nlf_integ); }
581 
582  const double init_energy = a.GetParGridFunctionEnergy(x);
583 
584  // 16. Visualize the starting mesh and metric values.
585  if (visualization)
586  {
587  char title[] = "Initial metric values";
588  vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
589  }
590 
591  // 17. Fix all boundary nodes, or fix only a given component depending on the
592  // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
593  // fixed x/y/z components of the node. Attribute 4 corresponds to an
594  // entirely fixed node. Other boundary attributes do not affect the node
595  // movement boundary conditions.
596  if (move_bnd == false)
597  {
598  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
599  ess_bdr = 1;
600  a.SetEssentialBC(ess_bdr);
601  }
602  else
603  {
604  const int nd = pfespace->GetBE(0)->GetDof();
605  int n = 0;
606  for (int i = 0; i < pmesh->GetNBE(); i++)
607  {
608  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
609  MFEM_VERIFY(!(dim == 2 && attr == 3),
610  "Boundary attribute 3 must be used only for 3D meshes. "
611  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
612  "components, rest for free nodes), or use -fix-bnd.");
613  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
614  if (attr == 4) { n += nd * dim; }
615  }
616  Array<int> ess_vdofs(n), vdofs;
617  n = 0;
618  for (int i = 0; i < pmesh->GetNBE(); i++)
619  {
620  const int attr = pmesh->GetBdrElement(i)->GetAttribute();
621  pfespace->GetBdrElementVDofs(i, vdofs);
622  if (attr == 1) // Fix x components.
623  {
624  for (int j = 0; j < nd; j++)
625  { ess_vdofs[n++] = vdofs[j]; }
626  }
627  else if (attr == 2) // Fix y components.
628  {
629  for (int j = 0; j < nd; j++)
630  { ess_vdofs[n++] = vdofs[j+nd]; }
631  }
632  else if (attr == 3) // Fix z components.
633  {
634  for (int j = 0; j < nd; j++)
635  { ess_vdofs[n++] = vdofs[j+2*nd]; }
636  }
637  else if (attr == 4) // Fix all components.
638  {
639  for (int j = 0; j < vdofs.Size(); j++)
640  { ess_vdofs[n++] = vdofs[j]; }
641  }
642  }
643  a.SetEssentialVDofs(ess_vdofs);
644  }
645 
646  // 18. As we use the Newton method to solve the resulting nonlinear system,
647  // here we setup the linear solver for the system's Jacobian.
648  Solver *S = NULL;
649  const double linsol_rtol = 1e-12;
650  if (lin_solver == 0)
651  {
652  S = new DSmoother(1, 1.0, max_lin_iter);
653  }
654  else if (lin_solver == 1)
655  {
656  CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
657  cg->SetMaxIter(max_lin_iter);
658  cg->SetRelTol(linsol_rtol);
659  cg->SetAbsTol(0.0);
660  cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
661  S = cg;
662  }
663  else
664  {
665  MINRESSolver *minres = new MINRESSolver(MPI_COMM_WORLD);
666  minres->SetMaxIter(max_lin_iter);
667  minres->SetRelTol(linsol_rtol);
668  minres->SetAbsTol(0.0);
669  minres->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
670  S = minres;
671  }
672 
673  // 19. Compute the minimum det(J) of the starting mesh.
674  tauval = infinity();
675  const int NE = pmesh->GetNE();
676  for (int i = 0; i < NE; i++)
677  {
679  for (int j = 0; j < ir->GetNPoints(); j++)
680  {
681  transf->SetIntPoint(&ir->IntPoint(j));
682  tauval = min(tauval, transf->Jacobian().Det());
683  }
684  }
685  double minJ0;
686  MPI_Allreduce(&tauval, &minJ0, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
687  tauval = minJ0;
688  if (myid == 0)
689  { cout << "Minimum det(J) of the original mesh is " << tauval << endl; }
690 
691  // 20. Finally, perform the nonlinear optimization.
692  NewtonSolver *newton = NULL;
693  if (tauval > 0.0)
694  {
695  tauval = 0.0;
696  TMOPNewtonSolver *tns = new TMOPNewtonSolver(pfespace->GetComm(), *ir);
697  if (target_id == 5)
698  {
699  tns->SetDiscreteAdaptTC(dynamic_cast<DiscreteAdaptTC *>(target_c));
700  }
701  newton = tns;
702  if (myid == 0)
703  { cout << "TMOPNewtonSolver is used (as all det(J) > 0)." << endl; }
704  }
705  else
706  {
707  if ( (dim == 2 && metric_id != 22 && metric_id != 252) ||
708  (dim == 3 && metric_id != 352) )
709  {
710  if (myid == 0)
711  { cout << "The mesh is inverted. Use an untangling metric.\n"; }
712  return 3;
713  }
714  double h0min = h0.Min(), h0min_all;
715  MPI_Allreduce(&h0min, &h0min_all, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
716  tauval -= 0.01 * h0min_all; // Slightly below minJ0 to avoid div by 0.
717  newton = new TMOPDescentNewtonSolver(pfespace->GetComm(), *ir);
718  if (myid == 0)
719  { cout << "TMOPDescentNewtonSolver is used (as some det(J) < 0).\n"; }
720  }
721  newton->SetPreconditioner(*S);
722  newton->SetMaxIter(newton_iter);
723  newton->SetRelTol(newton_rtol);
724  newton->SetAbsTol(0.0);
725  newton->SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
726  newton->SetOperator(a);
727  newton->Mult(b, x.GetTrueVector());
728  x.SetFromTrueVector();
729  if (myid == 0 && newton->GetConverged() == false)
730  {
731  cout << "NewtonIteration: rtol = " << newton_rtol << " not achieved."
732  << endl;
733  }
734  delete newton;
735 
736  // 21. Save the optimized mesh to a file. This output can be viewed later
737  // using GLVis: "glvis -m optimized -np num_mpi_tasks".
738  {
739  ostringstream mesh_name;
740  mesh_name << "optimized." << setfill('0') << setw(6) << myid;
741  ofstream mesh_ofs(mesh_name.str().c_str());
742  mesh_ofs.precision(8);
743  pmesh->Print(mesh_ofs);
744  }
745 
746  // 22. Compute the amount of energy decrease.
747  const double fin_energy = a.GetParGridFunctionEnergy(x);
748  double metric_part = fin_energy;
749  if (lim_const != 0.0)
750  {
751  lim_coeff.constant = 0.0;
752  metric_part = a.GetParGridFunctionEnergy(x);
753  lim_coeff.constant = lim_const;
754  }
755  if (myid == 0)
756  {
757  cout << "Initial strain energy: " << init_energy
758  << " = metrics: " << init_energy
759  << " + limiting term: " << 0.0 << endl;
760  cout << " Final strain energy: " << fin_energy
761  << " = metrics: " << metric_part
762  << " + limiting term: " << fin_energy - metric_part << endl;
763  cout << "The strain energy decreased by: " << setprecision(12)
764  << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
765  }
766 
767  // 23. Visualize the final mesh and metric values.
768  if (visualization)
769  {
770  char title[] = "Final metric values";
771  vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
772  }
773 
774  // 23. Visualize the mesh displacement.
775  if (visualization)
776  {
777  x0 -= x;
778  socketstream sock;
779  if (myid == 0)
780  {
781  sock.open("localhost", 19916);
782  sock << "solution\n";
783  }
784  pmesh->PrintAsOne(sock);
785  x0.SaveAsOne(sock);
786  if (myid == 0)
787  {
788  sock << "window_title 'Displacements'\n"
789  << "window_geometry "
790  << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
791  << "keys jRmclA" << endl;
792  }
793  }
794 
795  // 24. Free the used memory.
796  delete S;
797  delete target_c2;
798  delete metric2;
799  delete coeff1;
800  delete target_c;
801  delete adapt_coeff;
802  delete metric;
803  delete pfespace;
804  delete fec;
805  delete pmesh;
806 
807  MPI_Finalize();
808  return 0;
809 }
810 
811 // Defined with respect to the icf mesh.
812 double weight_fun(const Vector &x)
813 {
814  const double r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
815  const double den = 0.002;
816  double l2 = 0.2 + 0.5 * (std::tanh((r-0.16)/den) - std::tanh((r-0.17)/den)
817  + std::tanh((r-0.23)/den) - std::tanh((r-0.24)/den));
818  return l2;
819 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
Definition: tmop_tools.cpp:473
virtual void SetParDiscreteTargetSpec(ParGridFunction &tspec)
Definition: tmop.cpp:925
int Size() const
Logical size of the array.
Definition: array.hpp:124
Shifted barrier form of metric 56 (area, ideal barrier metric), 2D.
Definition: tmop.hpp:324
Shifted barrier form of 3D metric 16 (volume, ideal barrier metric), 3D.
Definition: tmop.hpp:442
Conjugate gradient method.
Definition: solvers.hpp:152
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:378
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:143
Shape &amp; volume, ideal barrier metric, 3D.
Definition: tmop.hpp:426
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:890
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:139
Subclass constant coefficient.
Definition: coefficient.hpp:67
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:696
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:376
double Det() const
Definition: densemat.cpp:449
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
Allows negative Jacobians. Used for untangling.
Definition: tmop_tools.hpp:106
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:56
void SaveAsOne(std::ostream &out=mfem::out)
Merge the local grid functions.
Definition: pgridfunc.cpp:501
Container class for integration rules.
Definition: intrules.hpp:307
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:360
void SetVolumeScale(double vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:635
Parallel non-linear operator on the true dofs.
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:151
Volume metric, 3D.
Definition: tmop.hpp:392
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:288
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
Definition: pgridfunc.cpp:87
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:301
MINRES method.
Definition: solvers.hpp:262
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:692
int main(int argc, char *argv[])
Definition: ex1.cpp:62
double weight_fun(const Vector &x)
Shape &amp; area, ideal barrier metric, 2D.
Definition: tmop.hpp:167
Shifted barrier form of metric 2 (shape, ideal barrier metric), 2D.
Definition: tmop.hpp:199
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:270
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:311
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:133
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
MPI_Comm GetComm() const
Definition: pfespace.hpp:239
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:1038
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:417
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Definition: tmop.hpp:784
Area metric, 2D.
Definition: tmop.hpp:234
Volume, ideal barrier metric, 3D.
Definition: tmop.hpp:408
double b
Definition: lissajous.cpp:42
void PrintAsOne(std::ostream &out=mfem::out)
Definition: pmesh.cpp:4231
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:7982
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:71
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:632
void SetMaxIter(int max_it)
Definition: solvers.hpp:65
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
Newton&#39;s method for solving F(x)=b for a given operator F.
Definition: solvers.hpp:300
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4102
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:394
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
double ind_values(const Vector &x)
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:3885
int GetConverged() const
Definition: solvers.hpp:69
int Dimension() const
Definition: mesh.hpp:744
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:1339
double GetParGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
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:76
void SetAbsTol(double atol)
Definition: solvers.hpp:64
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:189
void SetRelTol(double rtol)
Definition: solvers.hpp:63
Shape, ideal barrier metric, 2D.
Definition: tmop.hpp:218
Untangling metric, 2D.
Definition: tmop.hpp:305
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:314
Area, ideal barrier metric, 2D.
Definition: tmop.hpp:251
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
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, MatrixCoefficient *mspec)
Definition: tmop.cpp:882
Shape &amp; area metric, 2D.
Definition: tmop.hpp:183
double a
Definition: lissajous.cpp:41
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: pfespace.cpp:460
Class for integration point with weight.
Definition: intrules.hpp:25
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:338
void SetIntegrationRule(const IntegrationRule &irule)
Prescribe a fixed IntegrationRule to use.
Definition: nonlininteg.hpp:41
int dim
Definition: ex24.cpp:43
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
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.
int open(const char hostname[], int port)
class for C-function coefficient
Vector data type.
Definition: vector.hpp:48
virtual void Transform(const IntegrationPoint &, Vector &)=0
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6203
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:574
void SetDiscreteAdaptTC(DiscreteAdaptTC *tc)
Definition: tmop_tools.hpp:98
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
(DEPRECATED) Specify essential boundary conditions.
Base class for solvers.
Definition: operator.hpp:500
Class for parallel grid function.
Definition: pgridfunc.hpp:32
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i&#39;th boundary element.
Definition: fespace.cpp:1882
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:178
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:3913
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:570
Class for parallel meshes.
Definition: pmesh.hpp:32
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:376
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:1374
double GetElementVolume(int i)
Definition: mesh.cpp:103
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:779
aka closed Newton-Cotes
Definition: intrules.hpp:295
Shape, ideal barrier metric, 3D.
Definition: tmop.hpp:344
Metric without a type, 2D.
Definition: tmop.hpp:75
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:1328
bool Good() const
Definition: optparser.hpp:125
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:728