MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tmop_amr.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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  int myid = 0;
586  if (serial)
587  {
589  }
590  else
591  {
592 #ifdef MFEM_USE_MPI
593  myid = pnlf->ParFESpace()->GetMyRank();
595 #endif
596  }
597  if (!hradaptivity)
598  {
599  tmopns->Mult(b, x->GetTrueVector());
600  if (tmopns->GetConverged() == false)
601  {
602  if (myid == 0) { mfem::out << "Nonlinear solver: rtol not achieved.\n"; }
603  }
604  x->SetFromTrueVector();
605  return;
606  }
607 
608  bool radaptivity = true;
609 
610  tmop_dr->Reset();
611  tmop_r->Reset();
612 
613  if (serial)
614  {
615  for (int i_hr = 0; i_hr < hr_iter; i_hr++)
616  {
617  if (!radaptivity)
618  {
619  break;
620  }
621  mfem::out << i_hr << " r-adaptivity iteration.\n";
622 
624  tmopns->Mult(b, x->GetTrueVector());
625  x->SetFromTrueVector();
626 
627  mfem::out << "TMOP energy after r-adaptivity: " <<
629  ", Elements: " << mesh->GetNE() << std::endl;
630 
631  for (int i_h = 0; i_h < h_per_r_iter; i_h++)
632  {
633  // Derefinement step.
634  if (mesh->ncmesh)
635  {
636  tmop_dr->Apply(*mesh);
637  Update();
638  }
639  mfem::out << "TMOP energy after derefinement: " <<
641  ", Elements: " << mesh->GetNE() << std::endl;
642 
643  // Refinement step.
644  tmop_r->Apply(*mesh);
645  Update();
646  mfem::out << "TMOP energy after refinement: " <<
648  ", Elements: " << mesh->GetNE() << std::endl;
649 
650  if (!tmop_dr->Derefined() && tmop_r->Stop())
651  {
652  radaptivity = false;
653  mfem::out << "AMR stopping criterion satisfied. Stop.\n";
654  break;
655  }
656  } //n_h
657  } //n_hr
658  }
659  else
660  {
661 #ifdef MFEM_USE_MPI
662  int NEGlob;
663  double tmopenergy;
664  for (int i_hr = 0; i_hr < hr_iter; i_hr++)
665  {
666  if (!radaptivity)
667  {
668  break;
669  }
670  if (myid == 0) { mfem::out << i_hr << " r-adaptivity iteration.\n"; }
672  tmopns->Mult(b, x->GetTrueVector());
673  x->SetFromTrueVector();
674 
675  NEGlob = pmesh->GetGlobalNE();
676  tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
677  if (myid == 0)
678  {
679  mfem::out << "TMOP energy after r-adaptivity: " << tmopenergy <<
680  ", Elements: " << NEGlob << std::endl;
681  }
682 
683  for (int i_h = 0; i_h < h_per_r_iter; i_h++)
684  {
685  // Derefinement step.
686  if (pmesh->pncmesh)
687  {
689  ParUpdate();
690 
691  tmop_dr->Apply(*pmesh);
692  ParUpdate();
693  }
694  NEGlob = pmesh->GetGlobalNE();
695  tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
696  if (myid == 0)
697  {
698  mfem::out << "TMOP energy after derefinement: " << tmopenergy <<
699  ", Elements: " << NEGlob << std::endl;
700  }
701 
702  // Refinement step.
703  tmop_r->Apply(*pmesh);
704  ParUpdate();
705  NEGlob = pmesh->GetGlobalNE();
706  tmopenergy = pnlf->GetParGridFunctionEnergy(*x) / NEGlob;
707  if (myid == 0)
708  {
709  mfem::out << "TMOP energy after refinement: " << tmopenergy <<
710  ", Elements: " << NEGlob << std::endl;
711  }
712 
713  if (!tmop_dr->Derefined() && tmop_r->Stop())
714  {
715  radaptivity = false;
716  if (myid == 0)
717  {
718  mfem::out << "AMR stopping criterion satisfied. Stop.\n";
719  }
720  break;
721  }
722  } // n_r limit
723  } // n_hr
724 #endif
725  }
726 }
727 
728 #ifdef MFEM_USE_MPI
730 {
731  ParNCMesh *pncmesh = pmesh->pncmesh;
732  if (pncmesh)
733  {
734  const Table &dreftable = pncmesh->GetDerefinementTable();
735  Array<int> drefs, new_ranks;
736  for (int i = 0; i < dreftable.Size(); i++)
737  {
738  drefs.Append(i);
739  }
740  pncmesh->GetFineToCoarsePartitioning(drefs, new_ranks);
741  pmesh->Rebalance(new_ranks);
742  }
743 }
744 #endif
745 
747 {
748  // Update FESpace
749  for (int i = 0; i < fespacearr.Size(); i++)
750  {
751  fespacearr[i]->Update();
752  }
753  // Update nodal GF
754  for (int i = 0; i < gridfuncarr.Size(); i++)
755  {
756  gridfuncarr[i]->Update();
757  gridfuncarr[i]->SetTrueVector();
758  gridfuncarr[i]->SetFromTrueVector();
759  }
760 
761  // Update Discrete Indicator for all the TMOP_Integrators in NonLinearForm
763  TMOP_Integrator *ti = NULL;
764  TMOPComboIntegrator *co = NULL;
765  DiscreteAdaptTC *dtc = NULL;
766  for (int i = 0; i < integs.Size(); i++)
767  {
768  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
769  if (ti)
770  {
772  dtc = ti->GetDiscreteAdaptTC();
773  if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
774  }
775  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
776  if (co)
777  {
779  for (int j = 0; j < ati.Size(); j++)
780  {
781  ati[j]->UpdateAfterMeshTopologyChange();
782  dtc = ati[j]->GetDiscreteAdaptTC();
783  if (dtc) { dtc->UpdateAfterMeshTopologyChange(); }
784  }
785  }
786  }
787 
788  // Update the Nonlinear form and set Essential BC.
790 }
791 
792 #ifdef MFEM_USE_MPI
794 {
795  // Update FESpace
796  for (int i = 0; i < pfespacearr.Size(); i++)
797  {
798  pfespacearr[i]->Update();
799  }
800  // Update nodal GF
801  for (int i = 0; i < pgridfuncarr.Size(); i++)
802  {
803  pgridfuncarr[i]->Update();
804  pgridfuncarr[i]->SetTrueVector();
805  pgridfuncarr[i]->SetFromTrueVector();
806  }
807 
808  // Update Discrete Indicator
810  TMOP_Integrator *ti = NULL;
811  TMOPComboIntegrator *co = NULL;
812  DiscreteAdaptTC *dtc = NULL;
813  for (int i = 0; i < integs.Size(); i++)
814  {
815  ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
816  if (ti)
817  {
819  dtc = ti->GetDiscreteAdaptTC();
820  if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
821  }
822  co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
823  if (co)
824  {
826  for (int j = 0; j < ati.Size(); j++)
827  {
828  ati[j]->ParUpdateAfterMeshTopologyChange();
829  dtc = ati[j]->GetDiscreteAdaptTC();
830  if (dtc) { dtc->ParUpdateAfterMeshTopologyChange(); }
831  }
832  }
833  }
834 
835  // Update the Nonlinear form and set Essential BC.
837 }
838 #endif
839 
841 {
842  const FiniteElementSpace &fes = *mesh_->GetNodalFESpace();
843 
844  // Update Nonlinear form and Set Essential BC
845  nlf_->Update();
846  const int dim = fes.GetFE(0)->GetDim();
847  if (move_bnd == false)
848  {
849  Array<int> ess_bdr(mesh_->bdr_attributes.Max());
850  ess_bdr = 1;
851  nlf_->SetEssentialBC(ess_bdr);
852  }
853  else
854  {
855  const int nd = fes.GetBE(0)->GetDof();
856  int n = 0;
857  for (int i = 0; i < mesh_->GetNBE(); i++)
858  {
859  const int attr = mesh_->GetBdrElement(i)->GetAttribute();
860  MFEM_VERIFY(!(dim == 2 && attr == 3),
861  "Boundary attribute 3 must be used only for 3D meshes. "
862  "Adjust the attributes (1/2/3/4 for fixed x/y/z/all "
863  "components, rest for free nodes), or use -fix-bnd.");
864  if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
865  if (attr == 4) { n += nd * dim; }
866  }
867  Array<int> ess_vdofs(n), vdofs;
868  n = 0;
869  for (int i = 0; i < mesh_->GetNBE(); i++)
870  {
871  const int attr = mesh_->GetBdrElement(i)->GetAttribute();
872  fes.GetBdrElementVDofs(i, vdofs);
873  if (attr == 1) // Fix x components.
874  {
875  for (int j = 0; j < nd; j++)
876  { ess_vdofs[n++] = vdofs[j]; }
877  }
878  else if (attr == 2) // Fix y components.
879  {
880  for (int j = 0; j < nd; j++)
881  { ess_vdofs[n++] = vdofs[j+nd]; }
882  }
883  else if (attr == 3) // Fix z components.
884  {
885  for (int j = 0; j < nd; j++)
886  { ess_vdofs[n++] = vdofs[j+2*nd]; }
887  }
888  else if (attr == 4) // Fix all components.
889  {
890  for (int j = 0; j < vdofs.Size(); j++)
891  { ess_vdofs[n++] = vdofs[j]; }
892  }
893  }
894  nlf_->SetEssentialVDofs(ess_vdofs);
895  }
896 }
897 
898 }
Abstract class for all finite elements.
Definition: fe_base.hpp:235
virtual double GetDerefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Definition: tmop.cpp:2714
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
ThresholdRefiner * tmop_r
Definition: tmop_amr.hpp:215
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
ParFiniteElementSpace * ParFESpace() const
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void GetElementAverages(GridFunction &avgs) const
Definition: gridfunc.cpp:1800
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:703
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
const CoarseFineTransformations & GetDerefinementTransforms()
Definition: ncmesh.cpp:4399
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
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:734
long GetGlobalNE() const
Return the total (global) number of elements.
Definition: mesh.hpp:960
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:931
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:3711
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3284
void ParUpdateAfterMeshTopologyChange()
Definition: tmop.cpp:1343
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:534
int AddTriangle(int v1, int v2, int v3, int attr=1)
Definition: mesh.cpp:1648
Array< int > aniso_flags
Definition: tmop_amr.hpp:34
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
Parallel non-linear operator on the true dofs.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
bool Stop() const
Check if STOP action is requested, e.g. stopping criterion is satisfied.
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
bool GetConverged() const
Definition: solvers.hpp:247
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
A parallel extension of the NCMesh class.
Definition: pncmesh.hpp:64
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:297
ParNCMesh * pncmesh
Definition: pmesh.hpp:325
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
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:590
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
bool Apply(Mesh &mesh)
Perform the mesh operation.
long GetSequence() const
Definition: mesh.hpp:1639
void Rebalance()
Definition: pmesh.cpp:3903
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1601
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
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
void FinalizeTriMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a triangular Mesh.
Definition: mesh.cpp:1927
void SetEnergyScalingFactor(double scale)
Definition: tmop_amr.hpp:116
const Table & GetDerefinementTable()
Definition: ncmesh.cpp:1805
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:2765
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
virtual void GetFineToCoarsePartitioning(const Array< int > &derefs, Array< int > &new_ranks) const
Definition: pncmesh.cpp:1342
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
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:1676
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
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
GridFunction * spat_gf
Definition: tmop_amr.hpp:38
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
void Set2(const double x1, const double x2)
Definition: intrules.hpp:80
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5273
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
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:652
int Dimension() const
Definition: mesh.hpp:999
double GetParGridFunctionEnergy(const Vector &x) const
Compute the energy corresponding to the state x.
const Vector & GetTrueVector() const
Read only access to the (optional) internal true-dof Vector.
Definition: gridfunc.hpp:134
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:2628
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
bool Derefined() const
Check if the mesh was de-refined.
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
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:1195
Array< GridFunction * > gridfuncarr
Definition: tmop_amr.hpp:201
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1359
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int NumRows() const
Get the number of rows (size of output) of the Operator. Synonym with Height().
Definition: operator.hpp:69
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
void RebalanceParNCMesh()
Definition: tmop_amr.cpp:729
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:3701
TMOPRefinerEstimator * tmop_r_est
Definition: tmop_amr.hpp:214
DiscreteAdaptTC * GetDiscreteAdaptTC() const
Definition: tmop.hpp:1691
Array< IntegrationRule * > HexIntRule
Definition: tmop_amr.hpp:31
void MakeCoarseToFineTable(Table &coarse_to_fine, bool want_ghosts=false) const
Definition: ncmesh.cpp:4457
void UpdateAfterMeshTopologyChange()
Definition: tmop.cpp:2465
double GetGridFunctionEnergy(const Vector &x) const
Compute the enery corresponding to the state x.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
ParNonlinearForm * pnlf
Definition: tmop_amr.hpp:207
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:5307
bool DerefineByError(Array< double > &elem_error, double threshold, int nc_limit=0, int op=1)
Definition: mesh.cpp:9133
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:2478
TMOPDeRefinerEstimator * tmop_dr_est
Definition: tmop_amr.hpp:216
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
ParFiniteElementSpace * GetTSpecParFESpace()
Definition: tmop.hpp:1202
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:273
Array< ParFiniteElementSpace * > pfespacearr
Definition: tmop_amr.hpp:209
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:904
GridFunction * x
Definition: tmop_amr.hpp:200
void UpdateNonlinearFormAndBC(Mesh *mesh, NonlinearForm *nlf)
Definition: tmop_amr.cpp:840
virtual void Mult(const Vector &b, Vector &x) const
Solve the nonlinear system with right-hand side b.
Definition: tmop_tools.hpp:224
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
void SetTspecDataForDerefinement(FiniteElementSpace *fes)
Computes target specification data with respect to the coarse FE space.
Definition: tmop.cpp:1609
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
NonlinearForm * nlf
Definition: tmop_amr.hpp:198
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.
Class for parallel grid function.
Definition: pgridfunc.hpp:32
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:3102
void SetTotalErrorFraction(double fraction)
Set the total fraction used in the computation of the threshold. The default value is 1/2...
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 GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition: mesh.cpp:9326
Class for parallel meshes.
Definition: pmesh.hpp:32
const Array< TMOP_Integrator * > & GetTMOPIntegrators() const
Definition: tmop.hpp:1737
TMOPNewtonSolver * tmopns
Definition: tmop_amr.hpp:199
void UpdateAfterMeshTopologyChange()
Update all discrete fields based on tspec and update for AMR.
Definition: tmop.cpp:1516
void ResetDerefinementTspecData()
Definition: tmop.hpp:1283
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1037
bool GetDerefineEnergyForIntegrator(TMOP_Integrator &tmopi, Vector &fine_energy)
Definition: tmop_amr.cpp:357
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:284
De-refinement operator using an error threshold.
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1303