MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pmesh-fitting.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// Boundary and Interface Fitting Miniapp
14// --------------------------------------------------------------
15//
16// This miniapp performs mesh optimization for controlling mesh quality and
17// aligning a selected set of nodes to boundary and/or interface of interest
18// defined using a level-set function. The mesh quality aspect is based on a
19// variational formulation of the Target-Matrix Optimization Paradigm (TMOP).
20// Boundary/interface alignment is weakly enforced using a penalization term
21// that moves a selected set of nodes towards the zero level set of a signed
22// smooth discrete function.
23//
24// See the following papers for more details:
25//
26// [1] "Adaptive Surface Fitting and Tangential Relaxation for High-Order Mesh Optimization" by
27// Knupp, Kolev, Mittal, Tomov.
28// [2] "High-Order Mesh Morphing for Boundary and Interface Fitting to Implicit Geometries" by
29// Barrera, Kolev, Mittal, Tomov.
30// [3] "The target-matrix optimization paradigm for high-order meshes" by
31// Dobrev, Knupp, Kolev, Mittal, Tomov.
32//
33// Compile with: make pmesh-fitting
34//
35// Sample runs:
36// Interface fitting:
37// mpirun -np 4 pmesh-fitting -o 3 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 5e4 -rtol 1e-5
38// mpirun -np 4 pmesh-fitting -m square01-tri.mesh -o 3 -rs 0 -mid 58 -tid 1 -ni 200 -vl 1 -sfc 1e4 -rtol 1e-5
39// Surface fitting with weight adaptation and termination based on fitting error:
40// mpirun -np 4 pmesh-fitting -o 2 -mid 2 -tid 1 -ni 100 -vl 2 -sfc 10 -rtol 1e-20 -st 0 -sfa 10.0 -sft 1e-5
41// Fitting to Fischer-Tropsch reactor like domain (requires GSLIB):
42// * mpirun -np 6 pmesh-fitting -m ../../data/inline-tri.mesh -o 2 -rs 4 -mid 2 -tid 1 -vl 2 -sfc 100 -rtol 1e-12 -ni 100 -li 40 -ae 1 -bnd -sbgmesh -slstype 2 -smtype 0 -sfa 10.0 -sft 1e-4 -amriter 5 -dist -mod-bndr-attr
43
44#include "mesh-fitting.hpp"
45
46using namespace mfem;
47using namespace std;
48
49int main (int argc, char *argv[])
50{
51 // 0. Initialize MPI and HYPRE.
52 Mpi::Init(argc, argv);
53 int myid = Mpi::WorldRank();
55
56 // 1. Set the method's default parameters.
57 const char *mesh_file = "square01.mesh";
58 int mesh_poly_deg = 1;
59 int rs_levels = 1;
60 int rp_levels = 0;
61 int metric_id = 2;
62 int target_id = 1;
63 real_t surface_fit_const = 100.0;
64 int quad_type = 1;
65 int quad_order = 8;
66 int solver_type = 0;
67 int solver_iter = 20;
68#ifdef MFEM_USE_SINGLE
69 real_t solver_rtol = 1e-4;
70#else
71 real_t solver_rtol = 1e-10;
72#endif
73 int solver_art_type = 0;
74 int lin_solver = 2;
75 int max_lin_iter = 100;
76 bool move_bnd = true;
77 bool visualization = true;
78 int verbosity_level = 0;
79 int adapt_eval = 0;
80 const char *devopt = "cpu";
81 real_t surface_fit_adapt = 0.0;
82 real_t surface_fit_threshold = -10;
83 bool adapt_marking = false;
84 bool surf_bg_mesh = false;
85 bool comp_dist = false;
86 int surf_ls_type = 1;
87 int marking_type = 0;
88 bool mod_bndr_attr = false;
89 bool material = false;
90 int mesh_node_ordering = 0;
91 int amr_iters = 0;
92
93 // 2. Parse command-line options.
94 OptionsParser args(argc, argv);
95 args.AddOption(&mesh_file, "-m", "--mesh",
96 "Mesh file to use.");
97 args.AddOption(&mesh_poly_deg, "-o", "--order",
98 "Polynomial degree of mesh finite element space.");
99 args.AddOption(&rs_levels, "-rs", "--refine-serial",
100 "Number of times to refine the mesh uniformly in serial.");
101 args.AddOption(&rp_levels, "-rp", "--refine-parallel",
102 "Number of times to refine the mesh uniformly in parallel.");
103 args.AddOption(&metric_id, "-mid", "--metric-id",
104 "Mesh optimization metric. See list in mesh-optimizer.");
105 args.AddOption(&target_id, "-tid", "--target-id",
106 "Target (ideal element) type:\n\t"
107 "1: Ideal shape, unit size\n\t"
108 "2: Ideal shape, equal size\n\t"
109 "3: Ideal shape, initial size\n\t"
110 "4: Given full analytic Jacobian (in physical space)\n\t"
111 "5: Ideal shape, given size (in physical space)");
112 args.AddOption(&surface_fit_const, "-sfc", "--surface-fit-const",
113 "Surface preservation constant.");
114 args.AddOption(&quad_type, "-qt", "--quad-type",
115 "Quadrature rule type:\n\t"
116 "1: Gauss-Lobatto\n\t"
117 "2: Gauss-Legendre\n\t"
118 "3: Closed uniform points");
119 args.AddOption(&quad_order, "-qo", "--quad_order",
120 "Order of the quadrature rule.");
121 args.AddOption(&solver_type, "-st", "--solver-type",
122 " Type of solver: (default) 0: Newton, 1: LBFGS");
123 args.AddOption(&solver_iter, "-ni", "--newton-iters",
124 "Maximum number of Newton iterations.");
125 args.AddOption(&solver_rtol, "-rtol", "--newton-rel-tolerance",
126 "Relative tolerance for the Newton solver.");
127 args.AddOption(&solver_art_type, "-art", "--adaptive-rel-tol",
128 "Type of adaptive relative linear solver tolerance:\n\t"
129 "0: None (default)\n\t"
130 "1: Eisenstat-Walker type 1\n\t"
131 "2: Eisenstat-Walker type 2");
132 args.AddOption(&lin_solver, "-ls", "--lin-solver",
133 "Linear solver:\n\t"
134 "0: l1-Jacobi\n\t"
135 "1: CG\n\t"
136 "2: MINRES\n\t"
137 "3: MINRES + Jacobi preconditioner\n\t"
138 "4: MINRES + l1-Jacobi preconditioner");
139 args.AddOption(&max_lin_iter, "-li", "--lin-iter",
140 "Maximum number of iterations in the linear solve.");
141 args.AddOption(&move_bnd, "-bnd", "--move-boundary", "-fix-bnd",
142 "--fix-boundary",
143 "Enable motion along horizontal and vertical boundaries.");
144 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
145 "--no-visualization",
146 "Enable or disable GLVis visualization.");
147 args.AddOption(&verbosity_level, "-vl", "--verbosity-level",
148 "Set the verbosity level - 0, 1, or 2.");
149 args.AddOption(&adapt_eval, "-ae", "--adaptivity-evaluator",
150 "0 - Advection based (DEFAULT), 1 - GSLIB.");
151 args.AddOption(&devopt, "-d", "--device",
152 "Device configuration string, see Device::Configure().");
153 args.AddOption(&surface_fit_adapt, "-sfa", "--adaptive-surface-fit",
154 "Enable or disable adaptive surface fitting.");
155 args.AddOption(&surface_fit_threshold, "-sft", "--surf-fit-threshold",
156 "Set threshold for surface fitting. TMOP solver will"
157 "terminate when max surface fitting error is below this limit");
158 args.AddOption(&adapt_marking, "-marking", "--adaptive-marking", "-no-amarking",
159 "--no-adaptive-marking",
160 "Enable or disable adaptive marking surface fitting.");
161 args.AddOption(&surf_bg_mesh, "-sbgmesh", "--surf-bg-mesh",
162 "-no-sbgmesh","--no-surf-bg-mesh",
163 "Use background mesh for surface fitting.");
164 args.AddOption(&comp_dist, "-dist", "--comp-dist",
165 "-no-dist","--no-comp-dist",
166 "Compute distance from 0 level set or not.");
167 args.AddOption(&surf_ls_type, "-slstype", "--surf-ls-type",
168 "1 - Circle (DEFAULT), 2 - Squircle, 3 - Butterfly.");
169 args.AddOption(&marking_type, "-smtype", "--surf-marking-type",
170 "1 - Interface (DEFAULT), 2 - Boundary attribute.");
171 args.AddOption(&mod_bndr_attr, "-mod-bndr-attr", "--modify-boundary-attribute",
172 "-fix-bndr-attr", "--fix-boundary-attribute",
173 "Change boundary attribute based on alignment with Cartesian axes.");
174 args.AddOption(&material, "-mat", "--mat",
175 "-no-mat","--no-mat", "Use default material attributes.");
176 args.AddOption(&mesh_node_ordering, "-mno", "--mesh_node_ordering",
177 "Ordering of mesh nodes."
178 "0 (default): byNodes, 1: byVDIM");
179 args.AddOption(&amr_iters, "-amriter", "--amr-iter",
180 "Number of amr iterations on background mesh");
181 args.Parse();
182 if (!args.Good())
183 {
184 if (myid == 0) { args.PrintUsage(cout); }
185 return 1;
186 }
187 if (myid == 0) { args.PrintOptions(cout); }
188
189 Device device(devopt);
190 if (myid == 0) { device.Print();}
191
192 // 3. Initialize and refine the starting mesh.
193 Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
194 for (int lev = 0; lev < rs_levels; lev++)
195 {
196 mesh->UniformRefinement();
197 }
198 const int dim = mesh->Dimension();
199
200 // Define level-set coefficient
201 FunctionCoefficient *ls_coeff = NULL;
202 if (surf_ls_type == 1) //Circle
203 {
205 }
206 else if (surf_ls_type == 2) // reactor
207 {
208 ls_coeff = new FunctionCoefficient(reactor);
209 }
210 else if (surf_ls_type == 6) // 3D shape
211 {
212 ls_coeff = new FunctionCoefficient(csg_cubecylsph);
213 }
214 else
215 {
216 MFEM_ABORT("Surface fitting level set type not implemented yet.")
217 }
218
219 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
220 delete mesh;
221 for (int lev = 0; lev < rp_levels; lev++) { pmesh->UniformRefinement(); }
222
223 // 4. Setup background mesh for surface fitting
224 ParMesh *pmesh_surf_fit_bg = NULL;
225 if (surf_bg_mesh)
226 {
227 Mesh *mesh_surf_fit_bg = NULL;
228 if (dim == 2)
229 {
230 mesh_surf_fit_bg =
232 }
233 else if (dim == 3)
234 {
235 mesh_surf_fit_bg =
237 }
238 mesh_surf_fit_bg->EnsureNCMesh();
239 pmesh_surf_fit_bg = new ParMesh(MPI_COMM_WORLD, *mesh_surf_fit_bg);
240 delete mesh_surf_fit_bg;
241 }
242
243 // 5. Define a finite element space on the mesh. Here we use vector finite
244 // elements which are tensor products of quadratic finite elements. The
245 // number of components in the vector finite element space is specified by
246 // the last parameter of the FiniteElementSpace constructor.
248 if (mesh_poly_deg <= 0)
249 {
250 fec = new QuadraticPosFECollection;
251 mesh_poly_deg = 2;
252 }
253 else { fec = new H1_FECollection(mesh_poly_deg, dim); }
254 ParFiniteElementSpace *pfespace =
255 new ParFiniteElementSpace(pmesh, fec, dim, mesh_node_ordering);
256
257 // 6. Make the mesh curved based on the above finite element space. This
258 // means that we define the mesh elements through a fespace-based
259 // transformation of the reference element.
260 pmesh->SetNodalFESpace(pfespace);
261
262 // 7. Get the mesh nodes (vertices and other degrees of freedom in the finite
263 // element space) as a finite element grid function in fespace. Note that
264 // changing x automatically changes the shapes of the mesh elements.
265 ParGridFunction x(pfespace);
266 pmesh->SetNodalGridFunction(&x);
267 x.SetTrueVector();
268
269 // 10. Save the starting (prior to the optimization) mesh to a file. This
270 // output can be viewed later using GLVis: "glvis -m perturbed -np
271 // num_mpi_tasks".
272 {
273 ostringstream mesh_name;
274 mesh_name << "perturbed.mesh";
275 ofstream mesh_ofs(mesh_name.str().c_str());
276 mesh_ofs.precision(8);
277 pmesh->PrintAsSerial(mesh_ofs);
278 }
279
280 // 11. Store the starting (prior to the optimization) positions.
281 ParGridFunction x0(pfespace);
282 x0 = x;
283
284 // 12. Form the integrator that uses the chosen metric and target.
285 TMOP_QualityMetric *metric = NULL;
286 switch (metric_id)
287 {
288 // T-metrics
289 case 2: metric = new TMOP_Metric_002; break;
290 case 58: metric = new TMOP_Metric_058; break;
291 case 80: metric = new TMOP_Metric_080(0.5); break;
292 case 303: metric = new TMOP_Metric_303; break;
293 case 328: metric = new TMOP_Metric_328; break;
294 default:
295 if (myid == 0) { cout << "Unknown metric_id: " << metric_id << endl; }
296 return 3;
297 }
298
299 if (metric_id < 300)
300 {
301 MFEM_VERIFY(dim == 2, "Incompatible metric for 3D meshes");
302 }
303 if (metric_id >= 300)
304 {
305 MFEM_VERIFY(dim == 3, "Incompatible metric for 2D meshes");
306 }
307
309 TargetConstructor *target_c = NULL;
310 switch (target_id)
311 {
312 case 1: target_t = TargetConstructor::IDEAL_SHAPE_UNIT_SIZE; break;
313 case 2: target_t = TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE; break;
314 case 3: target_t = TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE; break;
315 case 4: target_t = TargetConstructor::GIVEN_SHAPE_AND_SIZE; break;
316 default:
317 if (myid == 0) { cout << "Unknown target_id: " << target_id << endl; }
318 return 3;
319 }
320
321 if (target_c == NULL)
322 {
323 target_c = new TargetConstructor(target_t, MPI_COMM_WORLD);
324 }
325 target_c->SetNodes(x0);
326 TMOP_Integrator *tmop_integ = new TMOP_Integrator(metric, target_c);
327
328 // Setup the quadrature rules for the TMOP integrator.
329 IntegrationRules *irules = NULL;
330 switch (quad_type)
331 {
332 case 1: irules = &IntRulesLo; break;
333 case 2: irules = &IntRules; break;
334 case 3: irules = &IntRulesCU; break;
335 default:
336 if (myid == 0) { cout << "Unknown quad_type: " << quad_type << endl; }
337 return 3;
338 }
339 tmop_integ->SetIntegrationRules(*irules, quad_order);
340 if (myid == 0 && dim == 2)
341 {
342 cout << "Triangle quadrature points: "
343 << irules->Get(Geometry::TRIANGLE, quad_order).GetNPoints()
344 << "\nQuadrilateral quadrature points: "
345 << irules->Get(Geometry::SQUARE, quad_order).GetNPoints() << endl;
346 }
347 if (myid == 0 && dim == 3)
348 {
349 cout << "Tetrahedron quadrature points: "
350 << irules->Get(Geometry::TETRAHEDRON, quad_order).GetNPoints()
351 << "\nHexahedron quadrature points: "
352 << irules->Get(Geometry::CUBE, quad_order).GetNPoints()
353 << "\nPrism quadrature points: "
354 << irules->Get(Geometry::PRISM, quad_order).GetNPoints() << endl;
355 }
356
357 // Modify boundary attribute for surface node movement
358 // Sets attributes of a boundary element to 1/2/3 if it is parallel to x/y/z.
359 if (mod_bndr_attr)
360 {
362 pmesh->SetAttributes();
363 }
364 pmesh->ExchangeFaceNbrData();
365
366 // Surface fitting.
367 L2_FECollection mat_coll(0, dim);
368 H1_FECollection surf_fit_fec(mesh_poly_deg, dim);
369 ParFiniteElementSpace surf_fit_fes(pmesh, &surf_fit_fec);
370 ParFiniteElementSpace mat_fes(pmesh, &mat_coll);
371 ParGridFunction mat(&mat_fes);
372 ParGridFunction surf_fit_mat_gf(&surf_fit_fes);
373 ParGridFunction surf_fit_gf0(&surf_fit_fes);
374 Array<bool> surf_fit_marker(surf_fit_gf0.Size());
375 ConstantCoefficient surf_fit_coeff(surface_fit_const);
376 AdaptivityEvaluator *adapt_surface = NULL;
377 AdaptivityEvaluator *adapt_grad_surface = NULL;
378 AdaptivityEvaluator *adapt_hess_surface = NULL;
379
380 // Background mesh FECollection, FESpace, and GridFunction
381 H1_FECollection *surf_fit_bg_fec = NULL;
382 ParFiniteElementSpace *surf_fit_bg_fes = NULL;
383 ParGridFunction *surf_fit_bg_gf0 = NULL;
384 ParFiniteElementSpace *surf_fit_bg_grad_fes = NULL;
385 ParGridFunction *surf_fit_bg_grad = NULL;
386 ParFiniteElementSpace *surf_fit_bg_hess_fes = NULL;
387 ParGridFunction *surf_fit_bg_hess = NULL;
388
389 // If a background mesh is used, we interpolate the Gradient and Hessian
390 // from that mesh to the current mesh being optimized.
391 ParFiniteElementSpace *surf_fit_grad_fes = NULL;
392 ParGridFunction *surf_fit_grad = NULL;
393 ParFiniteElementSpace *surf_fit_hess_fes = NULL;
394 ParGridFunction *surf_fit_hess = NULL;
395
396 if (surf_bg_mesh)
397 {
398 pmesh_surf_fit_bg->SetCurvature(mesh_poly_deg);
399
400 Vector p_min(dim), p_max(dim);
401 pmesh->GetBoundingBox(p_min, p_max);
402 GridFunction &x_bg = *pmesh_surf_fit_bg->GetNodes();
403 const int num_nodes = x_bg.Size() / dim;
404 for (int i = 0; i < num_nodes; i++)
405 {
406 for (int d = 0; d < dim; d++)
407 {
408 real_t length_d = p_max(d) - p_min(d),
409 extra_d = 0.2 * length_d;
410 x_bg(i + d*num_nodes) = p_min(d) - extra_d +
411 x_bg(i + d*num_nodes) * (length_d + 2*extra_d);
412 }
413 }
414 surf_fit_bg_fec = new H1_FECollection(mesh_poly_deg+1, dim);
415 surf_fit_bg_fes = new ParFiniteElementSpace(pmesh_surf_fit_bg, surf_fit_bg_fec);
416 surf_fit_bg_gf0 = new ParGridFunction(surf_fit_bg_fes);
417 }
418
419 Array<int> vdofs;
420 if (surface_fit_const > 0.0)
421 {
422 surf_fit_gf0.ProjectCoefficient(*ls_coeff);
423 if (surf_bg_mesh)
424 {
425 OptimizeMeshWithAMRAroundZeroLevelSet(*pmesh_surf_fit_bg, *ls_coeff,
426 amr_iters, *surf_fit_bg_gf0);
427 pmesh_surf_fit_bg->Rebalance();
428 surf_fit_bg_fes->Update();
429 surf_fit_bg_gf0->Update();
430
431 if (comp_dist)
432 {
433 ComputeScalarDistanceFromLevelSet(*pmesh_surf_fit_bg, *ls_coeff,
434 *surf_fit_bg_gf0);
435 }
436 else { surf_fit_bg_gf0->ProjectCoefficient(*ls_coeff); }
437
438 surf_fit_bg_grad_fes =
439 new ParFiniteElementSpace(pmesh_surf_fit_bg, surf_fit_bg_fec, dim);
440 surf_fit_bg_grad = new ParGridFunction(surf_fit_bg_grad_fes);
441
442 surf_fit_grad_fes =
443 new ParFiniteElementSpace(pmesh, &surf_fit_fec, dim);
444 surf_fit_grad = new ParGridFunction(surf_fit_grad_fes);
445
446 surf_fit_bg_hess_fes =
447 new ParFiniteElementSpace(pmesh_surf_fit_bg, surf_fit_bg_fec, dim * dim);
448 surf_fit_bg_hess = new ParGridFunction(surf_fit_bg_hess_fes);
449
450 surf_fit_hess_fes =
451 new ParFiniteElementSpace(pmesh, &surf_fit_fec, dim * dim);
452 surf_fit_hess = new ParGridFunction(surf_fit_hess_fes);
453
454 //Setup gradient of the background mesh
455 const int size_bg = surf_fit_bg_gf0->Size();
456 for (int d = 0; d < pmesh_surf_fit_bg->Dimension(); d++)
457 {
458 ParGridFunction surf_fit_bg_grad_comp(
459 surf_fit_bg_fes, surf_fit_bg_grad->GetData() + d * size_bg);
460 surf_fit_bg_gf0->GetDerivative(1, d, surf_fit_bg_grad_comp);
461 }
462
463 //Setup Hessian on background mesh
464 int id = 0;
465 for (int d = 0; d < pmesh_surf_fit_bg->Dimension(); d++)
466 {
467 for (int idir = 0; idir < pmesh_surf_fit_bg->Dimension(); idir++)
468 {
469 ParGridFunction surf_fit_bg_grad_comp(
470 surf_fit_bg_fes, surf_fit_bg_grad->GetData() + d * size_bg);
471 ParGridFunction surf_fit_bg_hess_comp(
472 surf_fit_bg_fes, surf_fit_bg_hess->GetData()+ id * size_bg);
473 surf_fit_bg_grad_comp.GetDerivative(1, idir,
474 surf_fit_bg_hess_comp);
475 id++;
476 }
477 }
478 }
479 else // !surf_bg_mesh
480 {
481 if (comp_dist)
482 {
483 ComputeScalarDistanceFromLevelSet(*pmesh, *ls_coeff, surf_fit_gf0);
484 }
485 }
486
487 // Set material gridfunction
488 for (int i = 0; i < pmesh->GetNE(); i++)
489 {
490 if (material)
491 {
492 mat(i) = pmesh->GetAttribute(i)-1;
493 }
494 else
495 {
496 mat(i) = material_id(i, surf_fit_gf0);
497 pmesh->SetAttribute(i, mat(i) + 1);
498 }
499 }
500
501 // Adapt attributes for marking such that if all but 1 face of an element
502 // are marked, the element attribute is switched.
503 if (adapt_marking)
504 {
505 ModifyAttributeForMarkingDOFS(pmesh, mat, 0);
506 ModifyAttributeForMarkingDOFS(pmesh, mat, 1);
507 }
508
509 GridFunctionCoefficient coeff_mat(&mat);
510 surf_fit_mat_gf.ProjectDiscCoefficient(coeff_mat,
512 surf_fit_mat_gf.SetTrueVector();
513 surf_fit_mat_gf.SetFromTrueVector();
514
515 // Set DOFs for fitting
516 // Strategy 1: Choose face between elements of different attributes.
517 if (marking_type == 0)
518 {
520 const Vector &FaceNbrData = mat.FaceNbrData();
521 for (int j = 0; j < surf_fit_marker.Size(); j++)
522 {
523 surf_fit_marker[j] = false;
524 }
525 surf_fit_mat_gf = 0.0;
526
527 Array<int> dof_list;
528 Array<int> dofs;
529 for (int i = 0; i < pmesh->GetNumFaces(); i++)
530 {
531 auto tr = pmesh->GetInteriorFaceTransformations(i);
532 if (tr != NULL)
533 {
534 int mat1 = mat(tr->Elem1No);
535 int mat2 = mat(tr->Elem2No);
536 if (mat1 != mat2)
537 {
538 surf_fit_gf0.ParFESpace()->GetFaceDofs(i, dofs);
539 dof_list.Append(dofs);
540 }
541 }
542 }
543 for (int i = 0; i < pmesh->GetNSharedFaces(); i++)
544 {
545 auto tr = pmesh->GetSharedFaceTransformations(i);
546 if (tr != NULL)
547 {
548 int faceno = pmesh->GetSharedFace(i);
549 int mat1 = mat(tr->Elem1No);
550 int mat2 = FaceNbrData(tr->Elem2No-pmesh->GetNE());
551 if (mat1 != mat2)
552 {
553 surf_fit_gf0.ParFESpace()->GetFaceDofs(faceno, dofs);
554 dof_list.Append(dofs);
555 }
556 }
557 }
558 for (int i = 0; i < dof_list.Size(); i++)
559 {
560 surf_fit_marker[dof_list[i]] = true;
561 surf_fit_mat_gf(dof_list[i]) = 1.0;
562 }
563 }
564 // Strategy 2: Mark all boundaries with attribute marking_type
565 else if (marking_type > 0)
566 {
567 for (int i = 0; i < pmesh->GetNBE(); i++)
568 {
569 const int attr = pmesh->GetBdrElement(i)->GetAttribute();
570 if (attr == marking_type)
571 {
572 surf_fit_fes.GetBdrElementVDofs(i, vdofs);
573 for (int j = 0; j < vdofs.Size(); j++)
574 {
575 surf_fit_marker[vdofs[j]] = true;
576 surf_fit_mat_gf(vdofs[j]) = 1.0;
577 }
578 }
579 }
580 }
581
582 // Set AdaptivityEvaluators for transferring information from initial
583 // mesh to current mesh as it moves during adaptivity.
584 if (adapt_eval == 0)
585 {
586 adapt_surface = new AdvectorCG;
587 MFEM_VERIFY(!surf_bg_mesh, "Background meshes require GSLIB.");
588 }
589 else if (adapt_eval == 1)
590 {
591#ifdef MFEM_USE_GSLIB
592 adapt_surface = new InterpolatorFP;
593 if (surf_bg_mesh)
594 {
595 adapt_grad_surface = new InterpolatorFP;
596 adapt_hess_surface = new InterpolatorFP;
597 }
598#else
599 MFEM_ABORT("MFEM is not built with GSLIB support!");
600#endif
601 }
602 else { MFEM_ABORT("Bad interpolation option."); }
603
604 if (!surf_bg_mesh)
605 {
606 tmop_integ->EnableSurfaceFitting(surf_fit_gf0, surf_fit_marker,
607 surf_fit_coeff, *adapt_surface);
608 }
609 else
610 {
612 *surf_fit_bg_gf0, surf_fit_gf0,
613 surf_fit_marker, surf_fit_coeff, *adapt_surface,
614 *surf_fit_bg_grad, *surf_fit_grad, *adapt_grad_surface,
615 *surf_fit_bg_hess, *surf_fit_hess, *adapt_hess_surface);
616 }
617
618 if (visualization)
619 {
620 socketstream vis1, vis2, vis3, vis4, vis5;
621 common::VisualizeField(vis1, "localhost", 19916, surf_fit_gf0,
622 "Level Set", 0, 0, 300, 300);
623 common::VisualizeField(vis2, "localhost", 19916, mat,
624 "Materials", 300, 0, 300, 300);
625 common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
626 "Surface DOFs", 600, 0, 300, 300);
627 if (surf_bg_mesh)
628 {
629 common::VisualizeField(vis4, "localhost", 19916, *surf_fit_bg_gf0,
630 "Level Set - Background",
631 0, 400, 300, 300);
632 }
633 }
634 }
635 pmesh->SetAttributes();
636
637 // 13. Setup the final NonlinearForm (which defines the integral of interest,
638 // its first and second derivatives). Here we can use a combination of
639 // metrics, i.e., optimize the sum of two integrals, where both are
640 // scaled by used-defined space-dependent weights. Note that there are
641 // no command-line options for the weights and the type of the second
642 // metric; one should update those in the code.
643 ParNonlinearForm a(pfespace);
644 ConstantCoefficient *metric_coeff1 = NULL;
645 a.AddDomainIntegrator(tmop_integ);
646
647 // Compute the minimum det(J) of the starting mesh.
648 real_t min_detJ = infinity();
649 const int NE = pmesh->GetNE();
650 for (int i = 0; i < NE; i++)
651 {
652 const IntegrationRule &ir =
653 irules->Get(pfespace->GetFE(i)->GetGeomType(), quad_order);
655 for (int j = 0; j < ir.GetNPoints(); j++)
656 {
657 transf->SetIntPoint(&ir.IntPoint(j));
658 min_detJ = min(min_detJ, transf->Jacobian().Det());
659 }
660 }
661 MPI_Allreduce(MPI_IN_PLACE, &min_detJ, 1,
662 MPITypeMap<real_t>::mpi_type, MPI_MIN, MPI_COMM_WORLD);
663 if (myid == 0)
664 { cout << "Minimum det(J) of the original mesh is " << min_detJ << endl; }
665
666 MFEM_VERIFY(min_detJ > 0, "The input mesh is inverted, use mesh-optimizer.");
667
668 const real_t init_energy = a.GetParGridFunctionEnergy(x);
669 real_t init_metric_energy = init_energy;
670 if (surface_fit_const > 0.0)
671 {
672 surf_fit_coeff.constant = 0.0;
673 init_metric_energy = a.GetParGridFunctionEnergy(x);
674 surf_fit_coeff.constant = surface_fit_const;
675 }
676
677 // 14. Fix all boundary nodes, or fix only a given component depending on the
678 // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
679 // fixed x/y/z components of the node. Attribute dim+1 corresponds to
680 // an entirely fixed node.
681 if (move_bnd == false)
682 {
683 Array<int> ess_bdr(pmesh->bdr_attributes.Max());
684 ess_bdr = 1;
685 if (marking_type > 0)
686 {
687 ess_bdr[marking_type-1] = 0;
688 }
689 a.SetEssentialBC(ess_bdr);
690 }
691 else
692 {
693 int n = 0;
694 for (int i = 0; i < pmesh->GetNBE(); i++)
695 {
696 const int nd = pfespace->GetBE(i)->GetDof();
697 const int attr = pmesh->GetBdrElement(i)->GetAttribute();
698 MFEM_VERIFY(!(dim == 2 && attr == 3),
699 "Boundary attribute 3 must be used only for 3D meshes. "
700 "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
701 "components, rest for free nodes), or use -fix-bnd.");
702 if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
703 if (attr == 4) { n += nd * dim; }
704 }
705 Array<int> ess_vdofs(n);
706 n = 0;
707 for (int i = 0; i < pmesh->GetNBE(); i++)
708 {
709 const int nd = pfespace->GetBE(i)->GetDof();
710 const int attr = pmesh->GetBdrElement(i)->GetAttribute();
711 pfespace->GetBdrElementVDofs(i, vdofs);
712 if (attr == 1) // Fix x components.
713 {
714 for (int j = 0; j < nd; j++)
715 { ess_vdofs[n++] = vdofs[j]; }
716 }
717 else if (attr == 2) // Fix y components.
718 {
719 for (int j = 0; j < nd; j++)
720 { ess_vdofs[n++] = vdofs[j+nd]; }
721 }
722 else if (attr == 3) // Fix z components.
723 {
724 for (int j = 0; j < nd; j++)
725 { ess_vdofs[n++] = vdofs[j+2*nd]; }
726 }
727 else if (attr == 4) // Fix all components.
728 {
729 for (int j = 0; j < vdofs.Size(); j++)
730 { ess_vdofs[n++] = vdofs[j]; }
731 }
732 }
733 a.SetEssentialVDofs(ess_vdofs);
734 }
735
736 // 15. As we use the Newton method to solve the resulting nonlinear system,
737 // here we setup the linear solver for the system's Jacobian.
738 Solver *S = NULL, *S_prec = NULL;
739#ifdef MFEM_USE_SINGLE
740 const real_t linsol_rtol = 1e-5;
741#else
742 const real_t linsol_rtol = 1e-12;
743#endif
744 if (lin_solver == 0)
745 {
746 S = new DSmoother(1, 1.0, max_lin_iter);
747 }
748 else if (lin_solver == 1)
749 {
750 CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
751 cg->SetMaxIter(max_lin_iter);
752 cg->SetRelTol(linsol_rtol);
753 cg->SetAbsTol(0.0);
754 cg->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
755 S = cg;
756 }
757 else
758 {
759 MINRESSolver *minres = new MINRESSolver(MPI_COMM_WORLD);
760 minres->SetMaxIter(max_lin_iter);
761 minres->SetRelTol(linsol_rtol);
762 minres->SetAbsTol(0.0);
763 if (verbosity_level > 2) { minres->SetPrintLevel(1); }
764 else { minres->SetPrintLevel(verbosity_level == 2 ? 3 : -1); }
765 if (lin_solver == 3 || lin_solver == 4)
766 {
767 auto hs = new HypreSmoother;
768 hs->SetType((lin_solver == 3) ? HypreSmoother::Jacobi
769 /* */ : HypreSmoother::l1Jacobi, 1);
770 hs->SetPositiveDiagonal(true);
771 S_prec = hs;
772 minres->SetPreconditioner(*S_prec);
773 }
774 S = minres;
775 }
776
777 // Perform the nonlinear optimization.
778 const IntegrationRule &ir =
779 irules->Get(pfespace->GetFE(0)->GetGeomType(), quad_order);
780 TMOPNewtonSolver solver(pfespace->GetComm(), ir, solver_type);
781 if (surface_fit_adapt > 0.0)
782 {
783 solver.SetAdaptiveSurfaceFittingScalingFactor(surface_fit_adapt);
784 }
785 if (surface_fit_threshold > 0)
786 {
787 solver.SetTerminationWithMaxSurfaceFittingError(surface_fit_threshold);
788 }
789 // Provide all integration rules in case of a mixed mesh.
790 solver.SetIntegrationRules(*irules, quad_order);
791 if (solver_type == 0)
792 {
793 // Specify linear solver when we use a Newton-based solver.
794 solver.SetPreconditioner(*S);
795 }
796 solver.SetMaxIter(solver_iter);
797 solver.SetRelTol(solver_rtol);
798 solver.SetAbsTol(0.0);
799 solver.SetMinimumDeterminantThreshold(0.001*min_detJ);
800 if (solver_art_type > 0)
801 {
802 solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
803 }
804 solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
805 solver.SetOperator(a);
806 Vector b(0);
807 solver.Mult(b, x.GetTrueVector());
809
810 // 16. Save the optimized mesh to a file. This output can be viewed later
811 // using GLVis: "glvis -m optimized -np num_mpi_tasks".
812 {
813 ostringstream mesh_name;
814 mesh_name << "optimized.mesh";
815 ofstream mesh_ofs(mesh_name.str().c_str());
816 mesh_ofs.precision(8);
817 pmesh->PrintAsSerial(mesh_ofs);
818 }
819
820 // Compute the final energy of the functional.
821 const real_t fin_energy = a.GetParGridFunctionEnergy(x);
822 real_t fin_metric_energy = fin_energy;
823 if (surface_fit_const > 0.0)
824 {
825 surf_fit_coeff.constant = 0.0;
826 fin_metric_energy = a.GetParGridFunctionEnergy(x);
827 surf_fit_coeff.constant = surface_fit_const;
828 }
829
830 if (myid == 0)
831 {
832 std::cout << std::scientific << std::setprecision(4);
833 cout << "Initial strain energy: " << init_energy
834 << " = metrics: " << init_metric_energy
835 << " + extra terms: " << init_energy - init_metric_energy << endl;
836 cout << " Final strain energy: " << fin_energy
837 << " = metrics: " << fin_metric_energy
838 << " + extra terms: " << fin_energy - fin_metric_energy << endl;
839 cout << "The strain energy decreased by: "
840 << (init_energy - fin_energy) * 100.0 / init_energy << " %." << endl;
841 }
842
843 if (surface_fit_const > 0.0)
844 {
845 if (visualization)
846 {
847 socketstream vis2, vis3;
848 common::VisualizeField(vis2, "localhost", 19916, mat,
849 "Materials", 300, 400, 300, 300);
850 common::VisualizeField(vis3, "localhost", 19916, surf_fit_mat_gf,
851 "Surface DOFs", 600, 400, 300, 300);
852 }
853 real_t err_avg, err_max;
854 tmop_integ->GetSurfaceFittingErrors(x, err_avg, err_max);
855 if (myid == 0)
856 {
857 std::cout << "Avg fitting error: " << err_avg << std::endl
858 << "Max fitting error: " << err_max << std::endl;
859 }
860 }
861
862 // 18. Visualize the mesh displacement.
863 if (visualization)
864 {
865 x0 -= x;
866 socketstream vis;
867 common::VisualizeField(vis, "localhost", 19916, x0,
868 "Displacements", 900, 400, 300, 300, "jRmclA");
869 }
870
871 delete S;
872 delete S_prec;
873 delete metric_coeff1;
874 delete adapt_surface;
875 delete adapt_grad_surface;
876 delete adapt_hess_surface;
877 delete ls_coeff;
878 delete surf_fit_hess;
879 delete surf_fit_hess_fes;
880 delete surf_fit_bg_hess;
881 delete surf_fit_bg_hess_fes;
882 delete surf_fit_grad;
883 delete surf_fit_grad_fes;
884 delete surf_fit_bg_grad;
885 delete surf_fit_bg_grad_fes;
886 delete surf_fit_bg_gf0;
887 delete surf_fit_bg_fes;
888 delete surf_fit_bg_fec;
889 delete target_c;
890 delete metric;
891 delete pfespace;
892 delete fec;
893 delete pmesh_surf_fit_bg;
894 delete pmesh;
895
896 return 0;
897}
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
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition: array.hpp:769
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.
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
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
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
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.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
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
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:130
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 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
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:787
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 GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:6250
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1333
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1336
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
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
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, real_t sx=1.0, real_t sy=1.0, real_t sz=1.0, bool sfc_ordering=true)
Creates a mesh for the parallelepiped [0,sx]x[0,sy]x[0,sz], divided into nx*ny*nz hexahedra if type =...
Definition: mesh.cpp:4253
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:6200
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8973
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1229
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
Definition: mesh.cpp:1079
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition: mesh.cpp:10626
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, real_t sx=1.0, real_t sy=1.0, bool sfc_ordering=true)
Creates mesh for the rectangle [0,sx]x[0,sy], divided into nx*ny quadrilaterals if type = QUADRILATER...
Definition: mesh.cpp:4243
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10970
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
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
int GetFaceDofs(int i, Array< int > &dofs, int variant=0) const override
Definition: pfespace.cpp:518
void Update(bool want_transform=true) override
Definition: pfespace.cpp:3414
const FiniteElement * GetFE(int i) const override
Definition: pfespace.cpp:534
Class for parallel grid function.
Definition: pgridfunc.hpp:33
void ProjectDiscCoefficient(VectorCoefficient &coeff) override
Project a discontinuous vector coefficient as a grid function on a continuous finite element space....
Definition: pgridfunc.cpp:584
Vector & FaceNbrData()
Definition: pgridfunc.hpp:208
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
void GetDerivative(int comp, int der_comp, ParGridFunction &der)
Parallel version of GridFunction::GetDerivative(); see its documentation.
Definition: pgridfunc.cpp:522
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 Update() override
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:90
Class for parallel meshes.
Definition: pmesh.hpp:34
void ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
Definition: pmesh.cpp:2159
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:3153
void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1) override
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: pmesh.cpp:2007
void SetAttributes() override
Determine the sets of unique attribute values in domain and boundary elements.
Definition: pmesh.cpp:1589
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition: pmesh.cpp:2028
void Rebalance()
Definition: pmesh.cpp:4009
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
Definition: pmesh.cpp:3172
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Get the FaceElementTransformations for the given shared face (edge 2D) using the shared face index sf...
Definition: pmesh.cpp:2923
void PrintAsSerial(std::ostream &out=mfem::out, const std::string &comments="") const
Definition: pmesh.cpp:5281
void GetBoundingBox(Vector &p_min, Vector &p_max, int ref=2)
Definition: pmesh.cpp:6154
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 SetAdaptiveSurfaceFittingScalingFactor(real_t factor)
Definition: tmop_tools.hpp:240
void SetMinimumDeterminantThreshold(real_t threshold)
Definition: tmop_tools.hpp:256
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop_tools.hpp:207
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: tmop_tools.hpp:261
void SetTerminationWithMaxSurfaceFittingError(real_t max_error)
Definition: tmop_tools.hpp:252
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: tmop_tools.hpp:286
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1738
void EnableSurfaceFittingFromSource(const ParGridFunction &s_bg, ParGridFunction &s0, const Array< bool > &smarker, Coefficient &coeff, AdaptivityEvaluator &ae, const ParGridFunction &s_bg_grad, ParGridFunction &s0_grad, AdaptivityEvaluator &age, const ParGridFunction &s_bg_hess, ParGridFunction &s0_hess, AdaptivityEvaluator &ahe)
Fitting of certain DOFs in the current mesh to the zero level set of a function defined on another (f...
Definition: tmop.cpp:3010
void GetSurfaceFittingErrors(const Vector &pos, real_t &err_avg, real_t &err_max)
Definition: tmop.cpp:3079
void SetIntegrationRules(IntegrationRules &irules, int order)
Prescribe a set of integration rules; relevant for mixed meshes.
Definition: tmop.hpp:2021
void EnableSurfaceFitting(const GridFunction &s0, const Array< bool > &smarker, Coefficient &coeff, AdaptivityEvaluator &ae)
Fitting of certain DOFs to the zero level set of a function.
Definition: tmop.cpp:2943
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:471
3D barrier Shape (S) metric, well-posed (polyconvex & invex).
Definition: tmop.hpp:708
3D compound barrier Shape+Size (VS) metric (polyconvex, balanced).
Definition: tmop.hpp:911
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 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
Vector data type.
Definition: vector.hpp:80
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:227
int dim
Definition: ex24.cpp:53
int main()
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
real_t reactor(const Vector &x)
real_t csg_cubecylsph(const Vector &x)
void ModifyBoundaryAttributesForNodeMovement(ParMesh *pmesh, ParGridFunction &x)
void ModifyAttributeForMarkingDOFS(ParMesh *pmesh, ParGridFunction &mat, int attr_to_switch)
void OptimizeMeshWithAMRAroundZeroLevelSet(ParMesh &pmesh, FunctionCoefficient &ls_coeff, int amr_iter, ParGridFunction &distance_s, const int quad_order=5, Array< ParGridFunction * > *pgf_to_update=NULL)
void ComputeScalarDistanceFromLevelSet(ParMesh &pmesh, FunctionCoefficient &ls_coeff, ParGridFunction &distance_s, const int nDiffuse=2, const int pLapOrder=5, const int pLapNewton=50)
real_t circle_level_set(const Vector &x)
int material_id(int el_id, const GridFunction &g)
IntegrationRules IntRulesCU(0, Quadrature1D::ClosedUniform)
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
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
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
int material(Vector &x, Vector &xmin, Vector &xmax)
Definition: shaper.cpp:53
Helper struct to convert a C++ type to an MPI type.