MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pmesh-optimizer.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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 shape:
36// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -bnd -qt 1 -qo 8
37// Adapted analytic size+orientation:
38// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 200 -bnd -qt 1 -qo 8
39// Adapted analytic shape+orientation:
40// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
41//
42// Adapted analytic shape and/or size with hr-adaptivity:
43// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 9 -ni 50 -li 20 -hmid 55 -mid 7 -hr
44// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 10 -ni 50 -li 20 -hmid 55 -mid 7 -hr
45// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -tid 11 -ni 50 -li 20 -hmid 58 -mid 7 -hr
46//
47// Adapted discrete size:
48// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 94 -tid 5 -ni 50 -qo 4 -nor
49// (requires GSLIB):
50// * mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 80 -tid 5 -ni 50 -qo 4 -nor -mno 1 -ae 1
51// Adapted discrete size NC mesh;
52// mpirun -np 4 pmesh-optimizer -m amr-quad-q2.mesh -o 2 -rs 2 -mid 94 -tid 5 -ni 50 -qo 4 -nor
53// Adapted discrete size 3D with PA:
54// mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 321 -tid 5 -ls 3 -nor -pa
55// Adapted discrete size 3D with PA on device (requires CUDA):
56// * mpirun -n 4 pmesh-optimizer -m cube.mesh -o 3 -rs 3 -mid 321 -tid 5 -ls 3 -nor -lc 0.1 -pa -d cuda
57// Adapted discrete size; explicit combo of metrics; mixed tri/quad mesh:
58// mpirun -np 4 pmesh-optimizer -m ../../data/square-mixed.mesh -o 2 -rs 2 -mid 2 -tid 5 -ni 200 -bnd -qo 6 -cmb 2 -nor
59// Adapted discrete size+aspect_ratio:
60// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100
61// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 7 -tid 6 -ni 100 -qo 6 -ex -st 1 -nor
62// Adapted discrete size+orientation:
63// mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 36 -tid 8 -qo 4 -fd -nor
64// Adapted discrete aspect ratio (3D):
65// mpirun -np 4 pmesh-optimizer -m cube.mesh -o 2 -rs 2 -mid 302 -tid 7 -ni 20 -bnd -qt 1 -qo 8
66//
67// Adaptive limiting:
68// mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5
69// Adaptive limiting through the L-BFGS solver:
70// mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 400 -qo 5 -nor -vl 1 -alc 0.5 -st 1
71// Adaptive limiting through FD (requires GSLIB):
72// * mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
73//
74// Blade shape:
75// mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
76// (requires CUDA):
77// * mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8 -d cuda
78// Blade shape with FD-based solver:
79// mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
80// Blade limited shape:
81// mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
82// ICF shape and equal size:
83// mpirun -np 4 pmesh-optimizer -o 3 -mid 80 -bec -tid 2 -ni 25 -ls 3 -art 2 -qo 5
84// ICF shape and initial size:
85// mpirun -np 4 pmesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
86// ICF shape:
87// mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
88// ICF limited shape:
89// mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8 -lc 10
90// ICF combo shape + size (rings, slow convergence):
91// mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 1000 -bnd -qt 1 -qo 8 -cmb 1
92// Mixed tet / cube / hex mesh with limiting:
93// mpirun -np 4 pmesh-optimizer -m ../../data/fichera-mixed-p2.mesh -o 4 -rs 1 -mid 301 -tid 1 -fix-bnd -qo 6 -nor -lc 0.25
94// 3D pinched sphere shape (the mesh is in the mfem/data GitHub repository):
95// * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/ball-pert.mesh -o 4 -mid 303 -tid 1 -ni 20 -li 500 -fix-bnd
96// 2D non-conforming shape and equal size:
97// mpirun -np 4 pmesh-optimizer -m ./amr-quad-q2.mesh -o 2 -rs 1 -mid 9 -tid 2 -ni 200 -bnd -qt 1 -qo 8
98//
99// 2D untangling:
100// mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 22 -tid 1 -ni 50 -li 50 -qo 4 -fd -vl 1
101// 2D untangling with shifted barrier metric:
102// mpirun -np 4 pmesh-optimizer -m jagged.mesh -o 2 -mid 4 -tid 1 -ni 50 -qo 4 -fd -vl 1 -btype 1
103// 3D untangling (the mesh is in the mfem/data GitHub repository):
104// * mpirun -np 4 pmesh-optimizer -m ../../../mfem_data/cube-holes-inv.mesh -o 3 -mid 313 -tid 1 -rtol 1e-5 -li 50 -qo 4 -fd -vl 1
105// Shape optimization for a Kershaw transformed mesh using partial assembly:
106// Mesh for Kershaw transformation must be a Cartesian mesh with nx % 6 = ny % 2 = nz % 2 = 0.
107// Kershaw transformation can be imposed using the transformation ('t') feature in the mesh-explorer miniapp.
108// * mpirun - np 6 pmesh-optimizer -m kershaw-24x24x24.mesh -mid 303 -tid 1 -bnd -ni 100 -art 1 -ls 3 -qo 8 -li 40 -o 2 -qo 8 -ker -pa
109
110#include "mfem.hpp"
111#include "../common/mfem-common.hpp"
112#include <iostream>
113#include <fstream>
114#include "mesh-optimizer.hpp"
115
116using namespace mfem;
117using namespace std;
118
119int main (int argc, char *argv[])
120{
121 // 0. Initialize MPI and HYPRE.
122 Mpi::Init(argc, argv);
123 int myid = Mpi::WorldRank();
124 Hypre::Init();
125
126 // 1. Set the method's default parameters.
127 const char *mesh_file = "icf.mesh";
128 int mesh_poly_deg = 1;
129 int rs_levels = 0;
130 int rp_levels = 0;
131 real_t jitter = 0.0;
132 int metric_id = 1;
133 int target_id = 1;
134 real_t lim_const = 0.0;
135 real_t adapt_lim_const = 0.0;
136 int quad_type = 1;
137 int quad_order = 8;
138 int solver_type = 0;
139 int solver_iter = 20;
140#ifdef MFEM_USE_SINGLE
141 real_t solver_rtol = 1e-4;
142#else
143 real_t solver_rtol = 1e-10;
144#endif
145 int solver_art_type = 0;
146 int lin_solver = 2;
147 int max_lin_iter = 100;
148 bool move_bnd = true;
149 int combomet = 0;
150 bool bal_expl_combo = false;
151 bool hradaptivity = false;
152 int h_metric_id = -1;
153 bool normalization = false;
154 bool visualization = true;
155 int verbosity_level = 0;
156 bool fdscheme = false;
157 int adapt_eval = 0;
158 bool exactaction = false;
159 bool integ_over_targ = true;
160 const char *devopt = "cpu";
161 bool pa = false;
162 int n_hr_iter = 5;
163 int n_h_iter = 1;
164 int mesh_node_ordering = 0;
165 int barrier_type = 0;
166 int worst_case_type = 0;
167
168 // 2. Parse command-line options.
169 OptionsParser args(argc, argv);
170 args.AddOption(&mesh_file, "-m", "--mesh",
171 "Mesh file to use.");
172 args.AddOption(&mesh_poly_deg, "-o", "--order",
173 "Polynomial degree of mesh finite element space.");
174 args.AddOption(&rs_levels, "-rs", "--refine-serial",
175 "Number of times to refine the mesh uniformly in serial.");
176 args.AddOption(&rp_levels, "-rp", "--refine-parallel",
177 "Number of times to refine the mesh uniformly in parallel.");
178 args.AddOption(&jitter, "-ji", "--jitter",
179 "Random perturbation scaling factor.");
180 args.AddOption(&metric_id, "-mid", "--metric-id",
181 "Mesh optimization metric:\n\t"
182 "T-metrics\n\t"
183 "1 : |T|^2 -- 2D no type\n\t"
184 "2 : 0.5|T|^2/tau-1 -- 2D shape (condition number)\n\t"
185 "7 : |T-T^-t|^2 -- 2D shape+size\n\t"
186 "9 : tau*|T-T^-t|^2 -- 2D shape+size\n\t"
187 "14 : |T-I|^2 -- 2D shape+size+orientation\n\t"
188 "22 : 0.5(|T|^2-2*tau)/(tau-tau_0) -- 2D untangling\n\t"
189 "50 : 0.5|T^tT|^2/tau^2-1 -- 2D shape\n\t"
190 "55 : (tau-1)^2 -- 2D size\n\t"
191 "56 : 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 2D size\n\t"
192 "58 : |T^tT|^2/(tau^2)-2*|T|^2/tau+2 -- 2D shape\n\t"
193 "77 : 0.5(tau-1/tau)^2 -- 2D size\n\t"
194 "80 : (1-gamma)mu_2 + gamma mu_77 -- 2D shape+size\n\t"
195 "85 : |T-|T|/sqrt(2)I|^2 -- 2D shape+orientation\n\t"
196 "90 : balanced combo mu_50 & mu_77 -- 2D shape+size\n\t"
197 "94 : balanced combo mu_2 & mu_56 -- 2D shape+size\n\t"
198 "98 : (1/tau)|T-I|^2 -- 2D shape+size+orientation\n\t"
199 // "211: (tau-1)^2-tau+sqrt(tau^2+eps) -- 2D untangling\n\t"
200 // "252: 0.5(tau-1)^2/(tau-tau_0) -- 2D untangling\n\t"
201 "301: (|T||T^-1|)/3-1 -- 3D shape\n\t"
202 "302: (|T|^2|T^-1|^2)/9-1 -- 3D shape\n\t"
203 "303: (|T|^2)/3/tau^(2/3)-1 -- 3D shape\n\t"
204 "304: (|T|^3)/3^{3/2}/tau-1 -- 3D shape\n\t"
205 // "311: (tau-1)^2-tau+sqrt(tau^2+eps)-- 3D untangling\n\t"
206 "313: (|T|^2)(tau-tau0)^(-2/3)/3 -- 3D untangling\n\t"
207 "315: (tau-1)^2 -- 3D no type\n\t"
208 "316: 0.5(sqrt(tau)-1/sqrt(tau))^2 -- 3D no type\n\t"
209 "321: |T-T^-t|^2 -- 3D shape+size\n\t"
210 "322: |T-adjT^-t|^2 -- 3D shape+size\n\t"
211 "323: |J|^3-3sqrt(3)ln(det(J))-3sqrt(3) -- 3D shape+size\n\t"
212 "328: balanced combo mu_301 & mu_316 -- 3D shape+size\n\t"
213 "332: (1-gamma) mu_302 + gamma mu_315 -- 3D shape+size\n\t"
214 "333: (1-gamma) mu_302 + gamma mu_316 -- 3D shape+size\n\t"
215 "334: (1-gamma) mu_303 + gamma mu_316 -- 3D shape+size\n\t"
216 "328: balanced combo mu_302 & mu_318 -- 3D shape+size\n\t"
217 "347: (1-gamma) mu_304 + gamma mu_316 -- 3D shape+size\n\t"
218 // "352: 0.5(tau-1)^2/(tau-tau_0) -- 3D untangling\n\t"
219 "360: (|T|^3)/3^{3/2}-tau -- 3D shape\n\t"
220 "A-metrics\n\t"
221 "11 : (1/4*alpha)|A-(adjA)^T(W^TW)/omega|^2 -- 2D shape\n\t"
222 "36 : (1/alpha)|A-W|^2 -- 2D shape+size+orientation\n\t"
223 "107: (1/2*alpha)|A-|A|/|W|W|^2 -- 2D shape+orientation\n\t"
224 "126: (1-gamma)nu_11 + gamma*nu_14a -- 2D shape+size\n\t"
225 );
226 args.AddOption(&target_id, "-tid", "--target-id",
227 "Target (ideal element) type:\n\t"
228 "1: Ideal shape, unit size\n\t"
229 "2: Ideal shape, equal size\n\t"
230 "3: Ideal shape, initial size\n\t"
231 "4: Given full analytic Jacobian (in physical space)\n\t"
232 "5: Ideal shape, given size (in physical space)");
233 args.AddOption(&lim_const, "-lc", "--limit-const", "Limiting constant.");
234 args.AddOption(&adapt_lim_const, "-alc", "--adapt-limit-const",
235 "Adaptive limiting coefficient constant.");
236 args.AddOption(&quad_type, "-qt", "--quad-type",
237 "Quadrature rule type:\n\t"
238 "1: Gauss-Lobatto\n\t"
239 "2: Gauss-Legendre\n\t"
240 "3: Closed uniform points");
241 args.AddOption(&quad_order, "-qo", "--quad_order",
242 "Order of the quadrature rule.");
243 args.AddOption(&solver_type, "-st", "--solver-type",
244 " Type of solver: (default) 0: Newton, 1: LBFGS");
245 args.AddOption(&solver_iter, "-ni", "--newton-iters",
246 "Maximum number of Newton iterations.");
247 args.AddOption(&solver_rtol, "-rtol", "--newton-rel-tolerance",
248 "Relative tolerance for the Newton solver.");
249 args.AddOption(&solver_art_type, "-art", "--adaptive-rel-tol",
250 "Type of adaptive relative linear solver tolerance:\n\t"
251 "0: None (default)\n\t"
252 "1: Eisenstat-Walker type 1\n\t"
253 "2: Eisenstat-Walker type 2");
254 args.AddOption(&lin_solver, "-ls", "--lin-solver",
255 "Linear solver:\n\t"
256 "0: l1-Jacobi\n\t"
257 "1: CG\n\t"
258 "2: MINRES\n\t"
259 "3: MINRES + Jacobi preconditioner\n\t"
260 "4: MINRES + l1-Jacobi preconditioner");
261 args.AddOption(&max_lin_iter, "-li", "--lin-iter",
262 "Maximum number of iterations in the linear solve.");
263 args.AddOption(&move_bnd, "-bnd", "--move-boundary", "-fix-bnd",
264 "--fix-boundary",
265 "Enable motion along horizontal and vertical boundaries.");
266 args.AddOption(&combomet, "-cmb", "--combo-type",
267 "Combination of metrics options:\n\t"
268 "0: Use single metric\n\t"
269 "1: Shape + space-dependent size given analytically\n\t"
270 "2: Shape + adapted size given discretely; shared target");
271 args.AddOption(&bal_expl_combo, "-bec", "--balance-explicit-combo",
272 "-no-bec", "--balance-explicit-combo",
273 "Automatic balancing of explicit combo metrics.");
274 args.AddOption(&hradaptivity, "-hr", "--hr-adaptivity", "-no-hr",
275 "--no-hr-adaptivity",
276 "Enable hr-adaptivity.");
277 args.AddOption(&h_metric_id, "-hmid", "--h-metric",
278 "Same options as metric_id. Used to determine refinement"
279 " type for each element if h-adaptivity is enabled.");
280 args.AddOption(&normalization, "-nor", "--normalization", "-no-nor",
281 "--no-normalization",
282 "Make all terms in the optimization functional unitless.");
283 args.AddOption(&fdscheme, "-fd", "--fd_approximation",
284 "-no-fd", "--no-fd-approx",
285 "Enable finite difference based derivative computations.");
286 args.AddOption(&exactaction, "-ex", "--exact_action",
287 "-no-ex", "--no-exact-action",
288 "Enable exact action of TMOP_Integrator.");
289 args.AddOption(&integ_over_targ, "-it", "--integrate-target",
290 "-ir", "--integrate-reference",
291 "Integrate over target (-it) or reference (-ir) element.");
292 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
293 "--no-visualization",
294 "Enable or disable GLVis visualization.");
295 args.AddOption(&verbosity_level, "-vl", "--verbosity-level",
296 "Verbosity level for the involved iterative solvers:\n\t"
297 "0: no output\n\t"
298 "1: Newton iterations\n\t"
299 "2: Newton iterations + linear solver summaries\n\t"
300 "3: newton iterations + linear solver iterations");
301 args.AddOption(&adapt_eval, "-ae", "--adaptivity-evaluator",
302 "0 - Advection based (DEFAULT), 1 - GSLIB.");
303 args.AddOption(&devopt, "-d", "--device",
304 "Device configuration string, see Device::Configure().");
305 args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
306 "--no-partial-assembly", "Enable Partial Assembly.");
307 args.AddOption(&n_hr_iter, "-nhr", "--n_hr_iter",
308 "Number of hr-adaptivity iterations.");
309 args.AddOption(&n_h_iter, "-nh", "--n_h_iter",
310 "Number of h-adaptivity iterations per r-adaptivity"
311 "iteration.");
312 args.AddOption(&mesh_node_ordering, "-mno", "--mesh_node_ordering",
313 "Ordering of mesh nodes."
314 "0 (default): byNodes, 1: byVDIM");
315 args.AddOption(&barrier_type, "-btype", "--barrier-type",
316 "0 - None,"
317 "1 - Shifted Barrier,"
318 "2 - Pseudo Barrier.");
319 args.AddOption(&worst_case_type, "-wctype", "--worst-case-type",
320 "0 - None,"
321 "1 - Beta,"
322 "2 - PMean.");
323
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 if (h_metric_id < 0) { h_metric_id = metric_id; }
332
333 if (hradaptivity)
334 {
335 MFEM_VERIFY(strcmp(devopt,"cpu")==0, "HR-adaptivity is currently only"
336 " supported on cpus.");
337 }
338 Device device(devopt);
339 if (myid == 0) { device.Print();}
340
341 // 3. Initialize and refine the starting mesh.
342 Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
343 for (int lev = 0; lev < rs_levels; lev++)
344 {
345 mesh->UniformRefinement();
346 }
347 const int dim = mesh->Dimension();
348
349 if (hradaptivity) { mesh->EnsureNCMesh(); }
350 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
351 delete mesh;
352 for (int lev = 0; lev < rp_levels; lev++)
353 {
354 pmesh->UniformRefinement();
355 }
356
357 // 4. Define a finite element space on the mesh. Here we use vector finite
358 // elements which are tensor products of quadratic finite elements. The
359 // number of components in the vector finite element space is specified by
360 // the last parameter of the FiniteElementSpace constructor.
362 if (mesh_poly_deg <= 0)
363 {
364 fec = new QuadraticPosFECollection;
365 mesh_poly_deg = 2;
366 }
367 else { fec = new H1_FECollection(mesh_poly_deg, dim); }
368 ParFiniteElementSpace *pfespace = new ParFiniteElementSpace(pmesh, fec, dim,
369 mesh_node_ordering);
370
371 // 5. Make the mesh curved based on the above finite element space. This
372 // means that we define the mesh elements through a fespace-based
373 // transformation of the reference element.
374 pmesh->SetNodalFESpace(pfespace);
375
376 // 7. Get the mesh nodes (vertices and other degrees of freedom in the finite
377 // element space) as a finite element grid function in fespace. Note that
378 // changing x automatically changes the shapes of the mesh elements.
379 ParGridFunction x(pfespace);
380 pmesh->SetNodalGridFunction(&x);
381
382 // 8. Define a vector representing the minimal local mesh size in the mesh
383 // nodes. We index the nodes using the scalar version of the degrees of
384 // freedom in pfespace. Note: this is partition-dependent.
385 //
386 // In addition, compute average mesh size and total volume.
387 Vector h0(pfespace->GetNDofs());
388 h0 = infinity();
389 real_t vol_loc = 0.0;
390 Array<int> dofs;
391 for (int i = 0; i < pmesh->GetNE(); i++)
392 {
393 // Get the local scalar element degrees of freedom in dofs.
394 pfespace->GetElementDofs(i, dofs);
395 // Adjust the value of h0 in dofs based on the local mesh size.
396 const real_t hi = pmesh->GetElementSize(i);
397 for (int j = 0; j < dofs.Size(); j++)
398 {
399 h0(dofs[j]) = min(h0(dofs[j]), hi);
400 }
401 vol_loc += pmesh->GetElementVolume(i);
402 }
403 real_t vol_glb;
404 MPI_Allreduce(&vol_loc, &vol_glb, 1, MPITypeMap<real_t>::mpi_type,
405 MPI_SUM, MPI_COMM_WORLD);
406 const real_t small_phys_size = pow(vol_glb, 1.0 / dim) / 100.0;
407
408 // 9. Add a random perturbation to the nodes in the interior of the domain.
409 // We define a random grid function of fespace and make sure that it is
410 // zero on the boundary and its values are locally of the order of h0.
411 // The latter is based on the DofToVDof() method which maps the scalar to
412 // the vector degrees of freedom in pfespace.
413 ParGridFunction rdm(pfespace);
414 rdm.Randomize();
415 rdm -= 0.25; // Shift to random values in [-0.5,0.5].
416 rdm *= jitter;
417 rdm.HostReadWrite();
418 // Scale the random values to be of order of the local mesh size.
419 for (int i = 0; i < pfespace->GetNDofs(); i++)
420 {
421 for (int d = 0; d < dim; d++)
422 {
423 rdm(pfespace->DofToVDof(i,d)) *= h0(i);
424 }
425 }
426 Array<int> vdofs;
427 for (int i = 0; i < pfespace->GetNBE(); i++)
428 {
429 // Get the vector degrees of freedom in the boundary element.
430 pfespace->GetBdrElementVDofs(i, vdofs);
431 // Set the boundary values to zero.
432 for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
433 }
434 x -= rdm;
435 // Set the perturbation of all nodes from the true nodes.
436 x.SetTrueVector();
438
439 // 10. Save the starting (prior to the optimization) mesh to a file. This
440 // output can be viewed later using GLVis: "glvis -m perturbed -np
441 // num_mpi_tasks".
442 {
443 ostringstream mesh_name;
444 mesh_name << "perturbed.mesh";
445 ofstream mesh_ofs(mesh_name.str().c_str());
446 mesh_ofs.precision(8);
447 pmesh->PrintAsOne(mesh_ofs);
448 }
449
450 // 11. Store the starting (prior to the optimization) positions.
451 ParGridFunction x0(pfespace);
452 x0 = x;
453
454 // 12. Form the integrator that uses the chosen metric and target.
455 real_t min_detJ = -0.1;
456 TMOP_QualityMetric *metric = NULL;
457 switch (metric_id)
458 {
459 // T-metrics
460 case 1: metric = new TMOP_Metric_001; break;
461 case 2: metric = new TMOP_Metric_002; break;
462 case 4: metric = new TMOP_Metric_004; break;
463 case 7: metric = new TMOP_Metric_007; break;
464 case 9: metric = new TMOP_Metric_009; break;
465 case 14: metric = new TMOP_Metric_014; break;
466 case 22: metric = new TMOP_Metric_022(min_detJ); break;
467 case 50: metric = new TMOP_Metric_050; break;
468 case 55: metric = new TMOP_Metric_055; break;
469 case 56: metric = new TMOP_Metric_056; break;
470 case 58: metric = new TMOP_Metric_058; break;
471 case 66: metric = new TMOP_Metric_066(0.5); break;
472 case 77: metric = new TMOP_Metric_077; break;
473 case 80: metric = new TMOP_Metric_080(0.5); break;
474 case 85: metric = new TMOP_Metric_085; break;
475 case 90: metric = new TMOP_Metric_090; break;
476 case 94: metric = new TMOP_Metric_094; break;
477 case 98: metric = new TMOP_Metric_098; break;
478 // case 211: metric = new TMOP_Metric_211; break;
479 // case 252: metric = new TMOP_Metric_252(min_detJ); break;
480 case 301: metric = new TMOP_Metric_301; break;
481 case 302: metric = new TMOP_Metric_302; break;
482 case 303: metric = new TMOP_Metric_303; break;
483 case 304: metric = new TMOP_Metric_304; break;
484 // case 311: metric = new TMOP_Metric_311; break;
485 case 313: metric = new TMOP_Metric_313(min_detJ); break;
486 case 315: metric = new TMOP_Metric_315; break;
487 case 316: metric = new TMOP_Metric_316; break;
488 case 321: metric = new TMOP_Metric_321; break;
489 case 322: metric = new TMOP_Metric_322; break;
490 case 323: metric = new TMOP_Metric_323; break;
491 case 328: metric = new TMOP_Metric_328; break;
492 case 332: metric = new TMOP_Metric_332(0.5); break;
493 case 333: metric = new TMOP_Metric_333(0.5); break;
494 case 334: metric = new TMOP_Metric_334(0.5); break;
495 case 338: metric = new TMOP_Metric_338; break;
496 case 347: metric = new TMOP_Metric_347(0.5); break;
497 // case 352: metric = new TMOP_Metric_352(min_detJ); break;
498 case 360: metric = new TMOP_Metric_360; break;
499 // A-metrics
500 case 11: metric = new TMOP_AMetric_011; break;
501 case 36: metric = new TMOP_AMetric_036; break;
502 case 107: metric = new TMOP_AMetric_107a; break;
503 case 126: metric = new TMOP_AMetric_126(0.9); break;
504 default:
505 if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
506 return 3;
507 }
508 TMOP_QualityMetric *h_metric = NULL;
509 if (hradaptivity)
510 {
511 switch (h_metric_id)
512 {
513 case 1: h_metric = new TMOP_Metric_001; break;
514 case 2: h_metric = new TMOP_Metric_002; break;
515 case 7: h_metric = new TMOP_Metric_007; break;
516 case 9: h_metric = new TMOP_Metric_009; break;
517 case 55: h_metric = new TMOP_Metric_055; break;
518 case 56: h_metric = new TMOP_Metric_056; break;
519 case 58: h_metric = new TMOP_Metric_058; break;
520 case 77: h_metric = new TMOP_Metric_077; break;
521 case 315: h_metric = new TMOP_Metric_315; break;
522 case 316: h_metric = new TMOP_Metric_316; break;
523 case 321: h_metric = new TMOP_Metric_321; break;
524 default: cout << "Metric_id not supported for h-adaptivity: " << h_metric_id <<
525 endl;
526 return 3;
527 }
528 }
529
531 switch (barrier_type)
532 {
533 case 0: btype = TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::None;
534 break;
535 case 1: btype = TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Shifted;
536 break;
537 case 2: btype = TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Pseudo;
538 break;
539 default: cout << "barrier_type not supported: " << barrier_type << endl;
540 return 3;
541 }
542
544 switch (worst_case_type)
545 {
546 case 0: wctype = TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::None;
547 break;
548 case 1: wctype = TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::Beta;
549 break;
550 case 2: wctype = TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::PMean;
551 break;
552 default: cout << "worst_case_type not supported: " << worst_case_type << endl;
553 return 3;
554 }
555
556 TMOP_QualityMetric *untangler_metric = NULL;
557 if (barrier_type > 0 || worst_case_type > 0)
558 {
559 if (barrier_type > 0)
560 {
561 MFEM_VERIFY(metric_id == 4 || metric_id == 14 || metric_id == 66,
562 "Metric not supported for shifted/pseudo barriers.");
563 }
564 untangler_metric = new TMOP_WorstCaseUntangleOptimizer_Metric(*metric,
565 2,
566 1.5,
567 0.001,//0.01 for pseudo barrier
568 0.001,
569 btype,
570 wctype);
571 }
572
573 if (metric_id < 300 || h_metric_id < 300)
574 {
575 MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
576 }
577 if (metric_id >= 300 || h_metric_id >= 300)
578 {
579 MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
580 }
581
583 TargetConstructor *target_c = NULL;
584 HessianCoefficient *adapt_coeff = NULL;
585 HRHessianCoefficient *hr_adapt_coeff = NULL;
586 H1_FECollection ind_fec(mesh_poly_deg, dim);
587 ParFiniteElementSpace ind_fes(pmesh, &ind_fec);
588 ParFiniteElementSpace ind_fesv(pmesh, &ind_fec, dim);
589 ParGridFunction size(&ind_fes), aspr(&ind_fes), ori(&ind_fes);
590 ParGridFunction aspr3d(&ind_fesv);
591
592 const AssemblyLevel al =
593 pa ? AssemblyLevel::PARTIAL : AssemblyLevel::LEGACY;
594
595 switch (target_id)
596 {
597 case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
598 case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
599 case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
600 case 4:
601 {
603 AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
604 adapt_coeff = new HessianCoefficient(dim, metric_id);
605 tc->SetAnalyticTargetSpec(NULL, NULL, adapt_coeff);
606 target_c = tc;
607 break;
608 }
609 case 5: // Discrete size 2D or 3D
610 {
612 DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
613 if (adapt_eval == 0)
614 {
616 }
617 else
618 {
619#ifdef MFEM_USE_GSLIB
621#else
622 MFEM_ABORT("MFEM is not built with GSLIB.");
623#endif
624 }
625 ConstructSizeGF(size);
626 tc->SetParDiscreteTargetSize(size);
627 target_c = tc;
628 break;
629 }
630 case 6: // material indicator 2D
631 {
632 ParGridFunction d_x(&ind_fes), d_y(&ind_fes), disc(&ind_fes);
633
635 DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
637 disc.ProjectCoefficient(mat_coeff);
638 if (adapt_eval == 0)
639 {
641 }
642 else
643 {
644#ifdef MFEM_USE_GSLIB
646#else
647 MFEM_ABORT("MFEM is not built with GSLIB.");
648#endif
649 }
650 // Diffuse the interface
652
653 // Get partials with respect to x and y of the grid function
654 disc.GetDerivative(1,0,d_x);
655 disc.GetDerivative(1,1,d_y);
656
657 // Compute the squared magnitude of the gradient
658 for (int i = 0; i < size.Size(); i++)
659 {
660 size(i) = std::pow(d_x(i),2)+std::pow(d_y(i),2);
661 }
662 const real_t max = size.Max();
663 real_t max_all;
664 MPI_Allreduce(&max, &max_all, 1, MPITypeMap<real_t>::mpi_type,
665 MPI_MAX, MPI_COMM_WORLD);
666
667 for (int i = 0; i < d_x.Size(); i++)
668 {
669 d_x(i) = std::abs(d_x(i));
670 d_y(i) = std::abs(d_y(i));
671 }
672 const real_t eps = 0.01;
673 const real_t aspr_ratio = 20.0;
674 const real_t size_ratio = 40.0;
675
676 for (int i = 0; i < size.Size(); i++)
677 {
678 size(i) = (size(i)/max_all);
679 aspr(i) = (d_x(i)+eps)/(d_y(i)+eps);
680 aspr(i) = 0.1 + 0.9*(1-size(i))*(1-size(i));
681 if (aspr(i) > aspr_ratio) {aspr(i) = aspr_ratio;}
682 if (aspr(i) < 1.0/aspr_ratio) {aspr(i) = 1.0/aspr_ratio;}
683 }
684 Vector vals;
685 const int NE = pmesh->GetNE();
686 real_t volume = 0.0, volume_ind = 0.0;
687
688 for (int i = 0; i < NE; i++)
689 {
691 const IntegrationRule &ir =
692 IntRules.Get(pmesh->GetElementBaseGeometry(i), Tr->OrderJ());
693 size.GetValues(i, ir, vals);
694 for (int j = 0; j < ir.GetNPoints(); j++)
695 {
696 const IntegrationPoint &ip = ir.IntPoint(j);
697 Tr->SetIntPoint(&ip);
698 volume += ip.weight * Tr->Weight();
699 volume_ind += vals(j) * ip.weight * Tr->Weight();
700 }
701 }
702 real_t volume_all, volume_ind_all;
703 MPI_Allreduce(&volume, &volume_all, 1, MPITypeMap<real_t>::mpi_type,
704 MPI_SUM, MPI_COMM_WORLD);
705 MPI_Allreduce(&volume_ind, &volume_ind_all, 1,
706 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
707 const int NE_ALL = pmesh->GetGlobalNE();
708
709 const real_t avg_zone_size = volume_all / NE_ALL;
710
711 const real_t small_avg_ratio =
712 (volume_ind_all + (volume_all - volume_ind_all) / size_ratio)
713 / volume_all;
714
715 const real_t small_zone_size = small_avg_ratio * avg_zone_size;
716 const real_t big_zone_size = size_ratio * small_zone_size;
717
718 for (int i = 0; i < size.Size(); i++)
719 {
720 const real_t val = size(i);
721 const real_t a = (big_zone_size - small_zone_size) / small_zone_size;
722 size(i) = big_zone_size / (1.0+a*val);
723 }
724
725 DiffuseField(size, 2);
726 DiffuseField(aspr, 2);
727
728 tc->SetParDiscreteTargetSize(size);
730 target_c = tc;
731 break;
732 }
733 case 7: // Discrete aspect ratio 3D
734 {
736 DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
737 if (adapt_eval == 0)
738 {
740 }
741 else
742 {
743#ifdef MFEM_USE_GSLIB
745#else
746 MFEM_ABORT("MFEM is not built with GSLIB.");
747#endif
748 }
750 aspr3d.ProjectCoefficient(fd_aspr3d);
752 target_c = tc;
753 break;
754 }
755 case 8: // shape/size + orientation 2D
756 {
758 DiscreteAdaptTC *tc = new DiscreteAdaptTC(target_t);
759 if (adapt_eval == 0)
760 {
762 }
763 else
764 {
765#ifdef MFEM_USE_GSLIB
767#else
768 MFEM_ABORT("MFEM is not built with GSLIB.");
769#endif
770 }
771
772 ConstantCoefficient size_coeff(0.1*0.1);
773 size.ProjectCoefficient(size_coeff);
774 tc->SetParDiscreteTargetSize(size);
775
777 ori.ProjectCoefficient(ori_coeff);
779 target_c = tc;
780 break;
781 }
782 // Targets used for hr-adaptivity tests.
783 case 9: // size target in an annular region.
784 case 10: // size+aspect-ratio in an annular region.
785 case 11: // size+aspect-ratio target for a rotate sine wave
786 {
788 AnalyticAdaptTC *tc = new AnalyticAdaptTC(target_t);
789 hr_adapt_coeff = new HRHessianCoefficient(dim, target_id - 9);
790 tc->SetAnalyticTargetSpec(NULL, NULL, hr_adapt_coeff);
791 target_c = tc;
792 break;
793 }
794 default:
795 if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
796 return 3;
797 }
798 if (target_c == NULL)
799 {
800 target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
801 }
802 target_c->SetNodes(x0);
803
804 // Automatically balanced gamma in composite metrics.
805 auto metric_combo = dynamic_cast<TMOP_Combo_QualityMetric *>(metric);
806 if (metric_combo && bal_expl_combo)
807 {
808 Vector bal_weights;
809 metric_combo->ComputeBalancedWeights(x, *target_c, bal_weights);
810 metric_combo->SetWeights(bal_weights);
811 }
812
813 TMOP_QualityMetric *metric_to_use = barrier_type > 0 || worst_case_type > 0
814 ? untangler_metric
815 : metric;
816 auto tmop_integ = new TMOP_Integrator(metric_to_use, target_c, h_metric);
817 tmop_integ->IntegrateOverTarget(integ_over_targ);
818 if (barrier_type > 0 || worst_case_type > 0)
819 {
820 tmop_integ->ComputeUntangleMetricQuantiles(x, *pfespace);
821 }
822
823 // Finite differences for computations of derivatives.
824 if (fdscheme)
825 {
826 MFEM_VERIFY(pa == false, "PA for finite differences is not implemented.");
827 tmop_integ->EnableFiniteDifferences(x);
828 }
829 tmop_integ->SetExactActionFlag(exactaction);
830
831 // Setup the quadrature rules for the TMOP integrator.
832 IntegrationRules *irules = NULL;
833 switch (quad_type)
834 {
835 case 1: irules = &IntRulesLo; break;
836 case 2: irules = &IntRules; break;
837 case 3: irules = &IntRulesCU; break;
838 default:
839 if (myid == 0) { cout << "Unknown quad_type: " << quad_type << endl; }
840 return 3;
841 }
842 tmop_integ->SetIntegrationRules(*irules, quad_order);
843 if (myid == 0 && dim == 2)
844 {
845 cout << "Triangle quadrature points: "
846 << irules->Get(Geometry::TRIANGLE, quad_order).GetNPoints()
847 << "\nQuadrilateral quadrature points: "
848 << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
849 }
850 if (myid == 0 && dim == 3)
851 {
852 cout << "Tetrahedron quadrature points: "
853 << irules->Get(Geometry::TETRAHEDRON, quad_order).GetNPoints()
854 << "\nHexahedron quadrature points: "
855 << irules->Get(Geometry::CUBE, quad_order).GetNPoints()
856 << "\nPrism quadrature points: "
857 << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
858 }
859
860 // Limit the node movement.
861 // The limiting distances can be given by a general function of space.
862 ParFiniteElementSpace dist_pfespace(pmesh, fec); // scalar space
863 ParGridFunction dist(&dist_pfespace);
864 dist = 1.0;
865 // The small_phys_size is relevant only with proper normalization.
866 if (normalization) { dist = small_phys_size; }
867 ConstantCoefficient lim_coeff(lim_const);
868 if (lim_const != 0.0) { tmop_integ->EnableLimiting(x0, dist, lim_coeff); }
869
870 // Adaptive limiting.
871 ParGridFunction adapt_lim_gf0(&ind_fes);
872 ConstantCoefficient adapt_lim_coeff(adapt_lim_const);
873 AdaptivityEvaluator *adapt_lim_eval = NULL;
874 if (adapt_lim_const > 0.0)
875 {
876 MFEM_VERIFY(pa == false, "PA is not implemented for adaptive limiting");
877
878 FunctionCoefficient adapt_lim_gf0_coeff(adapt_lim_fun);
879 adapt_lim_gf0.ProjectCoefficient(adapt_lim_gf0_coeff);
880
881 if (adapt_eval == 0) { adapt_lim_eval = new AdvectorCG(al); }
882 else if (adapt_eval == 1)
883 {
884#ifdef MFEM_USE_GSLIB
885 adapt_lim_eval = new InterpolatorFP;
886#else
887 MFEM_ABORT("MFEM is not built with GSLIB support!");
888#endif
889 }
890 else { MFEM_ABORT("Bad interpolation option."); }
891
892 tmop_integ->EnableAdaptiveLimiting(adapt_lim_gf0, adapt_lim_coeff,
893 *adapt_lim_eval);
894 if (visualization)
895 {
896 socketstream vis1;
897 common::VisualizeField(vis1, "localhost", 19916, adapt_lim_gf0, "Zeta 0",
898 300, 600, 300, 300);
899 }
900 }
901
902 // Has to be after the enabling of the limiting / alignment, as it computes
903 // normalization factors for these terms as well.
904 if (normalization) { tmop_integ->ParEnableNormalization(x0); }
905
906 // 13. Setup the final NonlinearForm (which defines the integral of interest,
907 // its first and second derivatives). Here we can use a combination of
908 // metrics, i.e., optimize the sum of two integrals, where both are
909 // scaled by used-defined space-dependent weights. Note that there are
910 // no command-line options for the weights and the type of the second
911 // metric; one should update those in the code.
912 ParNonlinearForm a(pfespace);
913 if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
914 ConstantCoefficient *metric_coeff1 = NULL;
915 TMOP_QualityMetric *metric2 = NULL;
916 TargetConstructor *target_c2 = NULL;
917 FunctionCoefficient metric_coeff2(weight_fun);
918
919 // Explicit combination of metrics.
920 if (combomet > 0)
921 {
922 // First metric.
923 metric_coeff1 = new ConstantCoefficient(1.0);
924 tmop_integ->SetCoefficient(*metric_coeff1);
925
926 // Second metric.
927 if (dim == 2) { metric2 = new TMOP_Metric_077; }
928 else { metric2 = new TMOP_Metric_315; }
929 TMOP_Integrator *tmop_integ2 = NULL;
930 if (combomet == 1)
931 {
932 target_c2 = new TargetConstructor(
934 target_c2->SetVolumeScale(0.01);
935 target_c2->SetNodes(x0);
936 tmop_integ2 = new TMOP_Integrator(metric2, target_c2, h_metric);
937 tmop_integ2->SetCoefficient(metric_coeff2);
938 }
939 else { tmop_integ2 = new TMOP_Integrator(metric2, target_c, h_metric); }
940 tmop_integ2->IntegrateOverTarget(integ_over_targ);
941 tmop_integ2->SetIntegrationRules(*irules, quad_order);
942 if (fdscheme) { tmop_integ2->EnableFiniteDifferences(x); }
943 tmop_integ2->SetExactActionFlag(exactaction);
944
946 combo->AddTMOPIntegrator(tmop_integ);
947 combo->AddTMOPIntegrator(tmop_integ2);
948 if (normalization) { combo->ParEnableNormalization(x0); }
949 if (lim_const != 0.0) { combo->EnableLimiting(x0, dist, lim_coeff); }
950
951 a.AddDomainIntegrator(combo);
952 }
953 else
954 {
955 a.AddDomainIntegrator(tmop_integ);
956 }
957
958 if (pa) { a.Setup(); }
959
960 // Compute the minimum det(J) of the starting mesh.
961 min_detJ = infinity();
962 const int NE = pmesh->GetNE();
963 for (int i = 0; i < NE; i++)
964 {
965 const IntegrationRule &ir =
966 irules->Get(pfespace->GetFE(i)->GetGeomType(), quad_order);
968 for (int j = 0; j < ir.GetNPoints(); j++)
969 {
970 transf->SetIntPoint(&ir.IntPoint(j));
971 min_detJ = min(min_detJ, transf->Jacobian().Det());
972 }
973 }
974 real_t minJ0;
975 MPI_Allreduce(&min_detJ, &minJ0, 1, MPITypeMap<real_t>::mpi_type,
976 MPI_MIN, MPI_COMM_WORLD);
977 min_detJ = minJ0;
978 if (myid == 0)
979 { cout << "Minimum det(J) of the original mesh is " << min_detJ << endl; }
980
981 if (min_detJ < 0.0 && barrier_type == 0
982 && metric_id != 22 && metric_id != 211 && metric_id != 252
983 && metric_id != 311 && metric_id != 313 && metric_id != 352)
984 {
985 MFEM_ABORT("The input mesh is inverted! Try an untangling metric.");
986 }
987 if (min_detJ < 0.0)
988 {
989 MFEM_VERIFY(target_t == TargetConstructor::IDEAL_SHAPE_UNIT_SIZE,
990 "Untangling is supported only for ideal targets.");
991
992 const DenseMatrix &Wideal =
994 min_detJ /= Wideal.Det();
995
996 real_t h0min = h0.Min(), h0min_all;
997 MPI_Allreduce(&h0min, &h0min_all, 1, MPITypeMap<real_t>::mpi_type,
998 MPI_MIN, MPI_COMM_WORLD);
999 // Slightly below minJ0 to avoid div by 0.
1000 min_detJ -= 0.01 * h0min_all;
1001 }
1002
1003 // For HR tests, the energy is normalized by the number of elements.
1004 const real_t init_energy = a.GetParGridFunctionEnergy(x) /
1005 (hradaptivity ? pmesh->GetGlobalNE() : 1);
1006 real_t init_metric_energy = init_energy;
1007 if (lim_const > 0.0 || adapt_lim_const > 0.0)
1008 {
1009 lim_coeff.constant = 0.0;
1010 adapt_lim_coeff.constant = 0.0;
1011 init_metric_energy = a.GetParGridFunctionEnergy(x) /
1012 (hradaptivity ? pmesh->GetGlobalNE() : 1);
1013 lim_coeff.constant = lim_const;
1014 adapt_lim_coeff.constant = adapt_lim_const;
1015 }
1016
1017 // Visualize the starting mesh and metric values.
1018 // Note that for combinations of metrics, this only shows the first metric.
1019 if (visualization)
1020 {
1021 char title[] = "Initial metric values";
1022 vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
1023 }
1024
1025 // 14. Fix all boundary nodes, or fix only a given component depending on the
1026 // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
1027 // fixed x/y/z components of the node. Attribute dim+1 corresponds to
1028 // an entirely fixed node.
1029 if (move_bnd == false)
1030 {
1031 Array<int> ess_bdr(pmesh->bdr_attributes.Max());
1032 ess_bdr = 1;
1033 a.SetEssentialBC(ess_bdr);
1034 }
1035 else
1036 {
1037 int n = 0;
1038 for (int i = 0; i < pmesh->GetNBE(); i++)
1039 {
1040 const int nd = pfespace->GetBE(i)->GetDof();
1041 const int attr = pmesh->GetBdrElement(i)->GetAttribute();
1042 MFEM_VERIFY(!(dim == 2 && attr == 3),
1043 "Boundary attribute 3 must be used only for 3D meshes. "
1044 "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
1045 "components, rest for free nodes), or use -fix-bnd.");
1046 if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
1047 if (attr == 4) { n += nd * dim; }
1048 }
1049 Array<int> ess_vdofs(n);
1050 n = 0;
1051 for (int i = 0; i < pmesh->GetNBE(); i++)
1052 {
1053 const int nd = pfespace->GetBE(i)->GetDof();
1054 const int attr = pmesh->GetBdrElement(i)->GetAttribute();
1055 pfespace->GetBdrElementVDofs(i, vdofs);
1056 if (attr == 1) // Fix x components.
1057 {
1058 for (int j = 0; j < nd; j++)
1059 { ess_vdofs[n++] = vdofs[j]; }
1060 }
1061 else if (attr == 2) // Fix y components.
1062 {
1063 for (int j = 0; j < nd; j++)
1064 { ess_vdofs[n++] = vdofs[j+nd]; }
1065 }
1066 else if (attr == 3) // Fix z components.
1067 {
1068 for (int j = 0; j < nd; j++)
1069 { ess_vdofs[n++] = vdofs[j+2*nd]; }
1070 }
1071 else if (attr == 4) // Fix all components.
1072 {
1073 for (int j = 0; j < vdofs.Size(); j++)
1074 { ess_vdofs[n++] = vdofs[j]; }
1075 }
1076 }
1077 a.SetEssentialVDofs(ess_vdofs);
1078 }
1079
1080 // As we use the inexact Newton method to solve the resulting nonlinear
1081 // system, here we setup the linear solver for the system's Jacobian.
1082 Solver *S = NULL, *S_prec = NULL;
1083#ifdef MFEM_USE_SINGLE
1084 const real_t linsol_rtol = 1e-5;
1085#else
1086 const real_t linsol_rtol = 1e-12;
1087#endif
1088 // Level of output.
1089 IterativeSolver::PrintLevel linsolver_print;
1090 if (verbosity_level == 2)
1091 { linsolver_print.Errors().Warnings().FirstAndLast(); }
1092 if (verbosity_level > 2)
1093 { linsolver_print.Errors().Warnings().Iterations(); }
1094 if (lin_solver == 0)
1095 {
1096 S = new DSmoother(1, 1.0, max_lin_iter);
1097 }
1098 else if (lin_solver == 1)
1099 {
1100 CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
1101 cg->SetMaxIter(max_lin_iter);
1102 cg->SetRelTol(linsol_rtol);
1103 cg->SetAbsTol(0.0);
1104 cg->SetPrintLevel(linsolver_print);
1105 S = cg;
1106 }
1107 else
1108 {
1109 MINRESSolver *minres = new MINRESSolver(MPI_COMM_WORLD);
1110 minres->SetMaxIter(max_lin_iter);
1111 minres->SetRelTol(linsol_rtol);
1112 minres->SetAbsTol(0.0);
1113 minres->SetPrintLevel(linsolver_print);
1114 if (lin_solver == 3 || lin_solver == 4)
1115 {
1116 if (pa)
1117 {
1118 MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
1119 auto js = new OperatorJacobiSmoother;
1120 js->SetPositiveDiagonal(true);
1121 S_prec = js;
1122 }
1123 else
1124 {
1125 auto hs = new HypreSmoother;
1126 hs->SetType((lin_solver == 3) ? HypreSmoother::Jacobi
1127 /* */ : HypreSmoother::l1Jacobi, 1);
1128 hs->SetPositiveDiagonal(true);
1129 S_prec = hs;
1130 }
1131 minres->SetPreconditioner(*S_prec);
1132 }
1133 S = minres;
1134 }
1135
1136 //
1137 // Perform the nonlinear optimization.
1138 //
1139 const IntegrationRule &ir =
1140 irules->Get(pfespace->GetFE(0)->GetGeomType(), quad_order);
1141 TMOPNewtonSolver solver(pfespace->GetComm(), ir, solver_type);
1142 // Provide all integration rules in case of a mixed mesh.
1143 solver.SetIntegrationRules(*irules, quad_order);
1144 // Specify linear solver when we use a Newton-based solver.
1145 if (solver_type == 0) { solver.SetPreconditioner(*S); }
1146 // For untangling, the solver will update the min det(T) values.
1147 solver.SetMinDetPtr(&min_detJ);
1148 solver.SetMaxIter(solver_iter);
1149 solver.SetRelTol(solver_rtol);
1150 solver.SetAbsTol(0.0);
1151 if (solver_art_type > 0)
1152 {
1153 solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
1154 }
1155 // Level of output.
1156 IterativeSolver::PrintLevel newton_print;
1157 if (verbosity_level > 0)
1158 { newton_print.Errors().Warnings().Iterations(); }
1159 solver.SetPrintLevel(newton_print);
1160 // hr-adaptivity solver.
1161 // If hr-adaptivity is disabled, r-adaptivity is done once using the
1162 // TMOPNewtonSolver.
1163 // Otherwise, "hr_iter" iterations of r-adaptivity are done followed by
1164 // "h_per_r_iter" iterations of h-adaptivity after each r-adaptivity.
1165 // The solver terminates if an h-adaptivity iteration does not modify
1166 // any element in the mesh.
1167 TMOPHRSolver hr_solver(*pmesh, a, solver,
1168 x, move_bnd, hradaptivity,
1169 mesh_poly_deg, h_metric_id,
1170 n_hr_iter, n_h_iter);
1171 hr_solver.AddGridFunctionForUpdate(&x0);
1172 if (adapt_lim_const > 0.)
1173 {
1174 hr_solver.AddGridFunctionForUpdate(&adapt_lim_gf0);
1175 hr_solver.AddFESpaceForUpdate(&ind_fes);
1176 }
1177 hr_solver.Mult();
1178
1179 // 16. Save the optimized mesh to a file. This output can be viewed later
1180 // using GLVis: "glvis -m optimized -np num_mpi_tasks".
1181 {
1182 ostringstream mesh_name;
1183 mesh_name << "optimized.mesh";
1184 ofstream mesh_ofs(mesh_name.str().c_str());
1185 mesh_ofs.precision(8);
1186 pmesh->PrintAsOne(mesh_ofs);
1187 }
1188
1189 // Report the final energy of the functional.
1190 const real_t fin_energy = a.GetParGridFunctionEnergy(x) /
1191 (hradaptivity ? pmesh->GetGlobalNE() : 1);
1192 real_t fin_metric_energy = fin_energy;
1193 if (lim_const > 0.0 || adapt_lim_const > 0.0)
1194 {
1195 lim_coeff.constant = 0.0;
1196 adapt_lim_coeff.constant = 0.0;
1197 fin_metric_energy = a.GetParGridFunctionEnergy(x) /
1198 (hradaptivity ? pmesh->GetGlobalNE() : 1);
1199 lim_coeff.constant = lim_const;
1200 adapt_lim_coeff.constant = adapt_lim_const;
1201 }
1202 if (myid == 0)
1203 {
1204 std::cout << std::scientific << std::setprecision(4);
1205 cout << "Initial strain energy: " << init_energy
1206 << " = metrics: " << init_metric_energy
1207 << " + extra terms: " << init_energy - init_metric_energy << endl;
1208 cout << " Final strain energy: " << fin_energy
1209 << " = metrics: " << fin_metric_energy
1210 << " + extra terms: " << fin_energy - fin_metric_energy << endl;
1211 cout << "The strain energy decreased by: "
1212 << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
1213 }
1214
1215 // Visualize the final mesh and metric values.
1216 if (visualization)
1217 {
1218 char title[] = "Final metric values";
1219 vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 600);
1220 }
1221
1222 if (adapt_lim_const > 0.0 && visualization)
1223 {
1224 socketstream vis0;
1225 common::VisualizeField(vis0, "localhost", 19916, adapt_lim_gf0, "Xi 0",
1226 600, 600, 300, 300);
1227 }
1228
1229 // Visualize the mesh displacement.
1230 if (visualization)
1231 {
1232 x0 -= x;
1233 socketstream sock;
1234 if (myid == 0)
1235 {
1236 sock.open("localhost", 19916);
1237 sock << "solution\n";
1238 }
1239 pmesh->PrintAsOne(sock);
1240 x0.SaveAsOne(sock);
1241 if (myid == 0)
1242 {
1243 sock << "window_title 'Displacements'\n"
1244 << "window_geometry "
1245 << 1200 << " " << 0 << " " << 600 << " " << 600 << "\n"
1246 << "keys jRmclA" << endl;
1247 }
1248 }
1249
1250 delete S;
1251 delete S_prec;
1252 delete target_c2;
1253 delete metric2;
1254 delete metric_coeff1;
1255 delete adapt_lim_eval;
1256 delete target_c;
1257 delete hr_adapt_coeff;
1258 delete adapt_coeff;
1259 delete h_metric;
1260 delete metric;
1261 delete untangler_metric;
1262 delete pfespace;
1263 delete fec;
1264 delete pmesh;
1265
1266 return 0;
1267}
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, TMOPMatrixCoefficient *mspec)
Definition: tmop.cpp:1733
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
Conjugate gradient method.
Definition: solvers.hpp:513
A coefficient that is constant across space and time.
Definition: coefficient.hpp:85
Data type for scaled Jacobi-type smoother of sparse matrix.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
real_t Det() const
Definition: densemat.cpp:535
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:286
virtual void SetParDiscreteTargetAspectRatio(const ParGridFunction &tspec_)
Definition: tmop.cpp:1899
virtual void SetParDiscreteTargetSize(const ParGridFunction &tspec_)
Definition: tmop.cpp:1879
void SetAdaptivityEvaluator(AdaptivityEvaluator *ae)
Definition: tmop.hpp:1662
virtual void SetParDiscreteTargetOrientation(const ParGridFunction &tspec_)
Definition: tmop.cpp:1909
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition: eltrans.hpp:131
virtual int OrderJ() const =0
Return the order of the elements of the Jacobian of the transformation.
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition: eltrans.hpp:119
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
int GetAttribute() const
Return element's attribute.
Definition: element.hpp:55
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:27
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition: fespace.cpp:3204
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:710
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:749
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
Definition: fespace.cpp:303
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition: fespace.cpp:249
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:326
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:329
A general function coefficient.
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:150
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:260
Parallel smoothers in hypre.
Definition: hypre.hpp:1020
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition: hypre.cpp:3550
@ l1Jacobi
l1-scaled Jacobi
Definition: hypre.hpp:1080
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition: hypre.hpp:74
Class for integration point with weight.
Definition: intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:259
Container class for integration rules.
Definition: intrules.hpp:416
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:1006
void SetRelTol(real_t rtol)
Definition: solvers.hpp:209
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void SetMaxIter(int max_it)
Definition: solvers.hpp:211
void SetAbsTol(real_t atol)
Definition: solvers.hpp:210
MINRES method.
Definition: solvers.hpp:628
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.hpp:640
Mesh data type.
Definition: mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:282
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition: mesh.hpp:1298
real_t GetElementSize(int i, int type=0)
Get the size of the i-th element relative to the perfect reference element.
Definition: mesh.cpp:107
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition: mesh.cpp:357
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:6200
real_t GetElementVolume(int i)
Definition: mesh.cpp:121
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1229
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:1255
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:10626
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10970
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1385
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void SetAdaptiveLinRtol(const int type=2, const real_t rtol0=0.5, const real_t rtol_max=0.9, const real_t alpha=0.5 *(1.0+sqrt(5.0)), const real_t gamma=1.0)
Enable adaptive linear solver relative tolerance algorithm.
Definition: solvers.cpp:1929
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:313
void SetPositiveDiagonal(bool pos_diag=true)
Replace diagonal entries with their absolute values.
Definition: solvers.hpp:349
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:29
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
void GetElementDofs(int i, Array< int > &dofs, DofTransformation &doftrans) const override
The same as GetElementDofs(), but with a user-allocated DofTransformation object. doftrans must be al...
Definition: pfespace.cpp:468
const FiniteElement * GetFE(int i) const override
Definition: pfespace.cpp:534
Class for parallel grid function.
Definition: pgridfunc.hpp:33
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:562
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:958
Class for parallel meshes.
Definition: pmesh.hpp:34
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition: pmesh.cpp:2028
void PrintAsOne(std::ostream &out=mfem::out, const std::string &comments="") const
Definition: pmesh.cpp:4968
Parallel non-linear operator on the true dofs.
Version of QuadraticFECollection with positive basis functions.
Definition: fe_coll.hpp:796
Base class for solvers.
Definition: operator.hpp:683
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:4839
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Adds the limiting term to the first integrator. Disables it for the rest.
Definition: tmop.cpp:4724
void AddTMOPIntegrator(TMOP_Integrator *ti)
Adds a new TMOP_Integrator to the combination.
Definition: tmop.hpp:2271
void AddGridFunctionForUpdate(GridFunction *gf)
Definition: tmop_amr.hpp:252
void AddFESpaceForUpdate(FiniteElementSpace *fes)
Definition: tmop_amr.hpp:253
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:207
void SetMinDetPtr(real_t *md_ptr)
Definition: tmop_tools.hpp:213
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: tmop_tools.hpp:286
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:1114
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:1132
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:1150
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1738
void SetExactActionFlag(bool flag_)
Flag to control if exact action of Integration is effected.
Definition: tmop.hpp:2241
void SetCoefficient(Coefficient &w1)
Sets a scaling Coefficient for the quality metric term of the integrator.
Definition: tmop.hpp:2045
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:4535
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop.hpp:2021
void IntegrateOverTarget(bool integ_over_target_)
Definition: tmop.hpp:2030
2D non-barrier metric without a type.
Definition: tmop.hpp:222
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:341
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:359
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:375
2D Shifted barrier form of shape metric (mu_2).
Definition: tmop.hpp:394
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:471
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:557
2D compound barrier Shape+Size (VS) metric (balanced).
Definition: tmop.hpp:572
2D compound barrier Shape+Size (VS) metric (balanced).
Definition: tmop.hpp:593
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:614
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:668
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:687
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:708
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:729
3D Shape (S) metric, untangling version of 303.
Definition: tmop.hpp:769
3D Size (V) metric.
Definition: tmop.hpp:790
3D Size (V) metric.
Definition: tmop.hpp:808
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:848
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:869
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:890
3D compound barrier Shape+Size (VS) metric (polyconvex, balanced).
Definition: tmop.hpp:911
3D compound barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:932
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:953
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:972
3D compound barrier Shape+Size (VS) metric (polyconvex, balanced).
Definition: tmop.hpp:994
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:1015
3D non-barrier Shape (S) metric.
Definition: tmop.hpp:1056
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P....
Definition: tmop.hpp:24
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1334
void SetVolumeScale(real_t vol_scale)
Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
Definition: tmop.hpp:1419
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:1413
TargetType
Target-matrix construction algorithms supported by this class.
Definition: tmop.hpp:1338
A general vector function coefficient.
Vector data type.
Definition: vector.hpp:80
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:816
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:494
real_t Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1212
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int dim
Definition: ex24.cpp:53
@ disc
Definition: ex25.cpp:151
int main()
real_t a
Definition: lissajous.cpp:41
real_t adapt_lim_fun(const Vector &x)
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
real_t weight_fun(const Vector &x)
real_t material_indicator_2d(const Vector &x)
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
void ConstructSizeGF(GridFunction &size)
real_t discrete_ori_2d(const Vector &x)
void discrete_aspr_3d(const Vector &x, Vector &v)
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:22
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
Definition: fem_extras.cpp:92
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:45
Geometry Geometries
Definition: fe.cpp:49
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
Definition: tmop_tools.cpp:893
float real_t
Definition: config.hpp:43
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:486
Settings for the output behavior of the IterativeSolver.
Definition: solvers.hpp:79
Helper struct to convert a C++ type to an MPI type.