MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
check-tmop-metric.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 // Checks evaluation / 1st derivative / 2nd derivative for TMOP metrics. Serial.
13 // ./check-tmop-metric -mid 360
14 //
15 
16 #include "mfem.hpp"
17 #include <iostream>
18 
19 using namespace mfem;
20 using namespace std;
21 
22 int main(int argc, char *argv[])
23 {
24  int metric_id = 2;
25  int convergence_iter = 10;
26  bool verbose = false;
27 
28  // Choose metric.
29  OptionsParser args(argc, argv);
30  args.AddOption(&metric_id, "-mid", "--metric-id", "Metric id");
31  args.AddOption(&verbose, "-v", "-verbose", "-no-v", "--no-verbose",
32  "Enable extra screen output.");
33  args.AddOption(&convergence_iter, "-i", "--iterations",
34  "Number of iterations to check convergence of derivatives.");
35 
36  args.Parse();
37  if (!args.Good())
38  {
39  args.PrintUsage(cout);
40  return 1;
41  }
42  args.PrintOptions(cout);
43 
44  // Setup metric.
45  double tauval = -0.1;
46  TMOP_QualityMetric *metric = NULL;
47  switch (metric_id)
48  {
49  // T-metrics
50  case 1: metric = new TMOP_Metric_001; break;
51  case 2: metric = new TMOP_Metric_002; break;
52  case 7: metric = new TMOP_Metric_007; break;
53  case 9: metric = new TMOP_Metric_009; break;
54  case 14: metric = new TMOP_Metric_014; break;
55  case 22: metric = new TMOP_Metric_022(tauval); break;
56  case 50: metric = new TMOP_Metric_050; break;
57  case 55: metric = new TMOP_Metric_055; break;
58  case 56: metric = new TMOP_Metric_056; break;
59  case 58: metric = new TMOP_Metric_058; break;
60  case 77: metric = new TMOP_Metric_077; break;
61  case 80: metric = new TMOP_Metric_080(0.5); break;
62  case 85: metric = new TMOP_Metric_085; break;
63  case 98: metric = new TMOP_Metric_098; break;
64  // case 211: metric = new TMOP_Metric_211; break;
65  // case 252: metric = new TMOP_Metric_252(tauval); break;
66  case 301: metric = new TMOP_Metric_301; break;
67  case 302: metric = new TMOP_Metric_302; break;
68  case 303: metric = new TMOP_Metric_303; break;
69  case 304: metric = new TMOP_Metric_304; break;
70  // case 311: metric = new TMOP_Metric_311; break;
71  case 313: metric = new TMOP_Metric_313(tauval); break;
72  case 315: metric = new TMOP_Metric_315; break;
73  case 316: metric = new TMOP_Metric_316; break;
74  case 321: metric = new TMOP_Metric_321; break;
75  case 322: metric = new TMOP_Metric_322; break;
76  case 323: metric = new TMOP_Metric_323; break;
77  case 328: metric = new TMOP_Metric_328(0.5); break;
78  case 332: metric = new TMOP_Metric_332(0.5); break;
79  case 333: metric = new TMOP_Metric_333(0.5); break;
80  case 334: metric = new TMOP_Metric_334(0.5); break;
81  case 347: metric = new TMOP_Metric_347(0.5); break;
82  // case 352: metric = new TMOP_Metric_352(tauval); break;
83  case 360: metric = new TMOP_Metric_360; break;
84  // A-metrics
85  case 11: metric = new TMOP_AMetric_011; break;
86  case 36: metric = new TMOP_AMetric_036; break;
87  case 107: metric = new TMOP_AMetric_107a; break;
88  case 126: metric = new TMOP_AMetric_126(0.9); break;
89  default: cout << "Unknown metric_id: " << metric_id << endl; return 3;
90  }
91 
92  const int dim = (metric_id < 300) ? 2 : 3;
93  DenseMatrix T(dim);
94  Vector T_vec(T.GetData(), dim * dim);
95 
96  // Test evaluation.
97  int valid_cnt = 0, bad_cnt = 0;
98  for (int i = 0; i < 1000; i++)
99  {
100  T_vec.Randomize(i);
101  // Increase probability of det(T) > 0.
102  T(0, 0) += T_vec.Max();
103  if (T.Det() <= 0.0) { continue; }
104 
105  const double i_form = metric->EvalW(T),
106  m_form = metric->EvalWMatrixForm(T);
107  const double diff = fabs(i_form - m_form) / fabs(m_form);
108  if (diff > 1e-8)
109  {
110  bad_cnt++;
111  if (verbose)
112  {
113  cout << "Wrong metric computation: "
114  << i_form << " (invariant), " << m_form << " (matrix form) "
115  << diff << " (normalized difference) " << endl;
116  }
117  }
118  valid_cnt++;
119  }
120  cout << "--- EvalW: " << bad_cnt << " errors out of "
121  << valid_cnt << " comparisons with det(T) > 0.\n";
122 
123  Mesh *mesh;
124  if (dim == 2)
125  {
127  }
128  else
129  {
130  mesh = new Mesh(Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON));
131  }
132  H1_FECollection fec(2, dim);
133  FiniteElementSpace fespace(mesh, &fec, dim);
134  NonlinearForm a(&fespace);
135  mesh->SetNodalFESpace(&fespace);
136  GridFunction x(&fespace);
137  mesh->SetNodalGridFunction(&x);
138  x(0) = 0.25;
139 
141  tc.SetNodes(x);
142  auto integ = new TMOP_Integrator(metric, &tc, NULL);
143  a.AddDomainIntegrator(integ);
144 
146  const FiniteElement &fe = *fespace.GetFE(0);
147  Array<int> vdofs;
148  fespace.GetElementVDofs(0, vdofs);
149  Vector x_loc(x.Size());
150  x.GetSubVector(vdofs, x_loc);
151 
152  //
153  // Test 1st derivative (assuming EvalW is correct). Should be 2nd order.
154  //
155  Vector dF_0;
156  const double F_0 = integ->GetElementEnergy(fe, Tr, x_loc);
157  integ->AssembleElementVector(fe, Tr, x_loc, dF_0);
158  if (verbose) { cout << "***\ndF = \n"; dF_0.Print(); cout << "***\n"; }
159  double dx = 0.1;
160  double rate_dF_sum = 0.0, err_old;
161  for (int k = 0; k < convergence_iter; k++)
162  {
163  double err_k = 0.0;
164  for (int i = 0; i < x_loc.Size(); i++)
165  {
166  x_loc(i) += dx;
167  err_k = fmax(err_k, fabs(F_0 + dF_0(i) * dx -
168  integ->GetElementEnergy(fe, Tr, x_loc)));
169  x_loc(i) -= dx;
170  }
171  dx *= 0.5;
172 
173  if (verbose && k == 0)
174  {
175  std::cout << "dF error " << k << ": " << err_k << endl;
176  }
177  if (k > 0)
178  {
179  double r = log2(err_old / err_k);
180  rate_dF_sum += r;
181  if (verbose)
182  {
183  std::cout << "dF error " << k << ": " << err_k << " " << r << endl;
184  }
185  }
186  err_old = err_k;
187  }
188  std::cout << "--- EvalP: avg rate of convergence (should be 2): "
189  << rate_dF_sum / (convergence_iter - 1) << endl;
190 
191  //
192  // Test 2nd derivative (assuming EvalP is correct).
193  //
194  double min_avg_rate = 7.0;
195  DenseMatrix ddF_0;
196  integ->AssembleElementGrad(fe, Tr, x_loc, ddF_0);
197  if (verbose) { cout << "***\nddF = \n"; ddF_0.Print(); cout << "***\n"; }
198  for (int i = 0; i < x_loc.Size(); i++)
199  {
200  double rate_sum = 0.0;
201  double dx = 0.1;
202  for (int k = 0; k < convergence_iter; k++)
203  {
204  double err_k = 0.0;
205 
206  for (int j = 0; j < x_loc.Size(); j++)
207  {
208  x_loc(j) += dx;
209  Vector dF_dx;
210  integ->AssembleElementVector(fe, Tr, x_loc, dF_dx);
211  err_k = fmax(err_k, fabs(dF_0(i) + ddF_0(i, j) * dx - dF_dx(i)));
212  x_loc(j) -= dx;
213  }
214  dx *= 0.5;
215 
216  if (verbose && k == 0)
217  {
218  cout << "ddF error for dof " << i << ", " << k << ": "
219  << err_k << endl;
220  }
221  if (k > 0)
222  {
223  double r = log2(err_old / err_k);
224  rate_sum += r;
225  if (verbose)
226  {
227  cout << "ddF error for dof " << i << ", " << k << ": " << err_k
228  << " " << r << endl;
229  }
230 
231  }
232  err_old = err_k;
233  }
234  min_avg_rate = fmin(min_avg_rate, rate_sum / (convergence_iter - 1));
235  }
236  std::cout << "--- AssembleH: avg rate of convergence (should be 2): "
237  << min_avg_rate << endl;
238 
239  delete metric;
240  delete mesh;
241  return 0;
242 }
Abstract class for all finite elements.
Definition: fe_base.hpp:235
3D non-barrier Shape (S) metric.
Definition: tmop.hpp:939
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:744
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:872
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:625
double Det() const
Definition: densemat.cpp:449
virtual double EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
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:3723
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:493
void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
Adds new Domain Integrator.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:604
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
2D barrier shape (S) metric (polyconvex).
Definition: tmop.hpp:275
3D non-barrier metric without a type.
Definition: tmop.hpp:707
2D barrier size (V) metric (polyconvex).
Definition: tmop.hpp:475
2D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:1033
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:786
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:785
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:314
3D barrier Shape+Size (VS) metric, well-posed (invex).
Definition: tmop.hpp:765
2D Shifted barrier form of shape metric (mu_2).
Definition: tmop.hpp:363
2D barrier shape (S) metric (not polyconvex).
Definition: tmop.hpp:433
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:2183
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:646
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
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 Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
2D non-barrier size (V) metric (not polyconvex).
Definition: tmop.hpp:398
3D barrier metric without a type.
Definition: tmop.hpp:725
void SetNodes(const GridFunction &n)
Set the nodes to be used in the target-matrix construction.
Definition: tmop.hpp:1246
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:896
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:531
virtual void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:5285
virtual double EvalWMatrixForm(const DenseMatrix &Jpt) const
Evaluates the metric in matrix form (opposed to invariant form). Used for validating the invariant ev...
Definition: tmop.hpp:47
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:1015
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:807
2D barrier (not a shape) metric (polyconvex).
Definition: tmop.hpp:382
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:724
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
2D barrier size (V) metric (polyconvex).
Definition: tmop.hpp:415
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
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:3713
2D barrier Shape+Size (VS) metric (not polyconvex).
Definition: tmop.hpp:332
double a
Definition: lissajous.cpp:41
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
2D barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:997
3D Shape (S) metric, untangling version of 303.
Definition: tmop.hpp:686
int dim
Definition: ex24.cpp:53
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:2781
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Vector data type.
Definition: vector.hpp:60
3D barrier Shape+Size (VS) metric, well-posed (polyconvex).
Definition: tmop.hpp:851
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:5332
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1166
3D barrier Shape+Size (VS) metric (polyconvex).
Definition: tmop.hpp:828
2D barrier Shape+Orientation (OS) metric (polyconvex).
Definition: tmop.hpp:516
int main()
3D barrier Shape (S) metric, well-posed (polyconvex &amp; invex).
Definition: tmop.hpp:585
2D non-barrier metric without a type.
Definition: tmop.hpp:197
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
A TMOP integrator class based on any given TMOP_QualityMetric and TargetConstructor.
Definition: tmop.hpp:1565
2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
Definition: tmop.hpp:348