MFEM  v4.5.2
Finite element discretization library
tmop_tools.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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"
15 #include "../general/osockstream.hpp"
16 
17 namespace mfem
18 {
19 
20 using namespace mfem;
21 
22 void AdvectorCG::SetInitialField(const Vector &init_nodes,
23  const Vector &init_field)
24 {
25  nodes0 = init_nodes;
26  field0 = init_field;
27 }
28 
30  Vector &new_field,
31  int new_nodes_ordering)
32 {
34 #ifdef MFEM_USE_MPI
35  if (pfes) { space = pfes; }
36 #endif
37  int fes_ordering = space->GetOrdering(),
38  ncomp = space->GetVDim();
39 
40  // TODO: Implement for AMR meshes.
41  const int pnt_cnt = field0.Size() / ncomp;
42 
43  new_field = field0;
44  Vector new_field_temp;
45  for (int i = 0; i < ncomp; i++)
46  {
47  if (fes_ordering == Ordering::byNODES)
48  {
49  new_field_temp.MakeRef(new_field, i*pnt_cnt, pnt_cnt);
50  }
51  else
52  {
53  new_field_temp.SetSize(pnt_cnt);
54  for (int j = 0; j < pnt_cnt; j++)
55  {
56  new_field_temp(j) = new_field(i + j*ncomp);
57  }
58  }
59  ComputeAtNewPositionScalar(new_nodes, new_field_temp);
60  if (fes_ordering == Ordering::byVDIM)
61  {
62  for (int j = 0; j < pnt_cnt; j++)
63  {
64  new_field(i + j*ncomp) = new_field_temp(j);
65  }
66  }
67  }
68 
69  field0 = new_field;
70  nodes0 = new_nodes;
71 }
72 
73 void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
74  Vector &new_field)
75 {
76  Mesh *m = mesh;
77 #ifdef MFEM_USE_MPI
78  if (pmesh) { m = pmesh; }
79 #endif
80 
81  MFEM_VERIFY(m != NULL, "No mesh has been given to the AdaptivityEvaluator.");
82 
83  // This will be used to move the positions.
84  GridFunction *mesh_nodes = m->GetNodes();
85  *mesh_nodes = nodes0;
86  double minv = new_field.Min(), maxv = new_field.Max();
87 
88  // Velocity of the positions.
89  GridFunction u(mesh_nodes->FESpace());
90  subtract(new_nodes, nodes0, u);
91 
92  // Define a scalar FE space for the solution, and the advection operator.
93  TimeDependentOperator *oper = NULL;
94  FiniteElementSpace *fess = NULL;
95 #ifdef MFEM_USE_MPI
96  ParFiniteElementSpace *pfess = NULL;
97 #endif
98  if (fes)
99  {
100  fess = new FiniteElementSpace(fes->GetMesh(), fes->FEColl(), 1);
101  oper = new SerialAdvectorCGOper(nodes0, u, *fess, al);
102  }
103 #ifdef MFEM_USE_MPI
104  else if (pfes)
105  {
106  pfess = new ParFiniteElementSpace(pfes->GetParMesh(), pfes->FEColl(), 1);
107  oper = new ParAdvectorCGOper(nodes0, u, *pfess, al, opt_mt);
108  }
109 #endif
110  MFEM_VERIFY(oper != NULL,
111  "No FE space has been given to the AdaptivityEvaluator.");
112  ode_solver.Init(*oper);
113 
114  // Compute some time step [mesh_size / speed].
115  double h_min = std::numeric_limits<double>::infinity();
116  for (int i = 0; i < m->GetNE(); i++)
117  {
118  h_min = std::min(h_min, m->GetElementSize(i));
119  }
120  double v_max = 0.0;
121  const int s = new_field.Size();
122 
123  u.HostReadWrite();
124  for (int i = 0; i < s; i++)
125  {
126  double vel = 0.;
127  for (int j = 0; j < m->Dimension(); j++)
128  {
129  vel += u(i+j*s)*u(i+j*s);
130  }
131  v_max = std::max(v_max, vel);
132  }
133 
134 #ifdef MFEM_USE_MPI
135  if (pfes)
136  {
137  double v_loc = v_max, h_loc = h_min;
138  MPI_Allreduce(&v_loc, &v_max, 1, MPI_DOUBLE, MPI_MAX, pfes->GetComm());
139  MPI_Allreduce(&h_loc, &h_min, 1, MPI_DOUBLE, MPI_MIN, pfes->GetComm());
140  }
141 #endif
142 
143  if (v_max == 0.0) // No need to change the field.
144  {
145  delete oper;
146  delete fess;
147 #ifdef MFEM_USE_MPI
148  delete pfess;
149 #endif
150  return;
151  }
152 
153  v_max = std::sqrt(v_max);
154  double dt = dt_scale * h_min / v_max;
155 
156  double t = 0.0;
157  bool last_step = false;
158  for (int ti = 1; !last_step; ti++)
159  {
160  if (t + dt >= 1.0)
161  {
162  dt = 1.0 - t;
163  last_step = true;
164  }
165  ode_solver.Step(new_field, t, dt);
166  }
167 
168  double glob_minv = minv,
169  glob_maxv = maxv;
170 #ifdef MFEM_USE_MPI
171  if (pfes)
172  {
173  MPI_Allreduce(&minv, &glob_minv, 1, MPI_DOUBLE, MPI_MIN, pfes->GetComm());
174  MPI_Allreduce(&maxv, &glob_maxv, 1, MPI_DOUBLE, MPI_MAX, pfes->GetComm());
175  }
176 #endif
177 
178  // Trim the overshoots and undershoots.
179  new_field.HostReadWrite();
180  for (int i = 0; i < s; i++)
181  {
182  if (new_field(i) < glob_minv) { new_field(i) = glob_minv; }
183  if (new_field(i) > glob_maxv) { new_field(i) = glob_maxv; }
184  }
185 
186  delete oper;
187  delete fess;
188 #ifdef MFEM_USE_MPI
189  delete pfess;
190 #endif
191 }
192 
194  GridFunction &vel,
195  FiniteElementSpace &fes,
196  AssemblyLevel al)
197  : TimeDependentOperator(fes.GetVSize()),
198  x0(x_start), x_now(*fes.GetMesh()->GetNodes()),
199  u(vel), u_coeff(&u), M(&fes), K(&fes), al(al)
200 {
202  K.AddDomainIntegrator(Kinteg);
204  K.Assemble(0);
205  K.Finalize(0);
206 
207  MassIntegrator *Minteg = new MassIntegrator;
208  M.AddDomainIntegrator(Minteg);
210  M.Assemble(0);
211  M.Finalize(0);
212 }
213 
214 void SerialAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
215 {
216  // Move the mesh.
217  const double t = GetTime();
218  add(x0, t, u, x_now);
219  K.FESpace()->GetMesh()->NodesUpdated();
220 
221  // Assemble on the new mesh.
222  K.BilinearForm::operator=(0.0);
223  K.Assemble();
224  Vector rhs(K.Size());
225  K.Mult(ind, rhs);
226  M.BilinearForm::operator=(0.0);
227  M.Assemble();
228 
229  di_dt = 0.0;
230  CGSolver lin_solver;
231  Solver *prec = nullptr;
232  Array<int> ess_tdof_list;
233  if (al == AssemblyLevel::PARTIAL)
234  {
235  prec = new OperatorJacobiSmoother(M, ess_tdof_list);
236  lin_solver.SetOperator(M);
237  }
238  else
239  {
240  prec = new DSmoother(M.SpMat());
241  lin_solver.SetOperator(M.SpMat());
242  }
243  lin_solver.SetPreconditioner(*prec);
244  lin_solver.SetRelTol(1e-12); lin_solver.SetAbsTol(0.0);
245  lin_solver.SetMaxIter(100);
246  lin_solver.SetPrintLevel(0);
247  lin_solver.Mult(rhs, di_dt);
248 
249  delete prec;
250 }
251 
252 #ifdef MFEM_USE_MPI
254  GridFunction &vel,
255  ParFiniteElementSpace &pfes,
256  AssemblyLevel al,
257  MemoryType mt)
258  : TimeDependentOperator(pfes.GetVSize()),
259  x0(x_start), x_now(*pfes.GetMesh()->GetNodes()),
260  u(vel), u_coeff(&u), M(&pfes), K(&pfes), al(al)
261 {
263  if (al == AssemblyLevel::PARTIAL)
264  {
265  Kinteg->SetPAMemoryType(mt);
266  }
267  K.AddDomainIntegrator(Kinteg);
269  K.Assemble(0);
270  K.Finalize(0);
271 
272  MassIntegrator *Minteg = new MassIntegrator;
273  if (al == AssemblyLevel::PARTIAL)
274  {
275  Minteg->SetPAMemoryType(mt);
276  }
277  M.AddDomainIntegrator(Minteg);
279  M.Assemble(0);
280  M.Finalize(0);
281 }
282 
283 void ParAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
284 {
285  // Move the mesh.
286  const double t = GetTime();
287  add(x0, t, u, x_now);
289 
290  // Assemble on the new mesh.
291  K.BilinearForm::operator=(0.0);
292  K.Assemble();
294  K.Mult(ind, rhs);
295  M.BilinearForm::operator=(0.0);
296  M.Assemble();
297 
298  HypreParVector *RHS = rhs.ParallelAssemble();
300  X = 0.0;
301 
302  OperatorHandle Mop;
303  Solver *prec = nullptr;
304  Array<int> ess_tdof_list;
305  if (al == AssemblyLevel::PARTIAL)
306  {
307  M.FormSystemMatrix(ess_tdof_list, Mop);
308  prec = new OperatorJacobiSmoother(M, ess_tdof_list);
309  }
310  else
311  {
312  Mop.Reset(M.ParallelAssemble());
313  prec = new HypreSmoother;
314  static_cast<HypreSmoother*>(prec)->SetType(HypreSmoother::Jacobi, 1);
315  }
316 
317  CGSolver lin_solver(M.ParFESpace()->GetParMesh()->GetComm());
318  lin_solver.SetPreconditioner(*prec);
319  lin_solver.SetOperator(*Mop);
320  lin_solver.SetRelTol(1e-8);
321  lin_solver.SetAbsTol(0.0);
322  lin_solver.SetMaxIter(100);
323  lin_solver.SetPrintLevel(0);
324  lin_solver.Mult(*RHS, X);
325  K.ParFESpace()->GetProlongationMatrix()->Mult(X, di_dt);
326 
327  delete RHS;
328  delete prec;
329 }
330 #endif
331 
332 #ifdef MFEM_USE_GSLIB
334  const Vector &init_field)
335 {
336  nodes0 = init_nodes;
337  Mesh *m = mesh;
338 #ifdef MFEM_USE_MPI
339  if (pmesh) { m = pmesh; }
340 #endif
341  m->SetNodes(nodes0);
342 
343  const double rel_bbox_el = 0.1;
344  const double newton_tol = 1.0e-12;
345  const int npts_at_once = 256;
346 
347  if (finder)
348  {
349  finder->FreeData();
350  delete finder;
351  }
352 
354 #ifdef MFEM_USE_MPI
355  if (pfes)
356  {
357  f = pfes;
358  finder = new FindPointsGSLIB(pfes->GetComm());
359  }
360  else { finder = new FindPointsGSLIB(); }
361 #else
362  finder = new FindPointsGSLIB();
363 #endif
364  finder->Setup(*m, rel_bbox_el, newton_tol, npts_at_once);
365 
366  field0_gf.SetSpace(f);
367  field0_gf = init_field;
368 }
369 
371  Vector &new_field,
372  int new_nodes_ordering)
373 {
374  finder->Interpolate(new_nodes, field0_gf, new_field, new_nodes_ordering);
375 }
376 
377 #endif
378 
380  const Vector &b) const
381 {
382  const FiniteElementSpace *fes = NULL;
383  double energy_in = 0.0;
384 #ifdef MFEM_USE_MPI
385  const ParNonlinearForm *p_nlf = dynamic_cast<const ParNonlinearForm *>(oper);
386  MFEM_VERIFY(!(parallel && p_nlf == NULL), "Invalid Operator subclass.");
387  if (parallel)
388  {
389  fes = p_nlf->FESpace();
390  energy_in = p_nlf->GetEnergy(x);
391  }
392 #endif
393  const bool serial = !parallel;
394  const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
395  MFEM_VERIFY(!(serial && nlf == NULL), "Invalid Operator subclass.");
396  if (serial)
397  {
398  fes = nlf->FESpace();
399  energy_in = nlf->GetEnergy(x);
400  }
401 
402  // Get the local prolongation of the solution vector.
403  Vector x_out_loc(fes->GetVSize(),
405  if (serial)
406  {
407  const SparseMatrix *cP = fes->GetConformingProlongation();
408  if (!cP) { x_out_loc = x; }
409  else { cP->Mult(x, x_out_loc); }
410  }
411 #ifdef MFEM_USE_MPI
412  else
413  {
414  fes->GetProlongationMatrix()->Mult(x, x_out_loc);
415  }
416 #endif
417 
418  double scale = 1.0;
419  if (surf_fit_max_threshold > 0.0)
420  {
421  double avg_err, max_err;
422  GetSurfaceFittingError(avg_err, max_err);
423  if (max_err < surf_fit_max_threshold)
424  {
426  {
427  mfem::out << "TMOPNewtonSolver converged "
428  "based on the surface fitting error.\n";
429  }
430  scale = 0.0;
431  return scale;
432  }
433  }
434 
435  // Check if the starting mesh (given by x) is inverted. Note that x hasn't
436  // been modified by the Newton update yet.
437  const double min_detT_in = ComputeMinDet(x_out_loc, *fes);
438  const bool untangling = (min_detT_in <= 0.0) ? true : false;
439  const double untangle_factor = 1.5;
440  if (untangling)
441  {
442  // Needed for the line search below. The untangling metrics see this
443  // reference to detect deteriorations.
444  MFEM_VERIFY(min_det_ptr != NULL, " Initial mesh was valid, but"
445  " intermediate mesh is invalid. Contact TMOP Developers.");
446  *min_det_ptr = untangle_factor * min_detT_in;
447  }
448 
449  const bool have_b = (b.Size() == Height());
450 
451  Vector x_out(x.Size());
452  bool x_out_ok = false;
453  double energy_out = 0.0, min_detT_out;
454  const double norm_in = Norm(r);
455 
456  const double detJ_factor = (solver_type == 1) ? 0.25 : 0.5;
458  // TODO:
459  // - Customized line search for worst-quality optimization.
460  // - What is the Newton exit criterion for worst-quality optimization?
461 
462  // Perform the line search.
463  for (int i = 0; i < 12; i++)
464  {
465  // Update the mesh and get the L-vector in x_out_loc.
466  add(x, -scale, c, x_out);
467  if (serial)
468  {
469  const SparseMatrix *cP = fes->GetConformingProlongation();
470  if (!cP) { x_out_loc = x_out; }
471  else { cP->Mult(x_out, x_out_loc); }
472  }
473 #ifdef MFEM_USE_MPI
474  else { fes->GetProlongationMatrix()->Mult(x_out, x_out_loc); }
475 #endif
476 
477  // Check the changes in detJ.
478  min_detT_out = ComputeMinDet(x_out_loc, *fes);
479  if (untangling == false && min_detT_out < 0.0)
480  {
481  // No untangling, and detJ got negative -- no good.
483  {
484  mfem::out << "Scale = " << scale << " Neg det(J) found.\n";
485  }
486  scale *= detJ_factor; continue;
487  }
488  if (untangling == true && min_detT_out < *min_det_ptr)
489  {
490  // Untangling, and detJ got even more negative -- no good.
492  {
493  mfem::out << "Scale = " << scale << " Neg det(J) decreased.\n";
494  }
495  scale *= detJ_factor; continue;
496  }
497 
498  // Skip the energy and residual checks when we're untangling. The
499  // untangling metrics change their denominators, which can affect the
500  // energy and residual, so their increase/decrease is not relevant.
501  if (untangling) { x_out_ok = true; break; }
502 
503  // Check the changes in total energy.
504  ProcessNewState(x_out);
505 
506  if (serial)
507  {
508  energy_out = nlf->GetGridFunctionEnergy(x_out_loc);
509  }
510 #ifdef MFEM_USE_MPI
511  else
512  {
513  energy_out = p_nlf->GetParGridFunctionEnergy(x_out_loc);
514  }
515 #endif
516  if (energy_out > energy_in + 0.2*fabs(energy_in) ||
517  std::isnan(energy_out) != 0)
518  {
520  {
521  mfem::out << "Scale = " << scale << " Increasing energy: "
522  << energy_in << " --> " << energy_out << '\n';
523  }
524  scale *= 0.5; continue;
525  }
526 
527  // Check the changes in the Newton residual.
528  oper->Mult(x_out, r);
529  if (have_b) { r -= b; }
530  double norm_out = Norm(r);
531 
532  if (norm_out > 1.2*norm_in)
533  {
535  {
536  mfem::out << "Scale = " << scale << " Norm increased: "
537  << norm_in << " --> " << norm_out << '\n';
538  }
539  scale *= 0.5; continue;
540  }
541  else { x_out_ok = true; break; }
542  } // end line search
543 
544  if (untangling)
545  {
546  // Update the global min detJ. Untangling metrics see this min_det_ptr.
547  if (min_detT_out > 0.0)
548  {
549  *min_det_ptr = 0.0;
552  { mfem::out << "The mesh has been untangled at the used points!\n"; }
553  }
554  else { *min_det_ptr = untangle_factor * min_detT_out; }
555  }
556 
559  {
560  if (untangling)
561  {
562  mfem::out << "Min det(T) change: "
563  << min_detT_in << " -> " << min_detT_out
564  << " with " << scale << " scaling.\n";
565  }
566  else
567  {
568  mfem::out << "Energy decrease: "
569  << energy_in << " --> " << energy_out << " or "
570  << (energy_in - energy_out) / energy_in * 100.0
571  << "% with " << scale << " scaling.\n";
572  }
573  }
574 
575  if (x_out_ok == false) { scale = 0.0; }
576 
579 
580  return scale;
581 }
582 
584 {
585  const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
586  const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
587  TMOP_Integrator *ti = NULL;
588  TMOPComboIntegrator *co = NULL;
589  for (int i = 0; i < integs.Size(); i++)
590  {
591  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
592  if (ti)
593  {
594  ti->UpdateSurfaceFittingWeight(factor);
595  }
596  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
597  if (co)
598  {
600  for (int j = 0; j < ati.Size(); j++)
601  {
602  ati[j]->UpdateSurfaceFittingWeight(factor);
603  }
604  }
605  }
606 }
607 
609 {
610  const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
611  const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
612  TMOP_Integrator *ti = NULL;
613  TMOPComboIntegrator *co = NULL;
614  weights.SetSize(0);
615  double weight;
616 
617  for (int i = 0; i < integs.Size(); i++)
618  {
619  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
620  if (ti)
621  {
622  weight = ti->GetSurfaceFittingWeight();
623  weights.Append(weight);
624  }
625  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
626  if (co)
627  {
629  for (int j = 0; j < ati.Size(); j++)
630  {
631  weight = ati[j]->GetSurfaceFittingWeight();
632  weights.Append(weight);
633  }
634  }
635  }
636 }
637 
639  double &err_max) const
640 {
641  const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
642  const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
643  TMOP_Integrator *ti = NULL;
644  TMOPComboIntegrator *co = NULL;
645 
646  err_avg = 0.0;
647  err_max = 0.0;
648  double err_avg_loc, err_max_loc;
649  for (int i = 0; i < integs.Size(); i++)
650  {
651  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
652  if (ti)
653  {
654  if (ti->IsSurfaceFittingEnabled())
655  {
656  ti->GetSurfaceFittingErrors(err_avg_loc, err_max_loc);
657  err_avg = std::fmax(err_avg_loc, err_avg);
658  err_max = std::fmax(err_max_loc, err_max);
659  }
660  }
661  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
662  if (co)
663  {
665  for (int j = 0; j < ati.Size(); j++)
666  {
667  if (ati[j]->IsSurfaceFittingEnabled())
668  {
669  ati[j]->GetSurfaceFittingErrors(err_avg_loc, err_max_loc);
670  err_avg = std::fmax(err_avg_loc, err_avg);
671  err_max = std::fmax(err_max_loc, err_max);
672  }
673  }
674  }
675  }
676 }
677 
679 {
680  const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
681  const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
682 
683  // Reset the update flags of all TargetConstructors. This is done to avoid
684  // repeated updates of shared TargetConstructors.
685  TMOP_Integrator *ti = NULL;
686  TMOPComboIntegrator *co = NULL;
687  DiscreteAdaptTC *dtc = NULL;
688  for (int i = 0; i < integs.Size(); i++)
689  {
690  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
691  if (ti)
692  {
693  dtc = ti->GetDiscreteAdaptTC();
694  if (dtc) { dtc->ResetUpdateFlags(); }
695  }
696  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
697  if (co)
698  {
700  for (int j = 0; j < ati.Size(); j++)
701  {
702  dtc = ati[j]->GetDiscreteAdaptTC();
703  if (dtc) { dtc->ResetUpdateFlags(); }
704  }
705  }
706  }
707 
708  if (parallel)
709  {
710 #ifdef MFEM_USE_MPI
711  const ParNonlinearForm *pnlf =
712  dynamic_cast<const ParNonlinearForm *>(oper);
713  const ParFiniteElementSpace *pfesc = pnlf->ParFESpace();
714  Vector x_loc(pfesc->GetVSize());
715  pfesc->GetProlongationMatrix()->Mult(x, x_loc);
716  for (int i = 0; i < integs.Size(); i++)
717  {
718  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
719  if (ti)
720  {
721  ti->UpdateAfterMeshPositionChange(x_loc, pfesc->GetOrdering());
722  ti->ComputeFDh(x_loc, *pfesc);
724  {
725  ti->ComputeUntangleMetricQuantiles(x_loc, *pfesc);
726  }
727  UpdateDiscreteTC(*ti, x_loc, pfesc->GetOrdering());
728  }
729  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
730  if (co)
731  {
733  for (int j = 0; j < ati.Size(); j++)
734  {
735  ati[j]->UpdateAfterMeshPositionChange(x_loc, pfesc->GetOrdering());
736  ati[j]->ComputeFDh(x_loc, *pfesc);
738  {
739  ati[j]->ComputeUntangleMetricQuantiles(x_loc, *pfesc);
740  }
741  UpdateDiscreteTC(*ati[j], x_loc, pfesc->GetOrdering());
742  }
743  }
744  }
745 #endif
746  }
747  else
748  {
749  const FiniteElementSpace *fesc = nlf->FESpace();
750  const Operator *P = nlf->GetProlongation();
751  Vector x_loc;
752  if (P)
753  {
754  x_loc.SetSize(P->Height());
755  P->Mult(x,x_loc);
756  }
757  else
758  {
759  x_loc = x;
760  }
761  for (int i = 0; i < integs.Size(); i++)
762  {
763  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
764  if (ti)
765  {
766  ti->UpdateAfterMeshPositionChange(x_loc, fesc->GetOrdering());
767  ti->ComputeFDh(x_loc, *fesc);
769  {
770  ti->ComputeUntangleMetricQuantiles(x_loc, *fesc);
771  }
772  UpdateDiscreteTC(*ti, x_loc, fesc->GetOrdering());
773  }
774  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
775  if (co)
776  {
778  for (int j = 0; j < ati.Size(); j++)
779  {
780  ati[j]->UpdateAfterMeshPositionChange(x_loc, fesc->GetOrdering());
781  ati[j]->ComputeFDh(x_loc, *fesc);
783  {
784  ati[j]->ComputeUntangleMetricQuantiles(x_loc, *fesc);
785  }
786  UpdateDiscreteTC(*ati[j], x_loc, fesc->GetOrdering());
787  }
788  }
789  }
790  }
791 
792  // Constant coefficient associated with the surface fitting terms if
793  // adaptive surface fitting is enabled. The idea is to increase the
794  // coefficient if the surface fitting error does not sufficiently
795  // decrease between subsequent TMOPNewtonSolver iterations.
797  {
798  double surf_fit_err_max = -10;
799  double surf_fit_err_avg = -10;
800  // Get surface fitting errors.
801  GetSurfaceFittingError(surf_fit_err_avg, surf_fit_err_max);
802  // Get array with surface fitting weights.
803  Array<double> weights;
804  GetSurfaceFittingWeight(weights);
805 
807  {
808  mfem::out << "Avg/Max surface fitting error: " <<
809  surf_fit_err_avg << " " <<
810  surf_fit_err_max << "\n";
811  mfem::out << "Min/Max surface fitting weight: " <<
812  weights.Min() << " " << weights.Max() << "\n";
813  }
814 
815  double change_surf_fit_err = surf_fit_err_avg_prvs-surf_fit_err_avg;
816  double rel_change_surf_fit_err = change_surf_fit_err/surf_fit_err_avg_prvs;
817  // Increase the surface fitting coefficient if the surface fitting error
818  // does not decrease sufficiently.
819  if (rel_change_surf_fit_err < 1.e-2)
820  {
822  }
823  surf_fit_err_avg_prvs = surf_fit_err_avg;
824  update_surf_fit_coeff = false;
825  }
826 }
827 
829  const Vector &x_new,
830  int x_ordering) const
831 {
832  const bool update_flag = true;
833  DiscreteAdaptTC *discrtc = ti.GetDiscreteAdaptTC();
834  if (discrtc)
835  {
836  discrtc->UpdateTargetSpecification(x_new, update_flag, x_ordering);
837  if (ti.GetFDFlag())
838  {
839  double dx = ti.GetFDh();
840  discrtc->UpdateGradientTargetSpecification(x_new, dx, update_flag, x_ordering);
841  discrtc->UpdateHessianTargetSpecification(x_new, dx, update_flag, x_ordering);
842  }
843  }
844 }
845 
847  const FiniteElementSpace &fes) const
848 {
849  double min_detJ = infinity();
850  const int NE = fes.GetNE(), dim = fes.GetMesh()->Dimension();
851  Array<int> xdofs;
852  DenseMatrix Jpr(dim);
853  const bool mixed_mesh = fes.GetMesh()->GetNumGeometries(dim) > 1;
854  if (dim == 1 || mixed_mesh || UsesTensorBasis(fes) == false)
855  {
856  for (int i = 0; i < NE; i++)
857  {
858  const int dof = fes.GetFE(i)->GetDof();
859  DenseMatrix dshape(dof, dim), pos(dof, dim);
860  Vector posV(pos.Data(), dof * dim);
861 
862  fes.GetElementVDofs(i, xdofs);
863  x_loc.GetSubVector(xdofs, posV);
864 
865  const IntegrationRule &irule = GetIntegrationRule(*fes.GetFE(i));
866  const int nsp = irule.GetNPoints();
867  for (int j = 0; j < nsp; j++)
868  {
869  fes.GetFE(i)->CalcDShape(irule.IntPoint(j), dshape);
870  MultAtB(pos, dshape, Jpr);
871  min_detJ = std::min(min_detJ, Jpr.Det());
872  }
873  }
874  }
875  else
876  {
877  min_detJ = dim == 2 ? MinDetJpr_2D(&fes, x_loc) :
878  dim == 3 ? MinDetJpr_3D(&fes, x_loc) : 0.0;
879  }
880  double min_detT_all = min_detJ;
881 #ifdef MFEM_USE_MPI
882  if (parallel)
883  {
884  auto p_nlf = dynamic_cast<const ParNonlinearForm *>(oper);
885  MPI_Allreduce(&min_detJ, &min_detT_all, 1, MPI_DOUBLE, MPI_MIN,
886  p_nlf->ParFESpace()->GetComm());
887  }
888 #endif
889  const DenseMatrix &Wideal =
891  min_detT_all /= Wideal.Det();
892 
893  return min_detT_all;
894 }
895 
896 #ifdef MFEM_USE_MPI
897 // Metric values are visualized by creating an L2 finite element functions and
898 // computing the metric values at the nodes.
900  const TargetConstructor &tc, ParMesh &pmesh,
901  char *title, int position)
902 {
904  ParFiniteElementSpace fes(&pmesh, &fec, 1);
905  ParGridFunction metric(&fes);
906  InterpolateTMOP_QualityMetric(qm, tc, pmesh, metric);
907  socketstream sock;
908  if (pmesh.GetMyRank() == 0)
909  {
910  sock.open("localhost", 19916);
911  sock << "solution\n";
912  }
913  pmesh.PrintAsOne(sock);
914  metric.SaveAsOne(sock);
915  if (pmesh.GetMyRank() == 0)
916  {
917  sock << "window_title '"<< title << "'\n"
918  << "window_geometry "
919  << position << " " << 0 << " " << 600 << " " << 600 << "\n"
920  << "keys jRmclA\n";
921  }
922 }
923 #endif
924 
925 // Metric values are visualized by creating an L2 finite element functions and
926 // computing the metric values at the nodes.
928  const TargetConstructor &tc, Mesh &mesh,
929  char *title, int position)
930 {
932  FiniteElementSpace fes(&mesh, &fec, 1);
933  GridFunction metric(&fes);
934  InterpolateTMOP_QualityMetric(qm, tc, mesh, metric);
935  osockstream sock(19916, "localhost");
936  sock << "solution\n";
937  mesh.Print(sock);
938  metric.Save(sock);
939  sock.send();
940  sock << "window_title '"<< title << "'\n"
941  << "window_geometry "
942  << position << " " << 0 << " " << 600 << " " << 600 << "\n"
943  << "keys jRmclA\n";
944 }
945 
946 }
VectorGridFunctionCoefficient u_coeff
Definition: tmop_tools.hpp:91
void vis_tmop_metric_p(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, ParMesh &pmesh, char *title, int position)
Definition: tmop_tools.cpp:899
void ComputeUntangleMetricQuantiles(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:4238
virtual void Interpolate(const GridFunction &field_in, Vector &field_out)
Definition: gslib.cpp:803
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
Definition: array.cpp:85
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
Conjugate gradient method.
Definition: solvers.hpp:493
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
AssemblyLevel
Enumeration defining the assembly level for bilinear and nonlinear form classes derived from Operator...
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1153
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
const Operator * oper
Definition: solvers.hpp:120
Data type for scaled Jacobi-type smoother of sparse matrix.
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field, int new_nodes_ordering=Ordering::byNODES)
Definition: tmop_tools.cpp:370
double GetSurfaceFittingWeight()
Get the surface fitting weight.
Definition: tmop.cpp:3915
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void Mult(const Vector &ind, Vector &di_dt) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: tmop_tools.cpp:214
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
void GetSurfaceFittingWeight(Array< double > &weights) const
Get the surface fitting weight for all the TMOP integrators.
Definition: tmop_tools.cpp:608
FiniteElementSpace * FESpace()
Return the FE space associated with the BilinearForm.
int Dimension() const
Definition: mesh.hpp:1047
double ComputeMinDet(const Vector &x_loc, const FiniteElementSpace &fes) const
Definition: tmop_tools.cpp:846
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:899
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric, const TargetConstructor &tc, const Mesh &mesh, GridFunction &metric_gf)
Interpolates the metric&#39;s values at the nodes of metric_gf.
Definition: tmop.cpp:4460
Base abstract class for first order time dependent operators.
Definition: operator.hpp:316
bool summary
A summary of the solver process will be reported after the last iteration to mfem::out.
Definition: solvers.hpp:91
const Array< TMOP_Integrator * > & GetTMOPIntegrators() const
Definition: tmop.hpp:2083
bool first_and_last
Information about the first and last iteration will be printed to mfem::out.
Definition: solvers.hpp:94
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
DiscreteAdaptTC * GetDiscreteAdaptTC() const
Definition: tmop.hpp:2031
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:712
double Det() const
Definition: densemat.cpp:488
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
ParFiniteElementSpace * pfes
Definition: tmop.hpp:1190
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void GetSurfaceFittingErrors(double &err_avg, double &err_max)
Definition: tmop.cpp:2882
double MinDetJpr_2D(const FiniteElementSpace *, const Vector &) const
Definition: tmop_pa_jp2.cpp:75
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
Parallel non-linear operator on the true dofs.
void Init(TimeDependentOperator &f_) override
Associate a TimeDependentOperator with the ODE solver.
Definition: ode.cpp:100
bool IsSurfaceFittingEnabled()
Definition: tmop.hpp:1973
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void vis_tmop_metric_s(int order, TMOP_QualityMetric &qm, const TargetConstructor &tc, Mesh &mesh, char *title, int position)
Definition: tmop_tools.cpp:927
void UpdateGradientTargetSpecification(const Vector &x, double dx, bool use_flag=false, int x_ordering=Ordering::byNODES)
Definition: tmop.cpp:2607
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:548
void UpdateSurfaceFittingWeight(double factor)
Update the surface fitting weight as surf_fit_coeff *= factor;.
Definition: tmop.cpp:3903
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition: fespace.hpp:957
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication: .
Geometry Geometries
Definition: fe.cpp:49
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:314
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:5921
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field, int new_nodes_ordering=Ordering::byNODES)
Definition: tmop_tools.cpp:29
double t
Current time.
Definition: operator.hpp:340
double Norm(const Vector &x) const
Definition: solvers.hpp:172
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:76
Data type sparse matrix.
Definition: sparsemat.hpp:50
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:756
virtual void FreeData()
Definition: gslib.cpp:267
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
void UpdateAfterMeshPositionChange(const Vector &new_x, int new_x_ordering=Ordering::byNODES)
Definition: tmop.cpp:4054
Array< NonlinearFormIntegrator * > * GetDNFI()
Access all integrators added with AddDomainIntegrator().
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:302
double b
Definition: lissajous.cpp:42
const AssemblyLevel al
Definition: tmop_tools.hpp:93
static MemoryType GetDeviceMemoryType()
Get the current Device MemoryType. This is the MemoryType used by most MFEM classes when allocating m...
Definition: device.hpp:273
virtual double ComputeScalingFactor(const Vector &x, const Vector &b) const
Definition: tmop_tools.cpp:379
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
void Assemble(int skip_zeros=1)
Assemble the local matrix.
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:107
Performs a single remap advection step in serial.
Definition: tmop_tools.hpp:85
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1232
virtual double GetEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
virtual void Mult(const Vector &ind, Vector &di_dt) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
Definition: tmop_tools.cpp:283
Parallel smoothers in hypre.
Definition: hypre.hpp:971
double GetGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
void SetPAMemoryType(MemoryType mt)
Definition: nonlininteg.hpp:50
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
double GetParGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
FiniteElementSpace * FESpace()
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
int GetMyRank() const
Definition: pmesh.hpp:353
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
void SetAbsTol(double atol)
Definition: solvers.hpp:200
virtual void ProcessNewState(const Vector &x) const
Definition: tmop_tools.cpp:678
PrintLevel print_options
Output behavior for the iterative solver.
Definition: solvers.hpp:137
void SetRelTol(double rtol)
Definition: solvers.hpp:199
double GetFDh() const
Definition: tmop.hpp:2048
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
ParFiniteElementSpace * ParFESpace() const
const AssemblyLevel al
Definition: tmop_tools.hpp:115
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
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:461
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1194
Performs a single remap advection step in parallel.
Definition: tmop_tools.hpp:107
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
bool GetFDFlag() const
Definition: tmop.hpp:2047
string space
double MinDetJpr_3D(const FiniteElementSpace *, const Vector &) const
Definition: tmop_pa_jp3.cpp:77
void PrintAsOne(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4983
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:894
void ResetUpdateFlags()
Used in combination with the Update methods to avoid extra computations.
Definition: tmop.hpp:1515
const IntegrationRule & GetIntegrationRule(const FiniteElement &el) const
Definition: tmop_tools.hpp:157
virtual double GetTime() const
Read the currently set time.
Definition: operator.hpp:357
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
virtual void FormSystemMatrix(const Array< int > &ess_tdof_list, OperatorHandle &A)
Form the linear system matrix A, see FormLinearSystem() for details.
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:479
void ComputeFDh(const Vector &x, const FiniteElementSpace &fes)
Determines the perturbation, h, for FD-based approximation.
Definition: tmop.cpp:4074
void UpdateHessianTargetSpecification(const Vector &x, double dx, bool use_flag=false, int x_ordering=Ordering::byNODES)
Definition: tmop.cpp:2642
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
int dim
Definition: ex24.cpp:53
virtual void SetInitialField(const Vector &init_nodes, const Vector &init_field)
Definition: tmop_tools.cpp:22
VectorGridFunctionCoefficient u_coeff
Definition: tmop_tools.hpp:113
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
SerialAdvectorCGOper(const Vector &x_start, GridFunction &vel, FiniteElementSpace &fes, AssemblyLevel al=AssemblyLevel::LEGACY)
Definition: tmop_tools.cpp:193
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:46
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:196
int open(const char hostname[], int port)
Open the socket stream on &#39;port&#39; at &#39;hostname&#39;.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:507
virtual void GetSurfaceFittingError(double &err_avg, double &err_max) const
Definition: tmop_tools.cpp:638
ParAdvectorCGOper(const Vector &x_start, GridFunction &vel, ParFiniteElementSpace &pfes, AssemblyLevel al=AssemblyLevel::LEGACY, MemoryType mt=MemoryType::DEFAULT)
Definition: tmop_tools.cpp:253
FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary collection of points...
Definition: gslib.hpp:52
Vector data type.
Definition: vector.hpp:60
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1749
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:576
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
void Step(Vector &x, double &t, double &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:109
void UpdateTargetSpecification(const Vector &new_x, bool use_flag=false, int new_x_ordering=Ordering::byNODES)
Definition: tmop.cpp:1960
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7949
virtual const Operator * GetProlongation() const
Get the finite element space prolongation matrix.
RefCoord s[3]
int Size() const
Get the size of the BilinearForm as a square matrix.
bool iterations
Detailed information about each iteration will be reported to mfem::out.
Definition: solvers.hpp:88
double GetEnergy(const ParGridFunction &x) const
Compute the energy of a ParGridFunction.
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
Definition: mesh.cpp:7961
Base class for solvers.
Definition: operator.hpp:682
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1234
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3673
Abstract operator.
Definition: operator.hpp:24
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...
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:468
void UpdateSurfaceFittingWeight(double factor) const
Update surface fitting weight as surf_fit_weight *= factor.
Definition: tmop_tools.cpp:583
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
virtual void SetInitialField(const Vector &init_nodes, const Vector &init_field)
Definition: tmop_tools.cpp:333
FiniteElementSpace * fes
Definition: tmop.hpp:1185
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
void NodesUpdated()
This function should be called after the mesh node coordinates have been updated externally, e.g. by modifying the internal nodal GridFunction returned by GetNodes().
Definition: mesh.hpp:1038
alpha (q . grad u, v)
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1638
void UpdateDiscreteTC(const TMOP_Integrator &ti, const Vector &x_new, int x_ordering=Ordering::byNODES) const
Definition: tmop_tools.cpp:828