MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
tmop_tools.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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#include "tmop_tools.hpp"
13#include "nonlinearform.hpp"
14#include "pnonlinearform.hpp"
16
17namespace mfem
18{
19
20using namespace mfem;
21
22void AdvectorCG::SetInitialField(const Vector &init_nodes,
23 const Vector &init_field)
24{
25 nodes0 = init_nodes;
26 field0 = init_field;
27}
28
29void AdvectorCG::ComputeAtNewPosition(const Vector &new_mesh_nodes,
30 Vector &new_field,
31 int nodes_ordering)
32{
33 MFEM_VERIFY(nodes0.Size() == new_mesh_nodes.Size(),
34 "AdvectorCG assumes fixed mesh topology!");
35
37#ifdef MFEM_USE_MPI
38 if (pfes) { space = pfes; }
39#endif
40 int fes_ordering = space->GetOrdering(),
41 ncomp = space->GetVDim();
42
43 const int dof_cnt = field0.Size() / ncomp;
44
45 new_field = field0;
46 Vector new_field_temp;
47 for (int i = 0; i < ncomp; i++)
48 {
49 if (fes_ordering == Ordering::byNODES)
50 {
51 new_field_temp.MakeRef(new_field, i*dof_cnt, dof_cnt);
52 }
53 else
54 {
55 new_field_temp.SetSize(dof_cnt);
56 for (int j = 0; j < dof_cnt; j++)
57 {
58 new_field_temp(j) = new_field(i + j*ncomp);
59 }
60 }
61 ComputeAtNewPositionScalar(new_mesh_nodes, new_field_temp);
62 if (fes_ordering == Ordering::byVDIM)
63 {
64 for (int j = 0; j < dof_cnt; j++)
65 {
66 new_field(i + j*ncomp) = new_field_temp(j);
67 }
68 }
69 }
70
71 field0 = new_field;
72 nodes0 = new_mesh_nodes;
73}
74
75void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_mesh_nodes,
76 Vector &new_field)
77{
78 Mesh *m = mesh;
79#ifdef MFEM_USE_MPI
80 if (pmesh) { m = pmesh; }
81#endif
82
83 MFEM_VERIFY(m != NULL, "No mesh has been given to the AdaptivityEvaluator.");
84
85 // This will be used to move the positions.
86 GridFunction *mesh_nodes = m->GetNodes();
87 *mesh_nodes = nodes0;
88 real_t minv = new_field.Min(), maxv = new_field.Max();
89
90 // Velocity of the positions.
91 GridFunction u(mesh_nodes->FESpace());
92 subtract(new_mesh_nodes, nodes0, u);
93
94 // Define a scalar FE space for the solution, and the advection operator.
95 TimeDependentOperator *oper = NULL;
96 FiniteElementSpace *fess = NULL;
97#ifdef MFEM_USE_MPI
98 ParFiniteElementSpace *pfess = NULL;
99#endif
100 if (fes)
101 {
102 fess = new FiniteElementSpace(fes->GetMesh(), fes->FEColl(), 1);
103 oper = new SerialAdvectorCGOper(nodes0, u, *fess, al);
104 }
105#ifdef MFEM_USE_MPI
106 else if (pfes)
107 {
108 pfess = new ParFiniteElementSpace(pfes->GetParMesh(), pfes->FEColl(), 1);
109 oper = new ParAdvectorCGOper(nodes0, u, *pfess, al, opt_mt);
110 }
111#endif
112 MFEM_VERIFY(oper != NULL,
113 "No FE space has been given to the AdaptivityEvaluator.");
114 ode_solver.Init(*oper);
115
116 // Compute some time step [mesh_size / speed].
117 real_t h_min = std::numeric_limits<real_t>::infinity();
118 for (int i = 0; i < m->GetNE(); i++)
119 {
120 h_min = std::min(h_min, m->GetElementSize(i));
121 }
122 real_t v_max = 0.0;
123 const int s = u.Size()/m->Dimension();
124
125 u.HostReadWrite();
126 for (int i = 0; i < s; i++)
127 {
128 real_t vel = 0.;
129 for (int j = 0; j < m->Dimension(); j++)
130 {
131 vel += u(i+j*s)*u(i+j*s);
132 }
133 v_max = std::max(v_max, vel);
134 }
135
136#ifdef MFEM_USE_MPI
137 if (pfes)
138 {
139 real_t v_loc = v_max, h_loc = h_min;
140 MPI_Allreduce(&v_loc, &v_max, 1, MPITypeMap<real_t>::mpi_type, MPI_MAX,
141 pfes->GetComm());
142 MPI_Allreduce(&h_loc, &h_min, 1, MPITypeMap<real_t>::mpi_type, MPI_MIN,
143 pfes->GetComm());
144 }
145#endif
146
147 if (v_max == 0.0) // No need to change the field.
148 {
149 delete oper;
150 delete fess;
151#ifdef MFEM_USE_MPI
152 delete pfess;
153#endif
154 return;
155 }
156
157 v_max = std::sqrt(v_max);
158 real_t dt = dt_scale * h_min / v_max;
159
160 real_t t = 0.0;
161 bool last_step = false;
162 while (!last_step)
163 {
164 if (t + dt >= 1.0)
165 {
166 dt = 1.0 - t;
167 last_step = true;
168 }
169 ode_solver.Step(new_field, t, dt);
170 }
171
172 real_t glob_minv = minv,
173 glob_maxv = maxv;
174#ifdef MFEM_USE_MPI
175 if (pfes)
176 {
177 MPI_Allreduce(&minv, &glob_minv, 1, MPITypeMap<real_t>::mpi_type, MPI_MIN,
178 pfes->GetComm());
179 MPI_Allreduce(&maxv, &glob_maxv, 1, MPITypeMap<real_t>::mpi_type, MPI_MAX,
180 pfes->GetComm());
181 }
182#endif
183
184 // Trim the overshoots and undershoots.
185 new_field.HostReadWrite();
186 for (int i = 0; i < new_field.Size(); i++)
187 {
188 if (new_field(i) < glob_minv) { new_field(i) = glob_minv; }
189 if (new_field(i) > glob_maxv) { new_field(i) = glob_maxv; }
190 }
191
192 delete oper;
193 delete fess;
194#ifdef MFEM_USE_MPI
195 delete pfess;
196#endif
197}
198
202 AssemblyLevel al)
203 : TimeDependentOperator(fes.GetVSize()),
204 x0(x_start), x_now(*fes.GetMesh()->GetNodes()),
205 u(vel), u_coeff(&u), M(&fes), K(&fes), al(al)
206{
208 K.AddDomainIntegrator(Kinteg);
210 K.Assemble(0);
211 K.Finalize(0);
212
213 MassIntegrator *Minteg = new MassIntegrator;
214 M.AddDomainIntegrator(Minteg);
216 M.Assemble(0);
217 M.Finalize(0);
218}
219
220void SerialAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
221{
222 // Move the mesh.
223 const real_t t = GetTime();
224 add(x0, t, u, x_now);
226
227 // Assemble on the new mesh.
228 K.BilinearForm::operator=(0.0);
229 K.Assemble();
230 Vector rhs(K.Size());
231 K.Mult(ind, rhs);
232 M.BilinearForm::operator=(0.0);
233 M.Assemble();
234
235#ifdef MFEM_USE_SINGLE
236 const real_t rtol = 1e-4;
237#else
238 const real_t rtol = 1e-12;
239#endif
240
241 // Solve.
242 di_dt = 0.0;
244 {
245 // Solve mat-free on the tdofs as the JacobiSmoother operates on tdofs.
246 OperatorPtr A;
247 Vector B, X;
248 Array<int> ess_tdof_list;
249 M.FormLinearSystem(ess_tdof_list, di_dt, rhs, A, X, B);
250 OperatorJacobiSmoother S(M, ess_tdof_list);
251 PCG(*A, S, B, X, 0, 100, rtol, 0.0);
252 M.RecoverFEMSolution(X, rhs, di_dt);
253 }
254 else
255 {
256 // Solve the SpMat directly on the ldofs.
257 DSmoother S(M.SpMat());
258 PCG(M.SpMat(), S, rhs, di_dt, 0, 100, rtol, 0.0);
259 }
260}
261
262#ifdef MFEM_USE_MPI
266 AssemblyLevel al,
267 MemoryType mt)
268 : TimeDependentOperator(pfes.GetVSize()),
269 x0(x_start), x_now(*pfes.GetMesh()->GetNodes()),
270 u(vel), u_coeff(&u), M(&pfes), K(&pfes), al(al)
271{
274 {
275 Kinteg->SetPAMemoryType(mt);
276 }
277 K.AddDomainIntegrator(Kinteg);
279 K.Assemble(0);
280 K.Finalize(0);
281
282 MassIntegrator *Minteg = new MassIntegrator;
284 {
285 Minteg->SetPAMemoryType(mt);
286 }
287 M.AddDomainIntegrator(Minteg);
289 M.Assemble(0);
290 M.Finalize(0);
291}
292
293void ParAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
294{
295 // Move the mesh.
296 const real_t t = GetTime();
297 add(x0, t, u, x_now);
299
300 // Assemble on the new mesh.
301 K.BilinearForm::operator=(0.0);
302 K.Assemble();
304 K.Mult(ind, rhs);
305 M.BilinearForm::operator=(0.0);
306 M.Assemble();
307
308 HypreParVector *RHS = rhs.ParallelAssemble();
310 X = 0.0;
311
312 OperatorHandle Mop;
313 Solver *prec = nullptr;
314 Array<int> ess_tdof_list;
316 {
317 M.FormSystemMatrix(ess_tdof_list, Mop);
318 prec = new OperatorJacobiSmoother(M, ess_tdof_list);
319 }
320 else
321 {
322 Mop.Reset(M.ParallelAssemble());
323 prec = new HypreSmoother;
324 static_cast<HypreSmoother*>(prec)->SetType(HypreSmoother::Jacobi, 1);
325 }
326
327 CGSolver lin_solver(M.ParFESpace()->GetParMesh()->GetComm());
328 lin_solver.SetPreconditioner(*prec);
329 lin_solver.SetOperator(*Mop);
330#ifdef MFEM_USE_SINGLE
331 const real_t rtol = 1e-4;
332#else
333 const real_t rtol = 1e-8;
334#endif
335 lin_solver.SetRelTol(rtol); lin_solver.SetAbsTol(0.0);
336 lin_solver.SetMaxIter(100);
337 lin_solver.SetPrintLevel(0);
338 lin_solver.Mult(*RHS, X);
339 K.ParFESpace()->GetProlongationMatrix()->Mult(X, di_dt);
340
341 delete RHS;
342 delete prec;
343}
344#endif
345
346#ifdef MFEM_USE_GSLIB
348 const Vector &init_field)
349{
350 nodes0 = init_nodes;
351 Mesh *m = mesh;
353#ifdef MFEM_USE_MPI
354 if (pmesh) { m = pmesh; }
355 if (pfes) { f = pfes; }
356#endif
357 m->SetNodes(nodes0);
358
359 if (m->GetNodes()->FESpace()->IsDGSpace())
360 {
361 MFEM_ABORT("InterpolatorFP is not supported for periodic meshes yet.");
362 }
363
364 const real_t rel_bbox_el = 0.1;
365 const real_t newton_tol = 1.0e-12;
366 const int npts_at_once = 256;
367
368 if (finder)
369 {
370 finder->FreeData();
371 delete finder;
372 }
373
374#ifdef MFEM_USE_MPI
375 if (pfes) { finder = new FindPointsGSLIB(pfes->GetComm()); }
376 else { finder = new FindPointsGSLIB(); }
377#else
378 finder = new FindPointsGSLIB();
379#endif
380 finder->Setup(*m, rel_bbox_el, newton_tol, npts_at_once);
381
382 field0_gf.SetSpace(f);
383 field0_gf = init_field;
384}
385
387 Vector &new_field, int nodes_ordering)
388{
389 // TODO - this is here only to prevent breaking user codes. To be removed.
390 // If the meshes are different, one has to call SetNewFieldFESpace().
391 // If only some positions are interpolated, use ComputeAtGivenPositions().
392 if (fes_new_field == nullptr && new_mesh_nodes.Size() != nodes0.Size())
393 {
394 MFEM_WARNING("Deprecated -- use ComputeAtGivenPositions() instead!");
395 ComputeAtGivenPositions(new_mesh_nodes, new_field, nodes_ordering);
396 return;
397 }
398
399 const FiniteElementSpace *fes_field =
400 (fes_new_field) ? fes_new_field : field0_gf.FESpace();
401 const int dim = fes_field->GetMesh()->Dimension();
402
403 if (new_mesh_nodes.Size() / dim != fes_field->GetNDofs())
404 {
405 // The nodes of the FE space don't coincide with the mesh nodes.
406 Vector mapped_nodes;
407 fes_field->GetNodePositions(new_mesh_nodes, mapped_nodes);
408 finder->Interpolate(mapped_nodes, field0_gf, new_field);
409 }
410 else
411 {
412 finder->Interpolate(new_mesh_nodes, field0_gf, new_field, nodes_ordering);
413 }
414}
415
417 Vector &values, int p_ordering)
418{
419 finder->Interpolate(positions, field0_gf, values, p_ordering);
420}
421
422#endif
423
425 const Vector &b) const
426{
427 const FiniteElementSpace *fes = NULL;
428 real_t energy_in = 0.0;
429#ifdef MFEM_USE_MPI
430 const ParNonlinearForm *p_nlf = dynamic_cast<const ParNonlinearForm *>(oper);
431 MFEM_VERIFY(!(parallel && p_nlf == NULL), "Invalid Operator subclass.");
432 if (parallel)
433 {
434 fes = p_nlf->FESpace();
435 energy_in = p_nlf->GetEnergy(d_in);
436 }
437#endif
438 const bool serial = !parallel;
439 const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
440 MFEM_VERIFY(!(serial && nlf == NULL), "Invalid Operator subclass.");
441 if (serial)
442 {
443 fes = nlf->FESpace();
444 energy_in = nlf->GetEnergy(d_in);
445 }
446
447 // Get the local prolongation of the solution vector.
448 Vector d_loc(fes->GetVSize(), (temp_mt == MemoryType::DEFAULT) ?
450 if (serial)
451 {
452 const SparseMatrix *cP = fes->GetConformingProlongation();
453 if (!cP) { d_loc = d_in; }
454 else { cP->Mult(d_in, d_loc); }
455 }
456#ifdef MFEM_USE_MPI
457 else
458 {
459 fes->GetProlongationMatrix()->Mult(d_in, d_loc);
460 }
461#endif
462
463 real_t scale = 1.0;
464 bool fitting = IsSurfaceFittingEnabled();
465 real_t init_fit_avg_err, init_fit_max_err = 0.0;
466 if (fitting && surf_fit_converge_error)
467 {
468 GetSurfaceFittingError(d_loc, init_fit_avg_err, init_fit_max_err);
469 // Check for convergence
470 if (init_fit_max_err < surf_fit_max_err_limit)
471 {
473 {
474 mfem::out << "TMOPNewtonSolver converged "
475 "based on the surface fitting error.\n";
476 }
477 scale = 0.0;
478 return scale;
479 }
480 }
481
483 {
485 {
486 mfem::out << "TMOPNewtonSolver terminated "
487 "based on max number of times surface fitting weight can"
488 "be increased. \n";
489 }
490 scale = 0.0;
491 return scale;
492 }
493
494 // Check if the starting mesh (given by x) is inverted. Note that x hasn't
495 // been modified by the Newton update yet.
496 const real_t min_detT_in = ComputeMinDet(d_loc, *fes);
497 const bool untangling = (min_detT_in <= 0.0) ? true : false;
498 const real_t untangle_factor = 1.5;
499 if (untangling)
500 {
501 // Needed for the line search below. The untangling metrics see this
502 // reference to detect deteriorations.
503 MFEM_VERIFY(min_det_ptr != NULL, " Initial mesh was valid, but"
504 " intermediate mesh is invalid. Contact TMOP Developers.");
505 MFEM_VERIFY(min_detJ_limit == 0.0,
506 "This setup is not supported. Contact TMOP Developers.");
507 *min_det_ptr = untangle_factor * min_detT_in;
508 }
509
510 const bool have_b = (b.Size() == Height());
511
512 Vector d_out(d_in.Size());
513 bool x_out_ok = false;
514 real_t energy_out = 0.0, min_detT_out;
515 const real_t norm_in = Norm(r);
516 real_t avg_fit_err, max_fit_err = 0.0;
517
518 const real_t detJ_factor = (solver_type == 1) ? 0.25 : 0.5;
520 // TODO:
521 // - Customized line search for worst-quality optimization.
522 // - What is the Newton exit criterion for worst-quality optimization?
523
524 // Perform the line search.
525 for (int i = 0; i < 12; i++)
526 {
527 avg_fit_err = 0.0;
528 max_fit_err = 0.0;
529
530 //
531 // Update the mesh and get the L-vector in x_out_loc.
532 //
533 // Form limited (line-search) displacement d_out = d_in - scale * c,
534 // and the corresponding mesh positions x_out = x_0 + d_out.
535 add(d_in, -scale, c, d_out);
536 if (serial)
537 {
538 const SparseMatrix *cP = fes->GetConformingProlongation();
539 if (!cP) { d_loc = d_out; }
540 else { cP->Mult(d_out, d_loc); }
541 }
542#ifdef MFEM_USE_MPI
543 else { fes->GetProlongationMatrix()->Mult(d_out, d_loc); }
544#endif
545
546 // Check the changes in detJ.
547 min_detT_out = ComputeMinDet(d_loc, *fes);
548 if (untangling == false && min_detT_out <= min_detJ_limit)
549 {
550 // No untangling, and detJ got negative (or small) -- no good.
552 {
553 mfem::out << "Scale = " << scale << " Neg det(J) found.\n";
554 }
555 scale *= detJ_factor; continue;
556 }
557 if (untangling == true && min_detT_out < *min_det_ptr)
558 {
559 // Untangling, and detJ got even more negative -- no good.
561 {
562 mfem::out << "Scale = " << scale << " Neg det(J) decreased.\n";
563 }
564 scale *= detJ_factor; continue;
565 }
566
567 // Skip the energy and residual checks when we're untangling. The
568 // untangling metrics change their denominators, which can affect the
569 // energy and residual, so their increase/decrease is not relevant.
570 if (untangling) { x_out_ok = true; break; }
571
572 // Update mesh-dependent quantities.
573 ProcessNewState(d_out);
574
575 // Ensure sufficient decrease in fitting error if we are trying to
576 // converge based on error.
577 if (fitting && surf_fit_converge_error)
578 {
579 GetSurfaceFittingError(d_loc, avg_fit_err, max_fit_err);
580 if (max_fit_err >= 1.2*init_fit_max_err)
581 {
583 {
584 mfem::out << "Scale = " << scale << " Surf fit err increased.\n";
585 }
586 scale *= 0.5; continue;
587 }
588 }
589
590 // Check the changes in total energy.
591 if (serial)
592 {
593 energy_out = nlf->GetEnergy(d_out);
594 }
595#ifdef MFEM_USE_MPI
596 else
597 {
598 energy_out = p_nlf->GetEnergy(d_out);
599 }
600#endif
601 if (energy_out > energy_in + 0.2*fabs(energy_in) ||
602 std::isnan(energy_out) != 0)
603 {
605 {
606 mfem::out << "Scale = " << scale << " Increasing energy: "
607 << energy_in << " --> " << energy_out << '\n';
608 }
609 scale *= 0.5; continue;
610 }
611
612 // Check the changes in the Newton residual.
613 oper->Mult(d_out, r);
614 if (have_b) { r -= b; }
615 real_t norm_out = Norm(r);
616
617 if (norm_out > 1.2*norm_in)
618 {
620 {
621 mfem::out << "Scale = " << scale << " Norm increased: "
622 << norm_in << " --> " << norm_out << '\n';
623 }
624 scale *= 0.5; continue;
625 }
626 else { x_out_ok = true; break; }
627 } // end line search
628
629 if (untangling)
630 {
631 // Update the global min detJ. Untangling metrics see this min_det_ptr.
632 if (min_detT_out > 0.0)
633 {
634 *min_det_ptr = 0.0;
637 { mfem::out << "The mesh has been untangled at the used points!\n"; }
638 }
639 else { *min_det_ptr = untangle_factor * min_detT_out; }
640 }
641
644 {
645 if (untangling)
646 {
647 mfem::out << "Min det(T) change: "
648 << min_detT_in << " -> " << min_detT_out
649 << " with " << scale << " scaling.\n";
650 }
651 else
652 {
653 mfem::out << "Energy decrease: "
654 << energy_in << " --> " << energy_out << " or "
655 << (energy_in - energy_out) / energy_in * 100.0
656 << "% with " << scale << " scaling.\n";
657 }
658 }
659
660 if (x_out_ok == false) { scale = 0.0; }
661
662 if (surf_fit_scale_factor > 0.0) { surf_fit_coeff_update = true; }
664
665 return scale;
666}
667
668void TMOPNewtonSolver::Mult(const Vector &b, Vector &x) const
669{
670 // Prolongate x to ldofs.
671 const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
672 auto fes_mesh_nodes = nlf->FESpace()->GetMesh()->GetNodes()->FESpace();
673 const Operator *P = fes_mesh_nodes->GetProlongationMatrix();
674 x_0.SetSpace(fes_mesh_nodes);
675 periodic = fes_mesh_nodes->IsDGSpace();
676 if (P)
677 {
678 MFEM_VERIFY(x.Size() == P->Width(),
679 "The input's size must be the tdof size of the mesh nodes.");
680 P->Mult(x, x_0);
681 }
682 else
683 {
684 MFEM_VERIFY(x.Size() == x_0.Size(),
685 "The input's size must match the size of the mesh nodes.");
686 x_0 = x;
687 }
688
689 // Pass down the initial position to the integrators.
690 const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
691 for (int i = 0; i < integs.Size(); i++)
692 {
693 auto ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
694 if (ti) { ti->SetInitialMeshPos(&x_0); }
695 auto co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
696 if (co) { co->SetInitialMeshPos(&x_0); }
697 }
698
699 // Solve for the displacement, which always starts from zero.
700 Vector dx(height); dx = 0.0;
701 if (solver_type == 0) { NewtonSolver::Mult(b, dx); }
702 else if (solver_type == 1) { LBFGSSolver::Mult(b, dx); }
703 else { MFEM_ABORT("Invalid solver_type"); }
704
705 // Form the final mesh using the computed displacement.
706 if (periodic)
707 {
708 Vector dx_loc(nlf->FESpace()->GetVSize());
709 const Operator *Pd = nlf->FESpace()->GetProlongationMatrix();
710 if (Pd) { Pd->Mult(dx, dx_loc); }
711 else { dx_loc = dx; }
712
713 GetPeriodicPositions(x_0, dx_loc, *fes_mesh_nodes, *nlf->FESpace(), x);
714 }
715 else { x += dx; }
716
717 // Make sure the pointers don't use invalid memory (x_0_loc is gone).
718 for (int i = 0; i < integs.Size(); i++)
719 {
720 auto ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
721 if (ti) { ti->SetInitialMeshPos(nullptr); }
722 auto co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
723 if (co) { co->SetInitialMeshPos(nullptr); }
724 }
725}
726
728{
729 const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
730 const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
731 TMOP_Integrator *ti = NULL;
732 TMOPComboIntegrator *co = NULL;
733 for (int i = 0; i < integs.Size(); i++)
734 {
735 ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
736 if (ti)
737 {
738 ti->UpdateSurfaceFittingWeight(factor);
739 }
740 co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
741 if (co)
742 {
744 for (int j = 0; j < ati.Size(); j++)
745 {
746 ati[j]->UpdateSurfaceFittingWeight(factor);
747 }
748 }
749 }
750}
751
753{
754 const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
755 const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
756 TMOP_Integrator *ti = NULL;
757 TMOPComboIntegrator *co = NULL;
758 weights.SetSize(0);
759 real_t weight;
760
761 for (int i = 0; i < integs.Size(); i++)
762 {
763 ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
764 if (ti && ti->IsSurfaceFittingEnabled())
765 {
766 weight = ti->GetSurfaceFittingWeight();
767 weights.Append(weight);
768 }
769 co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
770 if (co)
771 {
773 for (int j = 0; j < ati.Size(); j++)
774 {
775 if (ati[j]->IsSurfaceFittingEnabled())
776 {
777 weight = ati[j]->GetSurfaceFittingWeight();
778 weights.Append(weight);
779 }
780 }
781 }
782 }
783}
784
786 real_t &err_avg,
787 real_t &err_max) const
788{
789 const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
790 const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
791 TMOP_Integrator *ti = NULL;
792 TMOPComboIntegrator *co = NULL;
793
794 err_avg = 0.0;
795 err_max = 0.0;
796 real_t err_avg_loc, err_max_loc;
797 for (int i = 0; i < integs.Size(); i++)
798 {
799 ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
800 if (ti)
801 {
802 if (ti->IsSurfaceFittingEnabled())
803 {
804 ti->GetSurfaceFittingErrors(d_loc, err_avg_loc, err_max_loc);
805 err_avg = std::max(err_avg_loc, err_avg);
806 err_max = std::max(err_max_loc, err_max);
807 }
808 }
809 co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
810 if (co)
811 {
813 for (int j = 0; j < ati.Size(); j++)
814 {
815 if (ati[j]->IsSurfaceFittingEnabled())
816 {
817 ati[j]->GetSurfaceFittingErrors(d_loc, err_avg_loc, err_max_loc);
818 err_avg = std::max(err_avg_loc, err_avg);
819 err_max = std::max(err_max_loc, err_max);
820 }
821 }
822 }
823 }
824}
825
827{
828 const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
829 const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
830 TMOP_Integrator *ti = NULL;
831 TMOPComboIntegrator *co = NULL;
832
833 for (int i = 0; i < integs.Size(); i++)
834 {
835 ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
836 if (ti)
837 {
838 if (ti->IsSurfaceFittingEnabled())
839 {
840 return true;
841 }
842 }
843 co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
844 if (co)
845 {
847 for (int j = 0; j < ati.Size(); j++)
848 {
849 if (ati[j]->IsSurfaceFittingEnabled())
850 {
851 return true;
852 }
853 }
854 }
855 }
856 return false;
857}
858
860{
861 const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
862 const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
863
864 // Reset the update flags of all TargetConstructors. This is done to avoid
865 // repeated updates of shared TargetConstructors.
866 TMOP_Integrator *ti = NULL;
867 TMOPComboIntegrator *co = NULL;
868 DiscreteAdaptTC *dtc = NULL;
869 for (int i = 0; i < integs.Size(); i++)
870 {
871 ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
872 if (ti)
873 {
874 dtc = ti->GetDiscreteAdaptTC();
875 if (dtc) { dtc->ResetUpdateFlags(); }
876 }
877 co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
878 if (co)
879 {
881 for (int j = 0; j < ati.Size(); j++)
882 {
883 dtc = ati[j]->GetDiscreteAdaptTC();
884 if (dtc) { dtc->ResetUpdateFlags(); }
885 }
886 }
887 }
888
889 Vector dx_loc;
890 const Operator *P = nlf->GetProlongation();
891 if (P)
892 {
893 dx_loc.SetSize(P->Height());
894 P->Mult(dx, dx_loc);
895 }
896 else { dx_loc = dx; }
897
898 const FiniteElementSpace *dx_fes = nlf->FESpace();
899 for (int i = 0; i < integs.Size(); i++)
900 {
901 ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
902 if (ti)
903 {
904 ti->UpdateAfterMeshPositionChange(dx_loc, *dx_fes);
906 {
907 ti->ComputeUntangleMetricQuantiles(dx_loc, *dx_fes);
908 }
909 }
910 co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
911 if (co)
912 {
914 for (int j = 0; j < ati.Size(); j++)
915 {
916 ati[j]->UpdateAfterMeshPositionChange(dx_loc, *dx_fes);
918 {
919 ati[j]->ComputeUntangleMetricQuantiles(dx_loc, *dx_fes);
920 }
921 }
922 }
923 }
924
925 // Constant coefficient associated with the surface fitting terms if
926 // adaptive surface fitting is enabled. The idea is to increase the
927 // coefficient if the surface fitting error does not sufficiently
928 // decrease between subsequent TMOPNewtonSolver iterations.
930 {
931 // Get surface fitting errors.
933 // Get array with surface fitting weights.
934 Array<real_t> fitweights;
935 GetSurfaceFittingWeight(fitweights);
936
938 {
939 mfem::out << "Avg/Max surface fitting error: " <<
940 surf_fit_avg_err << " " <<
941 surf_fit_max_err << "\n";
942 mfem::out << "Min/Max surface fitting weight: " <<
943 fitweights.Min() << " " << fitweights.Max() << "\n";
944 }
945
946 real_t change_surf_fit_err = surf_fit_avg_err_prvs-surf_fit_avg_err;
947 real_t rel_change_surf_fit_err = change_surf_fit_err/surf_fit_avg_err_prvs;
948
949 // Increase the surface fitting coefficient if the surface fitting error
950 // does not decrease sufficiently. If we are converging based on residual,
951 // also make sure we have not reached the maximum fitting weight and
952 // error threshold.
953 if (rel_change_surf_fit_err < surf_fit_err_rel_change_limit &&
955 (fitweights.Max() < surf_fit_weight_limit &&
957 {
958 real_t scale_factor = std::min(surf_fit_scale_factor,
959 surf_fit_weight_limit/fitweights.Max());
960 UpdateSurfaceFittingWeight(scale_factor);
962 }
963 else
964 {
966 }
968 surf_fit_coeff_update = false;
969 }
970}
971
973 const FiniteElementSpace &fes) const
974{
975 real_t min_detJ = infinity();
976 const int NE = fes.GetNE(), dim = fes.GetMesh()->Dimension();
977 Array<int> xdofs;
978 DenseMatrix Jpr(dim);
979 const bool mixed_mesh = fes.GetMesh()->GetNumGeometries(dim) > 1;
980 if (dim == 1 || mixed_mesh ||
981 UsesTensorBasis(fes) == false || fes.IsVariableOrder())
982 {
983 for (int i = 0; i < NE; i++)
984 {
985 const int dof = fes.GetFE(i)->GetDof();
986 DenseMatrix dshape(dof, dim), pos(dof, dim);
987 Vector posV(pos.Data(), dof * dim);
988
989 x_0.GetElementDofValues(i, posV);
990 if (periodic)
991 {
992 auto n_el = dynamic_cast<const NodalFiniteElement *>(fes.GetFE(i));
993 n_el->ReorderLexToNative(dim, posV);
994 }
995
996 Vector d_loc_el;
997 fes.GetElementVDofs(i, xdofs);
998 d_loc.GetSubVector(xdofs, d_loc_el);
999 posV += d_loc_el;
1000
1001 const IntegrationRule &irule = GetIntegrationRule(*fes.GetFE(i));
1002 const int nsp = irule.GetNPoints();
1003 for (int j = 0; j < nsp; j++)
1004 {
1005 fes.GetFE(i)->CalcDShape(irule.IntPoint(j), dshape);
1006 MultAtB(pos, dshape, Jpr);
1007 min_detJ = std::min(min_detJ, Jpr.Det());
1008 }
1009 }
1010 }
1011 else
1012 {
1013 min_detJ = dim == 2 ? MinDetJpr_2D(&fes, d_loc) :
1014 dim == 3 ? MinDetJpr_3D(&fes, d_loc) : 0.0;
1015 }
1016#ifdef MFEM_USE_MPI
1017 if (parallel)
1018 {
1019 auto p_nlf = dynamic_cast<const ParNonlinearForm *>(oper);
1020 MPI_Allreduce(MPI_IN_PLACE, &min_detJ, 1, MPITypeMap<real_t>::mpi_type,
1021 MPI_MIN, p_nlf->ParFESpace()->GetComm());
1022 }
1023#endif
1024 const DenseMatrix &Wideal =
1026 min_detJ /= Wideal.Det();
1027
1028 return min_detJ;
1029}
1030
1031#ifdef MFEM_USE_MPI
1032// Metric values are visualized by creating an L2 finite element functions and
1033// computing the metric values at the nodes.
1035 const TargetConstructor &tc, ParMesh &pmesh,
1036 char *title, int position)
1037{
1039 ParFiniteElementSpace fes(&pmesh, &fec, 1);
1040 ParGridFunction metric(&fes);
1041 InterpolateTMOP_QualityMetric(qm, tc, pmesh, metric);
1042 socketstream sock;
1043 if (pmesh.GetMyRank() == 0)
1044 {
1045 sock.open("localhost", 19916);
1046 sock << "solution\n";
1047 }
1048 pmesh.PrintAsOne(sock);
1049 metric.SaveAsOne(sock);
1050 if (pmesh.GetMyRank() == 0)
1051 {
1052 sock << "window_title '"<< title << "'\n"
1053 << "window_geometry "
1054 << position << " " << 0 << " " << 600 << " " << 600 << "\n"
1055 << "keys jRmclA\n";
1056 }
1057}
1058#endif
1059
1060// Metric values are visualized by creating an L2 finite element functions and
1061// computing the metric values at the nodes.
1063 const TargetConstructor &tc, Mesh &mesh,
1064 char *title, int position)
1065{
1067 FiniteElementSpace fes(&mesh, &fec, 1);
1068 GridFunction metric(&fes);
1069 InterpolateTMOP_QualityMetric(qm, tc, mesh, metric);
1070 osockstream sock(19916, "localhost");
1071 sock << "solution\n";
1072 mesh.Print(sock);
1073 metric.Save(sock);
1074 sock.send();
1075 sock << "window_title '"<< title << "'\n"
1076 << "window_geometry "
1077 << position << " " << 0 << " " << 600 << " " << 600 << "\n"
1078 << "keys jRmclA\n";
1079}
1080
1081void GetPeriodicPositions(const Vector &x_0, const Vector &dx,
1082 const FiniteElementSpace &fesL2,
1083 const FiniteElementSpace &fesH1, Vector &x)
1084{
1085 x = x_0;
1086 Vector dx_r(x.Size());
1088 auto R_H1 = fesH1.GetElementRestriction(ord);
1089 auto R_L2 = fesL2.GetElementRestriction(ord);
1090 R_H1->Mult(dx, dx_r);
1091 R_L2->AddMultTranspose(dx_r, x);
1092}
1093
1094}
FiniteElementSpace * fes
Definition tmop.hpp:1409
ParFiniteElementSpace * pfes
Definition tmop.hpp:1414
void SetInitialField(const Vector &init_nodes, const Vector &init_field) override
void ComputeAtNewPosition(const Vector &new_mesh_nodes, Vector &new_field, int nodes_ordering=Ordering::byNODES) override
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:758
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
Definition array.cpp:85
int Size() const
Return the logical size of the array.
Definition array.hpp:147
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
@ GaussLobatto
Closed type.
Definition fe_base.hpp:36
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void Finalize(int skip_zeros=1) override
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x) override
Recover the solution of a linear system formed with FormLinearSystem().
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(....
int Size() const
Get the size of the BilinearForm as a square matrix.
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication: .
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
Conjugate gradient method.
Definition solvers.hpp:538
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:751
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Definition solvers.hpp:551
Data type for scaled Jacobi-type smoother of sparse matrix.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:114
real_t Det() const
Definition densemat.cpp:516
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition device.hpp:274
void ResetUpdateFlags()
Used in combination with the Update methods to avoid extra computations.
Definition tmop.hpp:1761
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points.
Definition gslib.hpp:67
void Setup(Mesh &m, const double bb_t=0.1, const double newt_tol=1.0e-12, const int npt_max=256)
Definition gslib.cpp:163
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition gslib.cpp:1721
virtual void FreeData()
Definition gslib.cpp:1128
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition fespace.hpp:709
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition fespace.hpp:845
virtual const Operator * GetProlongationMatrix() const
Definition fespace.hpp:727
void GetNodePositions(const Vector &mesh_nodes, Vector &fes_node_pos, int fes_nodes_ordering=Ordering::byNODES) const
Compute the space's node positions w.r.t. given mesh positions. The function uses FiniteElement::GetN...
Definition fespace.cpp:4332
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:332
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3835
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:876
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:891
const ElementRestrictionOperator * GetElementRestriction(ElementDofOrdering e_ordering) const
Return an Operator that converts L-vectors to E-vectors.
Definition fespace.cpp:1510
const SparseMatrix * GetConformingProlongation() const
Definition fespace.cpp:1456
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:878
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:848
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition geom.hpp:102
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual void GetElementDofValues(int el, Vector &dof_vals) const
FiniteElementSpace * FESpace()
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition gridfunc.cpp:225
Wrapper for hypre's parallel vector class.
Definition hypre.hpp:219
Parallel smoothers in hypre.
Definition hypre.hpp:1046
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
void ComputeAtGivenPositions(const Vector &positions, Vector &values, int p_ordering=Ordering::byNODES) override
Direct interpolation of field0_gf at the given positions.
void SetInitialField(const Vector &init_nodes, const Vector &init_field) override
void ComputeAtNewPosition(const Vector &new_mesh_nodes, Vector &new_field, int nodes_ordering=Ordering::byNODES) override
PrintLevel print_options
Output behavior for the iterative solver.
Definition solvers.hpp:157
const Operator * oper
Definition solvers.hpp:140
void SetRelTol(real_t rtol)
Definition solvers.hpp:229
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:174
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:72
void SetMaxIter(int max_it)
Definition solvers.hpp:231
void SetAbsTol(real_t atol)
Definition solvers.hpp:230
real_t Norm(const Vector &x) const
Return the inner product norm of x, using the inner product defined by Dot()
Definition solvers.hpp:200
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:346
void Mult(const Vector &b, Vector &x) const override
Solve the nonlinear system with right-hand side b.
Definition solvers.cpp:2103
Mesh data type.
Definition mesh.hpp:64
void NodesUpdated()
This function should be called after the mesh node coordinates have been updated externally,...
Definition mesh.hpp:2211
Geometry::Type GetTypicalElementGeometry() const
If the local mesh is not empty, return GetElementGeometry(0); otherwise, return a typical Geometry pr...
Definition mesh.cpp:1528
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2433
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
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 GetNodes(Vector &node_coord) const
Definition mesh.cpp:9294
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
Definition mesh.cpp:9306
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:7243
void Mult(const Vector &b, Vector &x) const override
Solve the nonlinear system with right-hand side b.
Definition solvers.cpp:1912
Class for standard nodal finite elements.
Definition fe_base.hpp:721
void ReorderLexToNative(int ncomp, Vector &dofs) const
Definition fe_base.cpp:976
void SetPAMemoryType(MemoryType mt)
FiniteElementSpace * FESpace()
virtual real_t GetEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
const Operator * GetProlongation() const override
Get the finite element space prolongation matrix.
Array< NonlinearFormIntegrator * > * GetDNFI()
Access all integrators added with AddDomainIntegrator().
Pointer to an Operator of a specified type.
Definition handle.hpp:34
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition handle.hpp:145
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:333
Abstract operator.
Definition operator.hpp:25
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
int height
Dimension of the output / number of rows in the matrix.
Definition operator.hpp:27
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Performs a single remap advection step in parallel.
ParAdvectorCGOper(const Vector &x_start, GridFunction &vel, ParFiniteElementSpace &pfes, AssemblyLevel al=AssemblyLevel::LEGACY, MemoryType mt=MemoryType::DEFAULT)
VectorGridFunctionCoefficient u_coeff
void Mult(const Vector &ind, Vector &di_dt) const override
Operator application: y=A(x).
const AssemblyLevel al
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A) override
Form the linear system matrix A, see FormLinearSystem() for details.
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:334
const Operator * GetProlongationMatrix() const override
ParMesh * GetParMesh() const
Definition pfespace.hpp:338
Class for parallel grid function.
Definition pgridfunc.hpp:50
void ParallelAssemble(Vector &tv) const
Returns the vector assembled on the true dofs.
void SaveAsOne(const char *fname, int precision=16) const
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
int GetMyRank() const
Definition pmesh.hpp:404
void PrintAsOne(std::ostream &out=mfem::out, const std::string &comments="") const
Definition pmesh.cpp:4967
Parallel non-linear operator on the true dofs.
real_t GetEnergy(const ParGridFunction &x) const
Compute the energy of a ParGridFunction.
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition ode.cpp:243
void Step(Vector &x, real_t &t, real_t &dt) override
Perform a time step from time t [in] to time t [out] based on the requested step size dt [in].
Definition ode.cpp:252
Performs a single remap advection step in serial.
VectorGridFunctionCoefficient u_coeff
SerialAdvectorCGOper(const Vector &x_start, GridFunction &vel, FiniteElementSpace &fes, AssemblyLevel al=AssemblyLevel::LEGACY)
const AssemblyLevel al
void Mult(const Vector &ind, Vector &di_dt) const override
Operator application: y=A(x).
Base class for solvers.
Definition operator.hpp:780
Data type sparse matrix.
Definition sparsemat.hpp:51
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
void SetInitialMeshPos(const GridFunction *x0)
Definition tmop.hpp:2438
const Array< TMOP_Integrator * > & GetTMOPIntegrators() const
Definition tmop.hpp:2455
real_t ComputeScalingFactor(const Vector &d, const Vector &b) const override
virtual void GetSurfaceFittingError(const Vector &d_loc, real_t &err_avg, real_t &err_max) const
void Mult(const Vector &b, Vector &x) const override
Optimizes the mesh positions given by x.
bool IsSurfaceFittingEnabled() const
Check if surface fitting is enabled.
real_t MinDetJpr_3D(const FiniteElementSpace *, const Vector &) const
const IntegrationRule & GetIntegrationRule(const FiniteElement &el) const
void UpdateSurfaceFittingWeight(real_t factor) const
Update surface fitting weight as surf_fit_weight *= factor.
real_t MinDetJpr_2D(const FiniteElementSpace *, const Vector &) const
void ProcessNewState(const Vector &dx) const override
real_t ComputeMinDet(const Vector &d_loc, const FiniteElementSpace &fes) const
void GetSurfaceFittingWeight(Array< real_t > &weights) const
Get the surface fitting weight for all the TMOP integrators.
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition tmop.hpp:1885
void UpdateAfterMeshPositionChange(const Vector &d, const FiniteElementSpace &d_fes)
Definition tmop.cpp:5364
void GetSurfaceFittingErrors(const Vector &d_loc, real_t &err_avg, real_t &err_max)
Definition tmop.cpp:3829
void ComputeUntangleMetricQuantiles(const Vector &d, const FiniteElementSpace &fes)
Definition tmop.cpp:5593
real_t GetSurfaceFittingWeight()
Get the surface fitting weight.
Definition tmop.cpp:5093
void SetInitialMeshPos(const GridFunction *x0)
Definition tmop.cpp:3496
bool IsSurfaceFittingEnabled()
Definition tmop.hpp:2332
DiscreteAdaptTC * GetDiscreteAdaptTC() const
Definition tmop.hpp:2395
void UpdateSurfaceFittingWeight(real_t factor)
Update the surface fitting weight as surf_fit_coeff *= factor;.
Definition tmop.cpp:5081
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:1481
Base abstract class for first order time dependent operators.
Definition operator.hpp:332
real_t t
Current time.
Definition operator.hpp:373
virtual real_t GetTime() const
Read the currently set time.
Definition operator.hpp:391
Vector data type.
Definition vector.hpp:82
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:1164
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:514
real_t Min() const
Returns the minimal element of the vector.
Definition vector.cpp:1123
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:653
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:622
int open(const char hostname[], int port)
Open the socket stream on 'port' at 'hostname'.
int dim
Definition ex24.cpp:53
Mesh * GetMesh(int type)
Definition ex29.cpp:218
real_t b
Definition lissajous.cpp:42
string space
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:47
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
Geometry Geometries
Definition fe.cpp:49
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:391
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric, const TargetConstructor &tc, const Mesh &mesh, GridFunction &metric_gf)
Interpolates the metric's values at the nodes of metric_gf.
Definition tmop.cpp:5823
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:949
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
void GetPeriodicPositions(const Vector &x_0, const Vector &dx, const FiniteElementSpace &fesL2, const FiniteElementSpace &fesH1, Vector &x)
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition fespace.hpp:1555
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition vector.cpp:547
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
float real_t
Definition config.hpp:43
MemoryType
Memory types supported by MFEM.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition fespace.hpp:83
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void vel(const Vector &x, real_t t, Vector &u)
bool iterations
Detailed information about each iteration will be reported to mfem::out.
Definition solvers.hpp:107
bool warnings
If a non-fatal problem has been detected some context-specific information will be reported to mfem::...
Definition solvers.hpp:104
bool first_and_last
Information about the first and last iteration will be printed to mfem::out.
Definition solvers.hpp:113
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Definition solvers.hpp:110
Helper struct to convert a C++ type to an MPI type.