MFEM  v4.6.0
Finite element discretization library
tmop_amr.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_amr.hpp"
13 
14 namespace mfem
15 {
16 
17 using namespace mfem;
18 
20 {
21  bool iso = false;
22  bool aniso = false;
23  if (amrmetric == 1 || amrmetric == 2 || amrmetric == 58)
24  {
25  aniso = true;
26  }
27  if (amrmetric == 55 || amrmetric == 56 || amrmetric == 77 ||
28  amrmetric == 315 || amrmetric == 316 || amrmetric == 321)
29  {
30  iso = true;
31  }
32  if (amrmetric == 7 || amrmetric == 9)
33  {
34  iso = true; aniso = true;
35  }
36 
37  MFEM_VERIFY(iso || aniso, "Metric type not supported in hr-adaptivity.");
38 
39  const int dim = mesh->Dimension();
40  const int num_ref_types = 3 + 4*(dim-2);
41  const int NEorig = mesh->GetNE();
42 
43  aniso_flags.SetSize(NEorig);
44  error_estimates.SetSize(NEorig);
45  Vector amr_base_energy(NEorig), amr_temp_energy(NEorig);
46  error_estimates = 1.*std::numeric_limits<float>::max();
47  aniso_flags = -1;
48  GetTMOPRefinementEnergy(0, amr_base_energy);
49 
50  for (int i = 1; i < num_ref_types+1; i++)
51  {
52  if ( dim == 2 && i < 3 && aniso != true ) { continue; }
53  if ( dim == 2 && i == 3 && iso != true ) { continue; }
54  if ( dim == 3 && i < 7 && aniso != true ) { continue; }
55  if ( dim == 3 && i == 7 && iso != true ) { continue; }
56 
57  GetTMOPRefinementEnergy(i, amr_temp_energy);
58 
59  for (int e = 0; e < NEorig; e++)
60  {
61  if ( amr_temp_energy(e) < error_estimates(e) )
62  {
63  error_estimates(e) = amr_temp_energy(e);
64  aniso_flags[e] = i;
65  }
66  }
67  }
69 
70  if (spat_gf)
71  {
72  L2_FECollection avg_fec(0, mesh->Dimension());
73  FiniteElementSpace avg_fes(spat_gf->FESpace()->GetMesh(), &avg_fec);
74  GridFunction elem_avg(&avg_fes);
75  spat_gf->GetElementAverages(elem_avg);
76  for (int i = 0; i < amr_base_energy.Size(); i++)
77  {
78  if (elem_avg(i) < spat_gf_critical) { amr_base_energy(i) = 0.; }
79  }
80  }
81 
82  error_estimates -= amr_base_energy;
83  error_estimates *= -1; // error = E(parent) - scaling_factor*mean(E(children))
85 }
86 
88  Vector &el_energy_vec)
89 {
90  const FiniteElementSpace *fes = mesh->GetNodalFESpace();
91  const int NE = fes->GetNE();
92  GridFunction *xdof = mesh->GetNodes();
93  xdof->SetTrueVector();
94  xdof->SetFromTrueVector();
95 
96  el_energy_vec.SetSize(NE);
97  el_energy_vec = std::numeric_limits<float>::max();
98 
99  for (int e = 0; e < NE; e++)
100  {
101  Geometry::Type gtype = fes->GetFE(e)->GetGeomType();
102  DenseMatrix tr, xsplit;
103  IntegrationRule *irule = NULL;
104 
105  if ( (gtype == Geometry::TRIANGLE && reftype > 0 && reftype < 3) ||
106  (gtype == Geometry::CUBE && reftype > 0 && reftype < 7) ||
107  (gtype == Geometry::TETRAHEDRON && reftype > 0 && reftype < 7) )
108  {
109  continue;
110  }
111 
112  switch (gtype)
113  {
114  case Geometry::TRIANGLE:
115  {
116  int ref_access = reftype == 0 ? 0 : 1;
117  xdof->GetVectorValues(e, *TriIntRule[ref_access], xsplit, tr);
118  irule = TriIntRule[ref_access];
119  break;
120  }
122  {
123  int ref_access = reftype == 0 ? 0 : 1;
124  xdof->GetVectorValues(e, *TetIntRule[ref_access], xsplit, tr);
125  irule = TetIntRule[ref_access];
126  break;
127  }
128  case Geometry::SQUARE:
129  {
130  MFEM_VERIFY(QuadIntRule[reftype], " Integration rule does not exist.");
131  xdof->GetVectorValues(e, *QuadIntRule[reftype], xsplit, tr);
132  irule = QuadIntRule[reftype];
133  break;
134  }
135  case Geometry::CUBE:
136  {
137  int ref_access = reftype == 0 ? 0 : 1;
138  xdof->GetVectorValues(e, *HexIntRule[ref_access], xsplit, tr);
139  irule = HexIntRule[ref_access];
140  break;
141  }
142  default:
143  MFEM_ABORT("Incompatible geometry type!");
144  }
145  xsplit.Transpose();
146 
147  el_energy_vec(e) = 0.; // Re-set to 0
148 
149  // The data format is xe1,xe2,..xen,ye1,ye2..yen.
150  // We will reformat it inside GetRefinementElementEnergy
151  Vector elfun(xsplit.GetData(), xsplit.NumCols()*xsplit.NumRows());
152 
154  TMOP_Integrator *ti = NULL;
155  TMOPComboIntegrator *co = NULL;
156  for (int i = 0; i < integs.Size(); i++)
157  {
158  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
159  if (ti)
160  {
161  el_energy_vec(e) = ti->GetRefinementElementEnergy(*fes->GetFE(e),
163  elfun,
164  *irule);
165  }
166  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
167  if (co)
168  {
170  for (int j = 0; j < ati.Size(); j++)
171  {
172  el_energy_vec(e) += ati[j]->GetRefinementElementEnergy(*fes->GetFE(e),
174  elfun,
175  *irule);
176  }
177  }
178  }
179  }
180 }
181 
183 {
184  HexIntRule.SetSize(1+1);
185  // Reftype = 0 -> original element
186  Mesh meshsplit = Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON);
187  Mesh base_mesh_copy(meshsplit);
188  HexIntRule[0] = SetIntRulesFromMesh(meshsplit);
189  meshsplit.Clear();
190 
191  // Reftype = 7
192  for (int i = 7; i < 8; i++)
193  {
194  Array<Refinement> marked_elements;
195  Mesh mesh_ref(base_mesh_copy);
196  for (int e = 0; e < mesh_ref.GetNE(); e++)
197  {
198  marked_elements.Append(Refinement(e, i));
199  }
200  mesh_ref.GeneralRefinement(marked_elements, 1, 0);
201  HexIntRule[1] = SetIntRulesFromMesh(mesh_ref);
202  mesh_ref.Clear();
203  }
204 }
205 
207 {
208  QuadIntRule.SetSize(3+1);
209 
210  // Reftype = 0 -> original element
212  Mesh base_mesh_copy(meshsplit);
213  QuadIntRule[0] = SetIntRulesFromMesh(meshsplit);
214  meshsplit.Clear();
215 
216  // Reftype = 1-3
217  for (int i = 1; i < 4; i++)
218  {
219  Array<Refinement> marked_elements;
220  Mesh mesh_ref(base_mesh_copy);
221  for (int e = 0; e < mesh_ref.GetNE(); e++)
222  {
223  marked_elements.Append(Refinement(e, i));
224  }
225  mesh_ref.GeneralRefinement(marked_elements, 1, 0);
226  QuadIntRule[i] = SetIntRulesFromMesh(mesh_ref);
227  mesh_ref.Clear();
228  }
229 }
230 
232 {
233  TriIntRule.SetSize(1+1);
234 
235  // Reftype = 0 // original element
236  const int Nvert = 3, NEsplit = 1;
237  Mesh meshsplit(2, Nvert, NEsplit, 0, 2);
238  const double tri_v[3][2] =
239  {
240  {0, 0}, {1, 0}, {0, 1}
241  };
242  const int tri_e[1][3] =
243  {
244  {0, 1, 2}
245  };
246 
247  for (int j = 0; j < Nvert; j++)
248  {
249  meshsplit.AddVertex(tri_v[j]);
250  }
251  meshsplit.AddTriangle(tri_e[0], 1);
252  meshsplit.FinalizeTriMesh(1, 1, true);
253 
254  Mesh base_mesh_copy(meshsplit);
255  TriIntRule[0] = SetIntRulesFromMesh(meshsplit);
256  meshsplit.Clear();
257 
258  // no anisotropic refinements for triangle
259  // Reftype = 3
260  for (int i = 1; i < 2; i++)
261  {
262  Array<Refinement> marked_elements;
263  Mesh mesh_ref(base_mesh_copy);
264  for (int e = 0; e < mesh_ref.GetNE(); e++)
265  {
266  marked_elements.Append(Refinement(e, i));
267  }
268  mesh_ref.GeneralRefinement(marked_elements, 1, 0);
269  TriIntRule[i] = SetIntRulesFromMesh(mesh_ref);
270  mesh_ref.Clear();
271  }
272 }
273 
275 {
276  TetIntRule.SetSize(1+1);
277 
278  // Reftype = 0 // original element
279  const int Nvert = 4, NEsplit = 1;
280  Mesh meshsplit(3, Nvert, NEsplit, 0, 3);
281  const double tet_v[4][3] =
282  {
283  {0, 0, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}
284  };
285  const int tet_e[1][4] =
286  {
287  {0, 1, 2, 3}
288  };
289 
290  for (int j = 0; j < Nvert; j++)
291  {
292  meshsplit.AddVertex(tet_v[j]);
293  }
294  meshsplit.AddTet(tet_e[0], 1);
295  meshsplit.FinalizeTetMesh(1, 1, true);
296 
297  Mesh base_mesh_copy(meshsplit);
298  TetIntRule[0] = SetIntRulesFromMesh(meshsplit);
299  meshsplit.Clear();
300 
301  // no anisotropic refinements for triangle
302  // Reftype = 7
303  for (int i = 1; i < 2; i++)
304  {
305  Array<Refinement> marked_elements;
306  Mesh mesh_ref(base_mesh_copy);
307  for (int e = 0; e < mesh_ref.GetNE(); e++)
308  {
309  marked_elements.Append(Refinement(e, i)); // ref_type will default to 7
310  }
311  mesh_ref.GeneralRefinement(marked_elements, 1, 0);
312  TetIntRule[i] = SetIntRulesFromMesh(mesh_ref);
313  mesh_ref.Clear();
314  }
315 }
316 
318 {
319  const int dim = meshsplit.Dimension();
320  H1_FECollection fec(order, dim);
321  FiniteElementSpace nodal_fes(&meshsplit, &fec, dim);
322  meshsplit.SetNodalFESpace(&nodal_fes);
323 
324  const int NEsplit = meshsplit.GetNE();
325  const int dof_cnt = nodal_fes.GetFE(0)->GetDof(),
326  pts_cnt = NEsplit * dof_cnt;
327 
328  DenseMatrix pos(dof_cnt, dim);
329  Vector posV(pos.Data(), dof_cnt * dim);
330  Array<int> xdofs(dof_cnt * dim);
331 
332  // Create an IntegrationRule on the nodes of the reference submesh.
333  IntegrationRule *irule = new IntegrationRule(pts_cnt);
334  GridFunction *nodesplit = meshsplit.GetNodes();
335 
336  int pt_id = 0;
337  for (int i = 0; i < NEsplit; i++)
338  {
339  nodal_fes.GetElementVDofs(i, xdofs);
340  nodesplit->GetSubVector(xdofs, posV);
341  for (int j = 0; j < dof_cnt; j++)
342  {
343  if (dim == 2)
344  {
345  irule->IntPoint(pt_id).Set2(pos(j, 0), pos(j, 1));
346  }
347  else if (dim == 3)
348  {
349  irule->IntPoint(pt_id).Set3(pos(j, 0), pos(j, 1), pos(j, 2));
350  }
351  pt_id++;
352  }
353  }
354  return irule;
355 }
356 
358  TMOP_Integrator &tmopi,
359  Vector &fine_energy)
360 {
361  DiscreteAdaptTC *tcd = tmopi.GetDiscreteAdaptTC();
362  fine_energy.SetSize(mesh->GetNE());
363 
364  if (serial)
365  {
366  Mesh meshcopy(*mesh);
367  FiniteElementSpace *tcdfes = NULL;
368  if (tcd)
369  {
370  tcdfes = new FiniteElementSpace(*tcd->GetTSpecFESpace(), &meshcopy);
371  }
372 
373  Vector local_err(meshcopy.GetNE());
374  local_err = 0.;
375  double threshold = std::numeric_limits<float>::max();
376  meshcopy.DerefineByError(local_err, threshold, 0, 1);
377 
378  if (meshcopy.GetGlobalNE() == mesh->GetGlobalNE())
379  {
380  delete tcdfes;
381  return false;
382  }
383 
384  if (tcd)
385  {
386  tcdfes->Update();
387  tcd->SetTspecDataForDerefinement(tcdfes);
388  }
389 
390  Vector coarse_energy(meshcopy.GetNE());
391  GetTMOPDerefinementEnergy(meshcopy, tmopi, coarse_energy);
392  if (tcd) { tcd->ResetDerefinementTspecData(); }
393  GetTMOPDerefinementEnergy(*mesh, tmopi, fine_energy);
394 
395  const CoarseFineTransformations &dtrans =
396  meshcopy.ncmesh->GetDerefinementTransforms();
397 
398  Table coarse_to_fine;
399  dtrans.MakeCoarseToFineTable(coarse_to_fine);
400 
401  Array<int> tabrow;
402  for (int pe = 0; pe < coarse_to_fine.Size(); pe++)
403  {
404  coarse_to_fine.GetRow(pe, tabrow);
405  int nchild = tabrow.Size();
406  double parent_energy = coarse_energy(pe);
407  for (int fe = 0; fe < nchild; fe++)
408  {
409  int child = tabrow[fe];
410  MFEM_VERIFY(child < mesh->GetNE(), " invalid coarse to fine mapping");
411  fine_energy(child) -= parent_energy;
412  }
413  }
414  delete tcdfes;
415  }
416  else
417  {
418 #ifdef MFEM_USE_MPI
419  ParMesh meshcopy(*pmesh);
420  ParFiniteElementSpace *tcdfes = NULL;
421  if (tcd)
422  {
423  tcdfes = new ParFiniteElementSpace(*tcd->GetTSpecParFESpace(), meshcopy);
424  }
425 
426  Vector local_err(meshcopy.GetNE());
427  local_err = 0.;
428  double threshold = std::numeric_limits<float>::max();
429  meshcopy.DerefineByError(local_err, threshold, 0, 1);
430 
431  if (meshcopy.GetGlobalNE() == pmesh->GetGlobalNE())
432  {
433  delete tcdfes;
434  return false;
435  }
436 
437  if (tcd)
438  {
439  tcdfes->Update();
440  tcd->SetTspecDataForDerefinement(tcdfes);
441  }
442 
443  Vector coarse_energy(meshcopy.GetNE());
444  GetTMOPDerefinementEnergy(meshcopy, tmopi, coarse_energy);
445  if (tcd) { tcd->ResetDerefinementTspecData(); }
446  GetTMOPDerefinementEnergy(*pmesh, tmopi, fine_energy);
447 
448  const CoarseFineTransformations &dtrans =
450 
451  Table coarse_to_fine;
452  dtrans.MakeCoarseToFineTable(coarse_to_fine);
453 
454  Array<int> tabrow;
455  for (int pe = 0; pe < meshcopy.GetNE(); pe++)
456  {
457  coarse_to_fine.GetRow(pe, tabrow);
458  int nchild = tabrow.Size();
459  double parent_energy = coarse_energy(pe);
460  for (int fe = 0; fe < nchild; fe++)
461  {
462  int child = tabrow[fe];
463  MFEM_VERIFY(child < pmesh->GetNE(), " invalid coarse to fine mapping");
464  fine_energy(child) -= parent_energy;
465  }
466  }
467  delete tcdfes;
468 #endif
469  }
470 
471  // error_estimate(e) = energy(parent_of_e)-energy(e)
472  // Negative energy means derefinement is desirable.
473  fine_energy *= -1;
474  return true;
475 }
476 
478 {
480  TMOP_Integrator *ti = NULL;
481  TMOPComboIntegrator *co = NULL;
483  error_estimates = 0.;
484  Vector fine_energy(mesh->GetNE());
485 
486  for (int i = 0; i < integs.Size(); i++)
487  {
488  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
489  if (ti)
490  {
491  bool deref = GetDerefineEnergyForIntegrator(*ti, fine_energy);
492  if (!deref) { error_estimates = 1; return; }
493  error_estimates += fine_energy;
494  }
495  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
496  if (co)
497  {
499  for (int j = 0; j < ati.Size(); j++)
500  {
501  bool deref = GetDerefineEnergyForIntegrator(*ati[j], fine_energy);
502  if (!deref) { error_estimates = 1; return; }
503  error_estimates += fine_energy;
504  }
505  }
506  }
507 }
508 
510  TMOP_Integrator &tmopi,
511  Vector &el_energy_vec)
512 {
513  const int cNE = cmesh.GetNE();
514  el_energy_vec.SetSize(cNE);
515  const FiniteElementSpace *fespace = cmesh.GetNodalFESpace();
516 
517  GridFunction *cxdof = cmesh.GetNodes();
518 
519  Array<int> vdofs;
520  Vector el_x;
521  const FiniteElement *fe;
523 
524  for (int j = 0; j < cNE; j++)
525  {
526  fe = fespace->GetFE(j);
527  fespace->GetElementVDofs(j, vdofs);
528  T = cmesh.GetElementTransformation(j);
529  cxdof->GetSubVector(vdofs, el_x);
530  el_energy_vec(j) = tmopi.GetDerefinementElementEnergy(*fe, *T, el_x);
531  }
532 }
533 
534 
536  TMOPNewtonSolver &tmopns_, GridFunction &x_,
537  bool move_bnd_, bool hradaptivity_,
538  int mesh_poly_deg_, int amr_metric_id_,
539  int hr_iter_, int h_per_r_iter_) :
540  mesh(&mesh_), nlf(&nlf_), tmopns(&tmopns_), x(&x_),
541  gridfuncarr(), fespacearr(),
542  move_bnd(move_bnd_), hradaptivity(hradaptivity_),
543  mesh_poly_deg(mesh_poly_deg_), amr_metric_id(amr_metric_id_),
544  serial(true), hr_iter(hr_iter_), h_per_r_iter(h_per_r_iter_)
545 {
546  if (!hradaptivity) { return; }
548  amr_metric_id);
555 }
556 
557 #ifdef MFEM_USE_MPI
559  TMOPNewtonSolver &tmopns_, ParGridFunction &px_,
560  bool move_bnd_, bool hradaptivity_,
561  int mesh_poly_deg_, int amr_metric_id_,
562  int hr_iter_, int h_per_r_iter_) :
563  mesh(&pmesh_), nlf(&pnlf_), tmopns(&tmopns_), x(&px_),
564  gridfuncarr(), fespacearr(),
565  move_bnd(move_bnd_), hradaptivity(hradaptivity_),
566  mesh_poly_deg(mesh_poly_deg_), amr_metric_id(amr_metric_id_),
567  pmesh(&pmesh_), pnlf(&pnlf_), pgridfuncarr(), pfespacearr(),
568  serial(false), hr_iter(hr_iter_), h_per_r_iter(h_per_r_iter_)
569 {
570  if (!hradaptivity) { return; }
572  amr_metric_id);
579 }
580 #endif
581 
583 {
584  Vector b(0);
585 #ifdef MFEM_USE_MPI
586  int myid = 0;
587 #endif
588  if (serial)
589  {
591  }
592  else
593  {
594 #ifdef MFEM_USE_MPI
595  myid = pnlf->ParFESpace()->GetMyRank();
597 #endif
598  }
599  if (!hradaptivity)
600  {
601  tmopns->Mult(b, x->GetTrueVector());
602  x->SetFromTrueVector();
603  return;
604  }
605 
606  bool radaptivity = true;
607 
608  tmop_dr->Reset();
609  tmop_r->Reset();
610 
611  if (serial)
612  {
613  for (int i_hr = 0; i_hr < hr_iter; i_hr++)
614  {
615  if (!radaptivity)
616  {
617  break;
618  }
619  mfem::out << i_hr << " r-adaptivity iteration.\n";
620 
622  tmopns->Mult(b, x->GetTrueVector());
623  x->SetFromTrueVector();
624 
625  mfem::out << "TMOP energy after r-adaptivity: " <<
627  ", Elements: " << mesh->GetNE() << std::endl;
628 
629  for (int i_h = 0; i_h < h_per_r_iter; i_h++)
630  {
631  // Derefinement step.
632  if (mesh->ncmesh)
633  {
634  tmop_dr->Apply(*mesh);
635  Update();
636  }
637  mfem::out << "TMOP energy after derefinement: " <<
639  ", Elements: " << mesh->GetNE() << std::endl;
640 
641  // Refinement step.
642  tmop_r->Apply(*mesh);
643  Update();
644  mfem::out << "TMOP energy after refinement: " <<
646  ", Elements: " << mesh->GetNE() << std::endl;
647 
648  if (!tmop_dr->Derefined() && tmop_r->Stop())
649  {
650  radaptivity = false;
651  mfem::out << "AMR stopping criterion satisfied. Stop.\n";
652  break;
653  }
654  } //n_h
655  } //n_hr
656  }
657  else
658  {
659 #ifdef MFEM_USE_MPI
660  int NEGlob;
661  double tmopenergy;
662  for (int i_hr = 0; i_hr < hr_iter; i_hr++)
663  {
664  if (!radaptivity)
665  {
666  break;
667  }
668  if (myid == 0) { mfem::out << i_hr << " r-adaptivity iteration.\n"; }
670  tmopns->Mult(b, x->GetTrueVector());
671  x->SetFromTrueVector();
672 
673  NEGlob = pmesh->GetGlobalNE();
674  tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
675  if (myid == 0)
676  {
677  mfem::out << "TMOP energy after r-adaptivity: " << tmopenergy <<
678  ", Elements: " << NEGlob << std::endl;
679  }
680 
681  for (int i_h = 0; i_h < h_per_r_iter; i_h++)
682  {
683  // Derefinement step.
684  if (pmesh->pncmesh)
685  {
687  ParUpdate();
688 
689  tmop_dr->Apply(*pmesh);
690  ParUpdate();
691  }
692  NEGlob = pmesh->GetGlobalNE();
693  tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
694  if (myid == 0)
695  {
696  mfem::out << "TMOP energy after derefinement: " << tmopenergy <<
697  ", Elements: " << NEGlob << std::endl;
698  }
699 
700  // Refinement step.
701  tmop_r->Apply(*pmesh);
702  ParUpdate();
703  NEGlob = pmesh->GetGlobalNE();
704  tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
705  if (myid == 0)
706  {
707  mfem::out << "TMOP energy after refinement: " << tmopenergy <<
708  ", Elements: " << NEGlob << std::endl;
709  }
710 
711  if (!tmop_dr->Derefined() && tmop_r->Stop())
712  {
713  radaptivity = false;
714  if (myid == 0)
715  {
716  mfem::out << "AMR stopping criterion satisfied. Stop.\n";
717  }
718  break;
719  }
720  } // n_r limit
721  } // n_hr
722 #endif
723  }
724 }
725 
726 #ifdef MFEM_USE_MPI
728 {
729  ParNCMesh *pncmesh = pmesh->pncmesh;
730  if (pncmesh)
731  {
732  const Table &dreftable = pncmesh->GetDerefinementTable();
733  Array<int> drefs, new_ranks;
734  for (int i = 0; i < dreftable.Size(); i++)
735  {
736  drefs.Append(i);
737  }
738  pncmesh->GetFineToCoarsePartitioning(drefs, new_ranks);
739  pmesh->Rebalance(new_ranks);
740  }
741 }
742 #endif
743 
745 {
746  // Update FESpace
747  for (int i = 0; i < fespacearr.Size(); i++)
748  {
749  fespacearr[i]->Update();
750  }
751  // Update nodal GF
752  for (int i = 0; i < gridfuncarr.Size(); i++)
753  {
754  gridfuncarr[i]->Update();
755  gridfuncarr[i]->SetTrueVector();
756  gridfuncarr[i]->SetFromTrueVector();
757  }
758 
759  // Update Discrete Indicator for all the TMOP_Integrators in NonLinearForm
761  TMOP_Integrator *ti = NULL;
762  TMOPComboIntegrator *co = NULL;
763  DiscreteAdaptTC *dtc = NULL;
764  for (int i = 0; i < integs.Size(); i++)
765  {
766  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
767  if (ti)
768  {
770  dtc = ti->GetDiscreteAdaptTC();
771  if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
772  }
773  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
774  if (co)
775  {
777  for (int j = 0; j < ati.Size(); j++)
778  {
779  ati[j]->UpdateAfterMeshTopologyChange();
780  dtc = ati[j]->GetDiscreteAdaptTC();
781  if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
782  }
783  }
784  }
785 
786  // Update the Nonlinear form and set Essential BC.
788 }
789 
790 #ifdef MFEM_USE_MPI
792 {
793  // Update FESpace
794  for (int i = 0; i < pfespacearr.Size(); i++)
795  {
796  pfespacearr[i]->Update();
797  }
798  // Update nodal GF
799  for (int i = 0; i < pgridfuncarr.Size(); i++)
800  {
801  pgridfuncarr[i]->Update();
802  pgridfuncarr[i]->SetTrueVector();
803  pgridfuncarr[i]->SetFromTrueVector();
804  }
805 
806  // Update Discrete Indicator
808  TMOP_Integrator *ti = NULL;
809  TMOPComboIntegrator *co = NULL;
810  DiscreteAdaptTC *dtc = NULL;
811  for (int i = 0; i < integs.Size(); i++)
812  {
813  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
814  if (ti)
815  {
817  dtc = ti->GetDiscreteAdaptTC();
818  if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
819  }
820  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
821  if (co)
822  {
824  for (int j = 0; j < ati.Size(); j++)
825  {
826  ati[j]->ParUpdateAfterMeshTopologyChange();
827  dtc = ati[j]->GetDiscreteAdaptTC();
828  if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
829  }
830  }
831  }
832 
833  // Update the Nonlinear form and set Essential BC.
835 }
836 #endif
837 
839 {
840  const FiniteElementSpace &fes = *mesh_->GetNodalFESpace();
841 
842  // Update Nonlinear form and Set Essential BC
843  nlf_->Update();
844  const int dim = fes.GetFE(0)->GetDim();
845  if (move_bnd == false)
846  {
847  Array<int> ess_bdr(mesh_->bdr_attributes.Max());
848  ess_bdr = 1;
849  nlf_->SetEssentialBC(ess_bdr);
850  }
851  else
852  {
853  const int nd = fes.GetBE(0)->GetDof();
854  int n = 0;
855  for (int i = 0; i < mesh_->GetNBE(); i++)
856  {
857  const int attr = mesh_->GetBdrElement(i)->GetAttribute();
858  MFEM_VERIFY(!(dim == 2 && attr == 3),
859  "Boundary attribute 3 must be used only for 3D meshes. "
860  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
861  "components, rest for free nodes), or use -fix-bnd.");
862  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
863  if (attr == 4) { n += nd * dim; }
864  }
865  Array<int> ess_vdofs(n), vdofs;
866  n = 0;
867  for (int i = 0; i < mesh_->GetNBE(); i++)
868  {
869  const int attr = mesh_->GetBdrElement(i)->GetAttribute();
870  fes.GetBdrElementVDofs(i, vdofs);
871  if (attr == 1) // Fix x components.
872  {
873  for (int j = 0; j < nd; j++)
874  { ess_vdofs[n++] = vdofs[j]; }
875  }
876  else if (attr == 2) // Fix y components.
877  {
878  for (int j = 0; j < nd; j++)
879  { ess_vdofs[n++] = vdofs[j+nd]; }
880  }
881  else if (attr == 3) // Fix z components.
882  {
883  for (int j = 0; j < nd; j++)
884  { ess_vdofs[n++] = vdofs[j+2*nd]; }
885  }
886  else if (attr == 4) // Fix all components.
887  {
888  for (int j = 0; j < vdofs.Size(); j++)
889  { ess_vdofs[n++] = vdofs[j]; }
890  }
891  }
892  nlf_->SetEssentialVDofs(ess_vdofs);
893  }
894 }
895 
896 }
Abstract class for all finite elements.
Definition: fe_base.hpp:233
void GetFineToCoarsePartitioning(const Array< int > &derefs, Array< int > &new_ranks) const
Definition: pncmesh.cpp:1361
virtual double GetDerefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Definition: tmop.cpp:3379
ThresholdRefiner * tmop_r
Definition: tmop_amr.hpp:215
bool Derefined() const
Check if the mesh was de-refined.
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:5630
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
const CoarseFineTransformations & GetDerefinementTransforms()
Definition: ncmesh.cpp:4623
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3402
void GetTMOPRefinementEnergy(int reftype, Vector &el_energy_vec)
Definition: tmop_amr.cpp:87
void SetFromTrueVector()
Shortcut for calling SetFromTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:150
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:762
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
static Mesh MakeCartesian3D(int nx, int ny, int nz, Element::Type type, double sx=1.0, double sy=1.0, double sz=1.0, bool sfc_ordering=true)
Definition: mesh.cpp:3783
const Array< TMOP_Integrator * > & GetTMOPIntegrators() const
Definition: tmop.hpp:2256
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3323
void ParUpdateAfterMeshTopologyChange()
Definition: tmop.cpp:1816
DiscreteAdaptTC * GetDiscreteAdaptTC() const
Definition: tmop.hpp:2204
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
long GetSequence() const
Definition: mesh.hpp:1982
int AddTriangle(int v1, int v2, int v3, int attr=1)
Definition: mesh.cpp:1680
Array< int > aniso_flags
Definition: tmop_amr.hpp:34
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void GetElementAverages(GridFunction &avgs) const
Definition: gridfunc.cpp:1798
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
Parallel non-linear operator on the true dofs.
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
Abstract parallel finite element space.
Definition: pfespace.hpp:28
TMOPHRSolver(Mesh &mesh_, NonlinearForm &nlf_, TMOPNewtonSolver &tmopns_, GridFunction &x_, bool move_bnd_, bool hradaptivity_, int mesh_poly_deg_, int amr_metric_id_, int hr_iter_=5, int h_per_r_iter_=1)
Definition: tmop_amr.cpp:535
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
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:2841
A parallel extension of the NCMesh class.
Definition: pncmesh.hpp:64
ParNCMesh * pncmesh
Definition: pmesh.hpp:388
Array< FiniteElementSpace * > fespacearr
Definition: tmop_amr.hpp:202
NonlinearForm * nlf
Definition: tmop_amr.hpp:28
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
bool Apply(Mesh &mesh)
Perform the mesh operation.
void Rebalance()
Definition: pmesh.cpp:4044
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1626
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
long long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:1115
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
Definition: mesh.cpp:1959
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
void SetEnergyScalingFactor(double scale)
Definition: tmop_amr.hpp:116
const Table & GetDerefinementTable()
Definition: ncmesh.cpp:1937
virtual void Reset()
Reset the associated estimator.
void FinalizeTetMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a tetrahedral Mesh.
Definition: mesh.cpp:2832
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
Array< NonlinearFormIntegrator * > * GetDNFI()
Access all integrators added with AddDomainIntegrator().
double b
Definition: lissajous.cpp:42
int AddTet(int v1, int v2, int v3, int v4, int attr=1)
Definition: mesh.cpp:1708
virtual void Reset()
Reset the associated estimator.
ThresholdDerefiner * tmop_dr
Definition: tmop_amr.hpp:217
Mesh refinement operator using an error threshold.
Array< ParGridFunction * > pgridfuncarr
Definition: tmop_amr.hpp:208
Array< IntegrationRule * > QuadIntRule
Definition: tmop_amr.hpp:31
GridFunction * spat_gf
Definition: tmop_amr.hpp:38
void Set2(const double x1, const double x2)
Definition: intrules.hpp:86
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5577
double GetGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
Array< IntegrationRule * > TetIntRule
Definition: tmop_amr.hpp:31
void GetTMOPDerefinementEnergy(Mesh &cmesh, TMOP_Integrator &tmopi, Vector &el_energy_vec)
Definition: tmop_amr.cpp:509
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
double GetParGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
virtual double GetRefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, const IntegrationRule &irule)
Computes the mean of the energies of the given element&#39;s children.
Definition: tmop.cpp:3293
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:76
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
void AddGridFunctionForUpdate(GridFunction *gf)
Definition: tmop_amr.hpp:252
virtual void Update()
Update the NonlinearForm to propagate updates of the associated FE space.
FiniteElementSpace * GetTSpecFESpace()
Get the FESpace associated with tspec.
Definition: tmop.hpp:1616
Array< GridFunction * > gridfuncarr
Definition: tmop_amr.hpp:201
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1419
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
ParFiniteElementSpace * ParFESpace() const
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 SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void RebalanceParNCMesh()
Definition: tmop_amr.cpp:727
const int mesh_poly_deg
Definition: tmop_amr.hpp:204
static Mesh MakeCartesian2D(int nx, int ny, Element::Type type, bool generate_edges=false, double sx=1.0, double sy=1.0, bool sfc_ordering=true)
Definition: mesh.cpp:3773
TMOPRefinerEstimator * tmop_r_est
Definition: tmop_amr.hpp:214
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:130
Array< IntegrationRule * > HexIntRule
Definition: tmop_amr.hpp:31
void UpdateAfterMeshTopologyChange()
Definition: tmop.cpp:3109
void MakeCoarseToFineTable(Table &coarse_to_fine, bool want_ghosts=false) const
Definition: ncmesh.cpp:4681
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:355
ParNonlinearForm * pnlf
Definition: tmop_amr.hpp:207
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i&#39;th boundary element. The returned indices are offsets int...
Definition: fespace.cpp:297
bool DerefineByError(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
Definition: mesh.cpp:9626
const int amr_metric_id
Definition: tmop_amr.hpp:204
int dim
Definition: ex24.cpp:53
IntegrationRule * SetIntRulesFromMesh(Mesh &meshsplit)
Definition: tmop_amr.cpp:317
void ParUpdateAfterMeshTopologyChange()
Definition: tmop.cpp:3122
TMOPDeRefinerEstimator * tmop_dr_est
Definition: tmop_amr.hpp:216
ParFiniteElementSpace * GetTSpecParFESpace()
Definition: tmop.hpp:1623
Array< IntegrationRule * > TriIntRule
Definition: tmop_amr.hpp:31
void SetEssentialBC(const Array< int > &bdr_attr_is_ess, Vector *rhs=NULL)
Specify essential boundary conditions.
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition: mesh.hpp:278
Array< ParFiniteElementSpace * > pfespacearr
Definition: tmop_amr.hpp:209
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:714
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:678
GridFunction * x
Definition: tmop_amr.hpp:200
void UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
Definition: tmop_amr.cpp:838
Vector data type.
Definition: vector.hpp:58
void SetTspecDataForDerefinement(FiniteElementSpace *fes)
Computes target specification data with respect to the coarse FE space.
Definition: tmop.cpp:2099
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
NonlinearForm * nlf
Definition: tmop_amr.hpp:198
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:70
void SetEssentialVDofs(const Array< int > &ess_vdofs_list)
Specify essential boundary conditions.
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9820
Class for parallel meshes.
Definition: pmesh.hpp:32
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: tmop_tools.hpp:261
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:3166
TMOPNewtonSolver * tmopns
Definition: tmop_amr.hpp:199
void UpdateAfterMeshTopologyChange()
Update all discrete fields based on tspec and update for AMR.
Definition: tmop.cpp:2005
void ResetDerefinementTspecData()
Definition: tmop.hpp:1713
bool GetDerefineEnergyForIntegrator(TMOP_Integrator &tmopi, Vector &fine_energy)
Definition: tmop_amr.cpp:357
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
De-refinement operator using an error threshold.
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1733