MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tmop.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.hpp"
13 #include "linearform.hpp"
14 #include "pgridfunc.hpp"
15 #include "tmop_tools.hpp"
16 #include "../general/forall.hpp"
17 
18 namespace mfem
19 {
20 
21 // Target-matrix optimization paradigm (TMOP) mesh quality metrics.
22 
23 double TMOP_Combo_QualityMetric::EvalWMatrixForm(const DenseMatrix &Jpt) const
24 {
25  double metric = 0.;
26  for (int i = 0; i < tmop_q_arr.Size(); i++)
27  {
28  metric += wt_arr[i]*tmop_q_arr[i]->EvalWMatrixForm(Jpt);
29  }
30  return metric;
31 }
32 
33 double TMOP_Combo_QualityMetric::EvalW(const DenseMatrix &Jpt) const
34 {
35  double metric = 0.;
36  for (int i = 0; i < tmop_q_arr.Size(); i++)
37  {
38  metric += wt_arr[i]*tmop_q_arr[i]->EvalW(Jpt);
39  }
40  return metric;
41 }
42 
43 void TMOP_Combo_QualityMetric::EvalP(const DenseMatrix &Jpt,
44  DenseMatrix &P) const
45 {
46  DenseMatrix Pt(P.Size());
47  P = 0.0;
48  for (int i = 0; i < tmop_q_arr.Size(); i++)
49  {
50  tmop_q_arr[i]->EvalP(Jpt, Pt);
51  Pt *= wt_arr[i];
52  P += Pt;
53  }
54 }
55 
56 void TMOP_Combo_QualityMetric::AssembleH(const DenseMatrix &Jpt,
57  const DenseMatrix &DS,
58  const double weight,
59  DenseMatrix &A) const
60 {
61  DenseMatrix At(A.Size());
62  for (int i = 0; i < tmop_q_arr.Size(); i++)
63  {
64  At = 0.0;
65  tmop_q_arr[i]->AssembleH(Jpt, DS, weight, At);
66  At *= wt_arr[i];
67  A += At;
68  }
69 }
70 
71 double TMOP_WorstCaseUntangleOptimizer_Metric::EvalW(const DenseMatrix &Jpt)
72 const
73 {
74  double metric_tilde = EvalWBarrier(Jpt);
75  double metric = metric_tilde;
76  if (wctype == WorstCaseType::PMean)
77  {
78  metric = std::pow(metric_tilde, exponent);
79  }
80  else if (wctype == WorstCaseType::Beta)
81  {
82  double beta = max_muT+muT_ep;
83  metric = metric_tilde/(beta-metric_tilde);
84  }
85  return metric;
86 }
87 
88 double TMOP_WorstCaseUntangleOptimizer_Metric::EvalWBarrier(
89  const DenseMatrix &Jpt) const
90 {
91  double denominator = 1.0;
92  if (btype == BarrierType::Shifted)
93  {
94  denominator = 2.0*(Jpt.Det()-std::min(alpha*min_detT-detT_ep, 0.0));
95  }
96  else if (btype == BarrierType::Pseudo)
97  {
98  double detT = Jpt.Det();
99  denominator = detT + std::sqrt(detT*detT + detT_ep*detT_ep);
100  }
101  return tmop_metric.EvalW(Jpt)/denominator;
102 }
103 
104 double TMOP_Metric_001::EvalW(const DenseMatrix &Jpt) const
105 {
106  ie.SetJacobian(Jpt.GetData());
107  return ie.Get_I1();
108 }
109 
110 void TMOP_Metric_001::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
111 {
112  ie.SetJacobian(Jpt.GetData());
113  P = ie.Get_dI1();
114 }
115 
116 void TMOP_Metric_001::AssembleH(const DenseMatrix &Jpt,
117  const DenseMatrix &DS,
118  const double weight,
119  DenseMatrix &A) const
120 {
121  ie.SetJacobian(Jpt.GetData());
122  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
123  ie.Assemble_ddI1(weight, A.GetData());
124 }
125 
126 double TMOP_Metric_skew2D::EvalW(const DenseMatrix &Jpt) const
127 {
128  MFEM_VERIFY(Jtr != NULL,
129  "Requires a target Jacobian, use SetTargetJacobian().");
130 
131  DenseMatrix Jpr(2, 2);
132  Mult(Jpt, *Jtr, Jpr);
133 
134  Vector col1, col2;
135  Jpr.GetColumn(0, col1);
136  Jpr.GetColumn(1, col2);
137  double norm_prod = col1.Norml2() * col2.Norml2();
138  const double cos_Jpr = (col1 * col2) / norm_prod,
139  sin_Jpr = fabs(Jpr.Det()) / norm_prod;
140 
141  Jtr->GetColumn(0, col1);
142  Jtr->GetColumn(1, col2);
143  norm_prod = col1.Norml2() * col2.Norml2();
144  const double cos_Jtr = (col1 * col2) / norm_prod,
145  sin_Jtr = fabs(Jtr->Det()) / norm_prod;
146 
147  return 0.5 * (1.0 - cos_Jpr * cos_Jtr - sin_Jpr * sin_Jtr);
148 }
149 
150 double TMOP_Metric_skew3D::EvalW(const DenseMatrix &Jpt) const
151 {
152  MFEM_VERIFY(Jtr != NULL,
153  "Requires a target Jacobian, use SetTargetJacobian().");
154 
155  DenseMatrix Jpr(3, 3);
156  Mult(Jpt, *Jtr, Jpr);
157 
158  Vector col1, col2, col3;
159  Jpr.GetColumn(0, col1);
160  Jpr.GetColumn(1, col2);
161  Jpr.GetColumn(2, col3);
162  double norm_c1 = col1.Norml2(),
163  norm_c2 = col2.Norml2(),
164  norm_c3 = col3.Norml2();
165  double cos_Jpr_12 = (col1 * col2) / (norm_c1 * norm_c2),
166  cos_Jpr_13 = (col1 * col3) / (norm_c1 * norm_c3),
167  cos_Jpr_23 = (col2 * col3) / (norm_c2 * norm_c3);
168  double sin_Jpr_12 = std::sqrt(1.0 - cos_Jpr_12 * cos_Jpr_12),
169  sin_Jpr_13 = std::sqrt(1.0 - cos_Jpr_13 * cos_Jpr_13),
170  sin_Jpr_23 = std::sqrt(1.0 - cos_Jpr_23 * cos_Jpr_23);
171 
172  Jtr->GetColumn(0, col1);
173  Jtr->GetColumn(1, col2);
174  Jtr->GetColumn(2, col3);
175  norm_c1 = col1.Norml2();
176  norm_c2 = col2.Norml2(),
177  norm_c3 = col3.Norml2();
178  double cos_Jtr_12 = (col1 * col2) / (norm_c1 * norm_c2),
179  cos_Jtr_13 = (col1 * col3) / (norm_c1 * norm_c3),
180  cos_Jtr_23 = (col2 * col3) / (norm_c2 * norm_c3);
181  double sin_Jtr_12 = std::sqrt(1.0 - cos_Jtr_12 * cos_Jtr_12),
182  sin_Jtr_13 = std::sqrt(1.0 - cos_Jtr_13 * cos_Jtr_13),
183  sin_Jtr_23 = std::sqrt(1.0 - cos_Jtr_23 * cos_Jtr_23);
184 
185  return (3.0 - cos_Jpr_12 * cos_Jtr_12 - sin_Jpr_12 * sin_Jtr_12
186  - cos_Jpr_13 * cos_Jtr_13 - sin_Jpr_13 * sin_Jtr_13
187  - cos_Jpr_23 * cos_Jtr_23 - sin_Jpr_23 * sin_Jtr_23) / 6.0;
188 }
189 
190 double TMOP_Metric_aspratio2D::EvalW(const DenseMatrix &Jpt) const
191 {
192  MFEM_VERIFY(Jtr != NULL,
193  "Requires a target Jacobian, use SetTargetJacobian().");
194 
195  DenseMatrix Jpr(2, 2);
196  Mult(Jpt, *Jtr, Jpr);
197 
198  Vector col1, col2;
199  Jpr.GetColumn(0, col1);
200  Jpr.GetColumn(1, col2);
201  const double ratio_Jpr = col2.Norml2() / col1.Norml2();
202 
203  Jtr->GetColumn(0, col1);
204  Jtr->GetColumn(1, col2);
205  const double ratio_Jtr = col2.Norml2() / col1.Norml2();
206 
207  return 0.5 * (ratio_Jpr / ratio_Jtr + ratio_Jtr / ratio_Jpr) - 1.0;
208 }
209 
210 double TMOP_Metric_aspratio3D::EvalW(const DenseMatrix &Jpt) const
211 {
212  MFEM_VERIFY(Jtr != NULL,
213  "Requires a target Jacobian, use SetTargetJacobian().");
214 
215  DenseMatrix Jpr(3, 3);
216  Mult(Jpt, *Jtr, Jpr);
217 
218  Vector col1, col2, col3;
219  Jpr.GetColumn(0, col1);
220  Jpr.GetColumn(1, col2);
221  Jpr.GetColumn(2, col3);
222  double norm_c1 = col1.Norml2(),
223  norm_c2 = col2.Norml2(),
224  norm_c3 = col3.Norml2();
225  double ratio_Jpr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
226  ratio_Jpr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
227  ratio_Jpr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
228 
229  Jtr->GetColumn(0, col1);
230  Jtr->GetColumn(1, col2);
231  Jtr->GetColumn(2, col3);
232  norm_c1 = col1.Norml2();
233  norm_c2 = col2.Norml2();
234  norm_c3 = col3.Norml2();
235  double ratio_Jtr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
236  ratio_Jtr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
237  ratio_Jtr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
238 
239  return ( 0.5 * (ratio_Jpr_1 / ratio_Jtr_1 + ratio_Jtr_1 / ratio_Jpr_1) +
240  0.5 * (ratio_Jpr_2 / ratio_Jtr_2 + ratio_Jtr_2 / ratio_Jpr_2) +
241  0.5 * (ratio_Jpr_3 / ratio_Jtr_3 + ratio_Jtr_3 / ratio_Jpr_3) - 3.0
242  ) / 3.0;
243 }
244 
245 double TMOP_Metric_002::EvalWMatrixForm(const DenseMatrix &Jpt) const
246 {
247  return 0.5 * Jpt.FNorm2() / Jpt.Det() - 1.0;
248 }
249 
250 double TMOP_Metric_002::EvalW(const DenseMatrix &Jpt) const
251 {
252  ie.SetJacobian(Jpt.GetData());
253  return 0.5 * ie.Get_I1b() - 1.0;
254 }
255 
256 void TMOP_Metric_002::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
257 {
258  ie.SetJacobian(Jpt.GetData());
259  P.Set(0.5, ie.Get_dI1b());
260 }
261 
262 void TMOP_Metric_002::AssembleH(const DenseMatrix &Jpt,
263  const DenseMatrix &DS,
264  const double weight,
265  DenseMatrix &A) const
266 {
267  ie.SetJacobian(Jpt.GetData());
268  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
269  ie.Assemble_ddI1b(0.5*weight, A.GetData());
270 }
271 
272 double TMOP_Metric_004::EvalW(const DenseMatrix &Jpt) const
273 {
274  ie.SetJacobian(Jpt.GetData());
275  return ie.Get_I1() - 2.0*ie.Get_I2b();
276 }
277 
278 void TMOP_Metric_004::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
279 {
280  ie.SetJacobian(Jpt.GetData());
281  Add(1.0, ie.Get_dI1(), -2.0, ie.Get_dI2b(), P);
282 }
283 
284 void TMOP_Metric_004::AssembleH(const DenseMatrix &Jpt,
285  const DenseMatrix &DS,
286  const double weight,
287  DenseMatrix &A) const
288 {
289  ie.SetJacobian(Jpt.GetData());
290  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
291 
292  ie.Assemble_ddI1(weight, A.GetData());
293  ie.Assemble_ddI2b(-2.0*weight, A.GetData());
294 }
295 
296 double TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
297 {
298  // mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
299  ie.SetJacobian(Jpt.GetData());
300  return ie.Get_I1()*(1. + 1./ie.Get_I2()) - 4.0;
301 }
302 
303 void TMOP_Metric_007::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
304 {
305  // P = d(I1*(1 + 1/I2)) = (1 + 1/I2) dI1 - I1/I2^2 dI2
306  ie.SetJacobian(Jpt.GetData());
307  const double I2 = ie.Get_I2();
308  Add(1. + 1./I2, ie.Get_dI1(), -ie.Get_I1()/(I2*I2), ie.Get_dI2(), P);
309 }
310 
311 void TMOP_Metric_007::AssembleH(const DenseMatrix &Jpt,
312  const DenseMatrix &DS,
313  const double weight,
314  DenseMatrix &A) const
315 {
316  // P = d(I1*(1 + 1/I2))
317  // = (1 + 1/I2) dI1 - I1/I2^2 dI2
318  //
319  // dP = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 -
320  // (dI2 x d(I1/I2^2)) - I1/I2^2 ddI2
321  // = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 +
322  // (-1/I2^2) (dI2 x [dI1 - 2 I1/I2 dI2]) - I1/I2^2 ddI2
323  // = (-1/I2^2) (dI1 x dI2 + dI2 x dI1) + (1 + 1/I2) ddI1 +
324  // (2 I1/I2^3) (dI2 x dI2) - I1/I2^2 ddI2
325  ie.SetJacobian(Jpt.GetData());
326  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
327  const double c1 = 1./ie.Get_I2();
328  const double c2 = weight*c1*c1;
329  const double c3 = ie.Get_I1()*c2;
330  ie.Assemble_ddI1(weight*(1. + c1), A.GetData());
331  ie.Assemble_ddI2(-c3, A.GetData());
332  ie.Assemble_TProd(-c2, ie.Get_dI1(), ie.Get_dI2(), A.GetData());
333  ie.Assemble_TProd(2*c1*c3, ie.Get_dI2(), A.GetData());
334 }
335 
336 double TMOP_Metric_009::EvalW(const DenseMatrix &Jpt) const
337 {
338  // mu_9 = det(J)*|J-J^{-t}|^2 = I1b * (I2b^2 + 1) - 4 * I2b
339  // = (I1 - 4)*I2b + I1b
340  ie.SetJacobian(Jpt.GetData());
341  return (ie.Get_I1() - 4.0)*ie.Get_I2b() + ie.Get_I1b();
342 }
343 
344 void TMOP_Metric_009::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
345 {
346  // mu_9 = (I1 - 4)*I2b + I1b
347  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
348  ie.SetJacobian(Jpt.GetData());
349  Add(ie.Get_I1() - 4.0, ie.Get_dI2b(), ie.Get_I2b(), ie.Get_dI1(), P);
350  P += ie.Get_dI1b();
351 }
352 
353 void TMOP_Metric_009::AssembleH(const DenseMatrix &Jpt,
354  const DenseMatrix &DS,
355  const double weight,
356  DenseMatrix &A) const
357 {
358  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
359  // dP = dI2b x dI1 + (I1-4)*ddI2b + dI1 x dI2b + I2b*ddI1 + ddI1b
360  // = (dI1 x dI2b + dI2b x dI1) + (I1-4)*ddI2b + I2b*ddI1 + ddI1b
361  ie.SetJacobian(Jpt.GetData());
362  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
363  ie.Assemble_TProd(weight, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
364  ie.Assemble_ddI2b(weight*(ie.Get_I1()-4.0), A.GetData());
365  ie.Assemble_ddI1(weight*ie.Get_I2b(), A.GetData());
366  ie.Assemble_ddI1b(weight, A.GetData());
367 }
368 
369 // mu_14 = |T-I|^2
370 double TMOP_Metric_014::EvalW(const DenseMatrix &Jpt) const
371 {
372  MFEM_VERIFY(Jtr != NULL,
373  "Requires a target Jacobian, use SetTargetJacobian().");
374 
375  DenseMatrix Id(2,2);
376 
377  Id(0,0) = 1; Id(0,1) = 0;
378  Id(1,0) = 0; Id(1,1) = 1;
379 
380  DenseMatrix Mat(2,2);
381  Mat = Jpt;
382  Mat.Add(-1,Id);
383  return Mat.FNorm2();
384 }
385 
386 double TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const
387 {
388  // mu_22 = (0.5*|J|^2 - det(J)) / (det(J) - tau0)
389  // = (0.5*I1 - I2b) / (I2b - tau0)
390  ie.SetJacobian(Jpt.GetData());
391  const double I2b = ie.Get_I2b();
392 
393  double d = I2b - min_detT;
394  if (d < 0.0 && min_detT == 0.0)
395  {
396  // The mesh has been untangled, but it's still possible to get negative
397  // detJ in FD calculations, as they move the nodes around with some small
398  // increments and can produce negative determinants. Thus we put a small
399  // value in the denominator. Note that here I2b < 0.
400  d = - I2b * 0.1;
401  }
402 
403  return (0.5*ie.Get_I1() - I2b) / d;
404 }
405 
406 void TMOP_Metric_022::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
407 {
408  // mu_22 = (0.5*I1 - I2b) / (I2b - tau0)
409  // P = 1/(I2b - tau0)*(0.5*dI1 - dI2b) - (0.5*I1 - I2b)/(I2b - tau0)^2*dI2b
410  // = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
411  ie.SetJacobian(Jpt.GetData());
412  const double c1 = 1.0/(ie.Get_I2b() - min_detT);
413  Add(c1/2, ie.Get_dI1(), (min_detT - ie.Get_I1()/2)*c1*c1, ie.Get_dI2b(), P);
414 }
415 
416 void TMOP_Metric_022::AssembleH(const DenseMatrix &Jpt,
417  const DenseMatrix &DS,
418  const double weight,
419  DenseMatrix &A) const
420 {
421  // P = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
422  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b) + 0.5/(I2b - tau0)*ddI1
423  // + (dI2b x dz) + z*ddI2b
424  //
425  // z = (tau0 - 0.5*I1)/(I2b - tau0)^2
426  // dz = -0.5/(I2b - tau0)^2*dI1 - 2*(tau0 - 0.5*I1)/(I2b - tau0)^3*dI2b
427  //
428  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b + dI2b x dI1)
429  // -2*z/(I2b - tau0)*(dI2b x dI2b)
430  // +0.5/(I2b - tau0)*ddI1 + z*ddI2b
431  ie.SetJacobian(Jpt.GetData());
432  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
433  const double c1 = 1.0/(ie.Get_I2b() - min_detT);
434  const double c2 = weight*c1/2;
435  const double c3 = c1*c2;
436  const double c4 = (2*min_detT - ie.Get_I1())*c3; // weight*z
437  ie.Assemble_TProd(-c3, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
438  ie.Assemble_TProd(-2*c1*c4, ie.Get_dI2b(), A.GetData());
439  ie.Assemble_ddI1(c2, A.GetData());
440  ie.Assemble_ddI2b(c4, A.GetData());
441 }
442 
443 double TMOP_Metric_050::EvalW(const DenseMatrix &Jpt) const
444 {
445  // mu_50 = 0.5*|J^t J|^2/det(J)^2 - 1
446  // = 0.5*(l1^4 + l2^4)/(l1*l2)^2 - 1
447  // = 0.5*((l1/l2)^2 + (l2/l1)^2) - 1 = 0.5*(l1/l2 - l2/l1)^2
448  // = 0.5*(l1/l2 + l2/l1)^2 - 2 = 0.5*I1b^2 - 2
449  ie.SetJacobian(Jpt.GetData());
450  const double I1b = ie.Get_I1b();
451  return 0.5*I1b*I1b - 2.0;
452 }
453 
454 void TMOP_Metric_050::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
455 {
456  // mu_50 = 0.5*I1b^2 - 2
457  // P = I1b*dI1b
458  ie.SetJacobian(Jpt.GetData());
459  P.Set(ie.Get_I1b(), ie.Get_dI1b());
460 }
461 
462 void TMOP_Metric_050::AssembleH(const DenseMatrix &Jpt,
463  const DenseMatrix &DS,
464  const double weight,
465  DenseMatrix &A) const
466 {
467  // P = I1b*dI1b
468  // dP = dI1b x dI1b + I1b*ddI1b
469  ie.SetJacobian(Jpt.GetData());
470  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
471  ie.Assemble_TProd(weight, ie.Get_dI1b(), A.GetData());
472  ie.Assemble_ddI1b(weight*ie.Get_I1b(), A.GetData());
473 }
474 
475 double TMOP_Metric_055::EvalW(const DenseMatrix &Jpt) const
476 {
477  // mu_55 = (det(J) - 1)^2 = (I2b - 1)^2
478  ie.SetJacobian(Jpt.GetData());
479  const double c1 = ie.Get_I2b() - 1.0;
480  return c1*c1;
481 }
482 
483 void TMOP_Metric_055::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
484 {
485  // mu_55 = (I2b - 1)^2
486  // P = 2*(I2b - 1)*dI2b
487  ie.SetJacobian(Jpt.GetData());
488  P.Set(2*(ie.Get_I2b() - 1.0), ie.Get_dI2b());
489 }
490 
491 void TMOP_Metric_055::AssembleH(const DenseMatrix &Jpt,
492  const DenseMatrix &DS,
493  const double weight,
494  DenseMatrix &A) const
495 {
496  // P = 2*(I2b - 1)*dI2b
497  // dP = 2*(dI2b x dI2b) + 2*(I2b - 1)*ddI2b
498  ie.SetJacobian(Jpt.GetData());
499  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
500  ie.Assemble_TProd(2*weight, ie.Get_dI2b(), A.GetData());
501  ie.Assemble_ddI2b(2*weight*(ie.Get_I2b() - 1.0), A.GetData());
502 }
503 
504 double TMOP_Metric_056::EvalW(const DenseMatrix &Jpt) const
505 {
506  // mu_56 = 0.5*(I2b + 1/I2b) - 1
507  ie.SetJacobian(Jpt.GetData());
508  const double I2b = ie.Get_I2b();
509  return 0.5*(I2b + 1.0/I2b) - 1.0;
510 }
511 
512 void TMOP_Metric_056::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
513 {
514  // mu_56 = 0.5*(I2b + 1/I2b) - 1
515  // P = 0.5*(1 - 1/I2b^2)*dI2b
516  ie.SetJacobian(Jpt.GetData());
517  P.Set(0.5 - 0.5/ie.Get_I2(), ie.Get_dI2b());
518 }
519 
520 void TMOP_Metric_056::AssembleH(const DenseMatrix &Jpt,
521  const DenseMatrix &DS,
522  const double weight,
523  DenseMatrix &A) const
524 {
525  // P = 0.5*(1 - 1/I2b^2)*dI2b
526  // dP = (1/I2b^3)*(dI2b x dI2b) + (0.5 - 0.5/I2)*ddI2b
527  ie.SetJacobian(Jpt.GetData());
528  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
529  ie.Assemble_TProd(weight/(ie.Get_I2()*ie.Get_I2b()),
530  ie.Get_dI2b(), A.GetData());
531  ie.Assemble_ddI2b(weight*(0.5 - 0.5/ie.Get_I2()), A.GetData());
532 }
533 
534 double TMOP_Metric_058::EvalWMatrixForm(const DenseMatrix &Jpt) const
535 {
536  // mu_58 = |J^t J|^2 / det(J)^2 - 2|J|^2 / det(J) + 2
537  DenseMatrix JtJ(2);
538  MultAAt(Jpt, JtJ);
539  JtJ.Transpose();
540  double det = Jpt.Det();
541 
542  return JtJ.FNorm2()/(det*det) - 2*Jpt.FNorm2()/det + 2.0;
543 }
544 
545 double TMOP_Metric_058::EvalW(const DenseMatrix &Jpt) const
546 {
547  // mu_58 = I1b*(I1b - 2)
548  ie.SetJacobian(Jpt.GetData());
549  const double I1b = ie.Get_I1b();
550  return I1b*(I1b - 2.0);
551 }
552 
553 void TMOP_Metric_058::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
554 {
555  // mu_58 = I1b*(I1b - 2)
556  // P = (2*I1b - 2)*dI1b
557  ie.SetJacobian(Jpt.GetData());
558  P.Set(2*ie.Get_I1b() - 2.0, ie.Get_dI1b());
559 }
560 
561 void TMOP_Metric_058::AssembleH(const DenseMatrix &Jpt,
562  const DenseMatrix &DS,
563  const double weight,
564  DenseMatrix &A) const
565 {
566  // P = (2*I1b - 2)*dI1b
567  // dP = 2*(dI1b x dI1b) + (2*I1b - 2)*ddI1b
568  ie.SetJacobian(Jpt.GetData());
569  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
570  ie.Assemble_TProd(2*weight, ie.Get_dI1b(), A.GetData());
571  ie.Assemble_ddI1b(weight*(2*ie.Get_I1b() - 2.0), A.GetData());
572 }
573 
574 double TMOP_Metric_077::EvalW(const DenseMatrix &Jpt) const
575 {
576  ie.SetJacobian(Jpt.GetData());
577  const double I2b = ie.Get_I2b();
578  return 0.5*(I2b*I2b + 1./(I2b*I2b) - 2.);
579 }
580 
581 void TMOP_Metric_077::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
582 {
583  // Using I2b^2 = I2.
584  // dmu77_dJ = 1/2 (1 - 1/I2^2) dI2_dJ.
585  ie.SetJacobian(Jpt.GetData());
586  const double I2 = ie.Get_I2();
587  P.Set(0.5 * (1.0 - 1.0 / (I2 * I2)), ie.Get_dI2());
588 }
589 
590 void TMOP_Metric_077::AssembleH(const DenseMatrix &Jpt,
591  const DenseMatrix &DS,
592  const double weight,
593  DenseMatrix &A) const
594 {
595  ie.SetJacobian(Jpt.GetData());
596  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
597  const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
598  ie.Assemble_ddI2(weight*0.5*(1.0 - I2inv_sq), A.GetData());
599  ie.Assemble_TProd(weight * I2inv_sq / I2, ie.Get_dI2(), A.GetData());
600 }
601 
602 // mu_85 = |T-T'|^2, where T'= |T|*I/sqrt(2)
603 double TMOP_Metric_085::EvalW(const DenseMatrix &Jpt) const
604 {
605  MFEM_VERIFY(Jtr != NULL,
606  "Requires a target Jacobian, use SetTargetJacobian().");
607 
608  DenseMatrix Id(2,2);
609  DenseMatrix Mat(2,2);
610  Mat = Jpt;
611 
612  Id(0,0) = 1; Id(0,1) = 0;
613  Id(1,0) = 0; Id(1,1) = 1;
614  Id *= Mat.FNorm()/pow(2,0.5);
615 
616  Mat.Add(-1.,Id);
617  return Mat.FNorm2();
618 }
619 
620 // mu_98 = 1/(tau)|T-I|^2
621 double TMOP_Metric_098::EvalW(const DenseMatrix &Jpt) const
622 {
623  MFEM_VERIFY(Jtr != NULL,
624  "Requires a target Jacobian, use SetTargetJacobian().");
625 
626  DenseMatrix Id(2,2);
627 
628  Id(0,0) = 1; Id(0,1) = 0;
629  Id(1,0) = 0; Id(1,1) = 1;
630 
631  DenseMatrix Mat(2,2);
632  Mat = Jpt;
633  Mat.Add(-1,Id);
634  return Mat.FNorm2()/Jtr->Det();
635 }
636 
637 double TMOP_Metric_211::EvalW(const DenseMatrix &Jpt) const
638 {
639  // mu_211 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
640  // = (I2b - 1)^2 - I2b + sqrt(I2b^2 + eps)
641  ie.SetJacobian(Jpt.GetData());
642  const double I2b = ie.Get_I2b();
643  return (I2b - 1.0)*(I2b - 1.0) - I2b + std::sqrt(I2b*I2b + eps);
644 }
645 
646 void TMOP_Metric_211::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
647 {
648  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
649 }
650 
651 void TMOP_Metric_211::AssembleH(const DenseMatrix &Jpt,
652  const DenseMatrix &DS,
653  const double weight,
654  DenseMatrix &A) const
655 {
656  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
657 }
658 
659 double TMOP_Metric_252::EvalW(const DenseMatrix &Jpt) const
660 {
661  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0).
662  ie.SetJacobian(Jpt.GetData());
663  const double I2b = ie.Get_I2b();
664  return 0.5*(I2b - 1.0)*(I2b - 1.0)/(I2b - tau0);
665 }
666 
667 void TMOP_Metric_252::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
668 {
669  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
670  // P = (c - 0.5*c*c ) * dI2b
671  //
672  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
673  ie.SetJacobian(Jpt.GetData());
674  const double I2b = ie.Get_I2b();
675  const double c = (I2b - 1.0)/(I2b - tau0);
676  P.Set(c - 0.5*c*c, ie.Get_dI2b());
677 }
678 
679 void TMOP_Metric_252::AssembleH(const DenseMatrix &Jpt,
680  const DenseMatrix &DS,
681  const double weight,
682  DenseMatrix &A) const
683 {
684  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
685  //
686  // P = (c - 0.5*c*c ) * dI2b
687  // dP = (1 - c)^2/(I2b - tau0)*(dI2b x dI2b) + (c - 0.5*c*c)*ddI2b
688  ie.SetJacobian(Jpt.GetData());
689  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
690  const double I2b = ie.Get_I2b();
691  const double c0 = 1.0/(I2b - tau0);
692  const double c = c0*(I2b - 1.0);
693  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI2b(), A.GetData());
694  ie.Assemble_ddI2b(weight*(c - 0.5*c*c), A.GetData());
695 }
696 
697 double TMOP_Metric_301::EvalWMatrixForm(const DenseMatrix &Jpt) const
698 {
699  // mu_301 = 1/3 |J| |J^-1| - 1.
700  ie.SetJacobian(Jpt.GetData());
701  DenseMatrix inv(3);
702  CalcInverse(Jpt, inv);
703  return Jpt.FNorm() * inv.FNorm() / 3.0 - 1.0;
704 }
705 
706 double TMOP_Metric_301::EvalW(const DenseMatrix &Jpt) const
707 {
708  // mu_301 = 1/3 sqrt(I1b * I2b) - 1
709  ie.SetJacobian(Jpt.GetData());
710  return std::sqrt(ie.Get_I1b()*ie.Get_I2b())/3. - 1.;
711 }
712 
713 void TMOP_Metric_301::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
714 {
715  // W = (1/3)*sqrt(I1b*I2b) - 1
716  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
717  ie.SetJacobian(Jpt.GetData());
718  const double a = 1./(6.*std::sqrt(ie.Get_I1b()*ie.Get_I2b()));
719  Add(a*ie.Get_I2b(), ie.Get_dI1b(), a*ie.Get_I1b(), ie.Get_dI2b(), P);
720 }
721 
722 void TMOP_Metric_301::AssembleH(const DenseMatrix &Jpt,
723  const DenseMatrix &DS,
724  const double weight,
725  DenseMatrix &A) const
726 {
727  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
728  // dW = (1/6)*[z2*dI1b + z1*dI2b], z1 = sqrt(I1b/I2b), z2 = sqrt(I2b/I1b)
729  // ddW = (1/6)*[dI1b x dz2 + z2*ddI1b + dI2b x dz1 + z1*ddI2b]
730  //
731  // dz1 = (1/2)*sqrt(I2b/I1b) [ (1/I2b)*dI1b + (I1b/(I2b*I2b))*dI2b ]
732  // = (1/2)/sqrt(I1b*I2b) [ dI1b + (I1b/I2b)*dI2b ]
733  // dz2 = (1/2)/sqrt(I1b*I2b) [ (I2b/I1b)*dI1b + dI2b ]
734  //
735  // dI1b x dz2 + dI2b x dz1 =
736  // (1/2)/sqrt(I1b*I2b) dI1b x [ (I2b/I1b)*dI1b + dI2b ] +
737  // (1/2)/sqrt(I1b*I2b) dI2b x [ dI1b + (I1b/I2b)*dI2b ] =
738  // (1/2)/sqrt(I1b*I2b) [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] x
739  // [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] =
740  // (1/2)/sqrt(I1b*I2b) [ 6*dW x 6*dW ] =
741  // (1/2)*(I1b*I2b)^{-3/2} (I2b*dI1b + I1b*dI2b) x (I2b*dI1b + I1b*dI2b)
742  //
743  // z1 = I1b/sqrt(I1b*I2b), z2 = I2b/sqrt(I1b*I2b)
744 
745  ie.SetJacobian(Jpt.GetData());
746  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
747  double d_I1b_I2b_data[9];
748  DenseMatrix d_I1b_I2b(d_I1b_I2b_data, 3, 3);
749  Add(ie.Get_I2b(), ie.Get_dI1b(), ie.Get_I1b(), ie.Get_dI2b(), d_I1b_I2b);
750  const double I1b_I2b = ie.Get_I1b()*ie.Get_I2b();
751  const double a = weight/(6*std::sqrt(I1b_I2b));
752  ie.Assemble_ddI1b(a*ie.Get_I2b(), A.GetData());
753  ie.Assemble_ddI2b(a*ie.Get_I1b(), A.GetData());
754  ie.Assemble_TProd(a/(2*I1b_I2b), d_I1b_I2b_data, A.GetData());
755 }
756 
757 double TMOP_Metric_302::EvalWMatrixForm(const DenseMatrix &Jpt) const
758 {
759  // mu_301 = |J|^2 |J^{-1}|^2 / 9 - 1.
760  ie.SetJacobian(Jpt.GetData());
761  DenseMatrix inv(3);
762  CalcInverse(Jpt, inv);
763  return Jpt.FNorm2() * inv.FNorm2() / 9.0 - 1.0;
764 }
765 
766 double TMOP_Metric_302::EvalW(const DenseMatrix &Jpt) const
767 {
768  // mu_2 = |J|^2 |J^{-1}|^2 / 9 - 1
769  // = (l1^2 + l2^2 + l3^3)*(l1^{-2} + l2^{-2} + l3^{-2}) / 9 - 1
770  // = I1*(l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/l1^2/l2^2/l3^2/9 - 1
771  // = I1*I2/det(J)^2/9 - 1 = I1b*I2b/9-1
772  ie.SetJacobian(Jpt.GetData());
773  return ie.Get_I1b()*ie.Get_I2b()/9. - 1.;
774 }
775 
776 void TMOP_Metric_302::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
777 {
778  // mu_2 = I1b*I2b/9-1
779  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
780  ie.SetJacobian(Jpt.GetData());
781  Add(ie.Get_I1b()/9, ie.Get_dI2b(), ie.Get_I2b()/9, ie.Get_dI1b(), P);
782 }
783 
784 void TMOP_Metric_302::AssembleH(const DenseMatrix &Jpt,
785  const DenseMatrix &DS,
786  const double weight,
787  DenseMatrix &A) const
788 {
789  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
790  // dP = (dI2b x dI1b)/9 + (I1b/9)*ddI2b + (dI1b x dI2b)/9 + (I2b/9)*ddI1b
791  // = (dI2b x dI1b + dI1b x dI2b)/9 + (I1b/9)*ddI2b + (I2b/9)*ddI1b
792  ie.SetJacobian(Jpt.GetData());
793  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
794  const double c1 = weight/9;
795  ie.Assemble_TProd(c1, ie.Get_dI1b(), ie.Get_dI2b(), A.GetData());
796  ie.Assemble_ddI2b(c1*ie.Get_I1b(), A.GetData());
797  ie.Assemble_ddI1b(c1*ie.Get_I2b(), A.GetData());
798 }
799 
800 double TMOP_Metric_303::EvalWMatrixForm(const DenseMatrix &Jpt) const
801 {
802  // mu_303 = |J|^2 / 3 / det(J)^(2/3) - 1.
803  ie.SetJacobian(Jpt.GetData());
804  return Jpt.FNorm2() / 3.0 / pow(Jpt.Det(), 2.0 / 3.0) - 1.0;
805 }
806 
807 double TMOP_Metric_303::EvalW(const DenseMatrix &Jpt) const
808 {
809  // mu_303 = |J|^2 / 3 / det(J)^(2/3) - 1 = I1b/3 - 1.
810  ie.SetJacobian(Jpt.GetData());
811  return ie.Get_I1b()/3.0 - 1.0;
812 }
813 
814 void TMOP_Metric_303::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
815 {
816  // mu_304 = I1b/3 - 1.
817  // P = dI1b/3.
818  ie.SetJacobian(Jpt.GetData());
819  P.Set(1./3., ie.Get_dI1b());
820 }
821 
822 void TMOP_Metric_303::AssembleH(const DenseMatrix &Jpt,
823  const DenseMatrix &DS,
824  const double weight,
825  DenseMatrix &A) const
826 {
827  // P = dI1b/3.
828  // dP = ddI1b/3.
829  ie.SetJacobian(Jpt.GetData());
830  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
831  ie.Assemble_ddI1b(weight/3., A.GetData());
832 }
833 
834 double TMOP_Metric_304::EvalWMatrixForm(const DenseMatrix &Jpt) const
835 {
836  // mu_304 = |J|^3 / 3^(3/2) / det(J) - 1
837  const double fnorm = Jpt.FNorm();
838  return fnorm * fnorm * fnorm / pow(3.0, 1.5) / Jpt.Det() - 1.0;
839 }
840 
841 double TMOP_Metric_304::EvalW(const DenseMatrix &Jpt) const
842 {
843  // mu_304 = (I1b/3)^3/2 - 1.
844  ie.SetJacobian(Jpt.GetData());
845  return pow(ie.Get_I1b()/3.0, 1.5) - 1.0;
846 }
847 
848 void TMOP_Metric_304::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
849 {
850  // mu_304 = (I1b/3)^3/2 - 1.
851  // P = 3/2 * (I1b/3)^1/2 * dI1b / 3 = 1/2 * (I1b/3)^1/2 * dI1b.
852  ie.SetJacobian(Jpt.GetData());
853  P.Set(0.5 * sqrt(ie.Get_I1b()/3.0), ie.Get_dI1b());
854 }
855 
856 void TMOP_Metric_304::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
857  const double weight, DenseMatrix &A) const
858 {
859  // P = 1/2 * (I1b/3)^1/2 * dI1b.
860  // dP = 1/12 * (I1b/3)^(-1/2) * (dI1b x dI1b) + 1/2 * (I1b/3)^1/2 * ddI1b.
861  ie.SetJacobian(Jpt.GetData());
862  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
863  ie.Assemble_TProd(weight / 12.0 / sqrt(ie.Get_I1b()/3.0),
864  ie.Get_dI1b(), A.GetData());
865  ie.Assemble_ddI1b(weight / 2.0 * sqrt(ie.Get_I1b()/3.0), A.GetData());
866 }
867 
868 double TMOP_Metric_311::EvalW(const DenseMatrix &Jpt) const
869 {
870  // mu_311 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
871  // = (I3b - 1)^2 - I3b + sqrt(I3b^2 + eps)
872  ie.SetJacobian(Jpt.GetData());
873  const double I3b = ie.Get_I3b();
874  return (I3b - 1.0)*(I3b - 1.0) - I3b + std::sqrt(I3b*I3b + eps);
875 }
876 
877 void TMOP_Metric_311::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
878 {
879  ie.SetJacobian(Jpt.GetData());
880  const double I3b = ie.Get_I3b();
881  const double c = 2*I3b-3+(I3b)/(std::pow((I3b*I3b+eps),0.5));
882  P.Set(c, ie.Get_dI3b());
883 }
884 
885 void TMOP_Metric_311::AssembleH(const DenseMatrix &Jpt,
886  const DenseMatrix &DS,
887  const double weight,
888  DenseMatrix &A) const
889 {
890  ie.SetJacobian(Jpt.GetData());
891  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
892  const double I3b = ie.Get_I3b();
893  const double c0 = I3b*I3b+eps;
894  const double c1 = 2 + 1/(pow(c0,0.5)) - I3b*I3b/(pow(c0,1.5));
895  const double c2 = 2*I3b - 3 + I3b/(pow(c0,0.5));
896  ie.Assemble_TProd(weight*c1, ie.Get_dI3b(), A.GetData());
897  ie.Assemble_ddI3b(c2*weight, A.GetData());
898 }
899 
900 double TMOP_Metric_313::EvalW(const DenseMatrix &Jpt) const
901 {
902  ie.SetJacobian(Jpt.GetData());
903 
904  const double I3b = ie.Get_I3b();
905  double d = I3b - min_detT;
906  if (d < 0.0 && min_detT == 0.0)
907  {
908  // The mesh has been untangled, but it's still possible to get negative
909  // detJ in FD calculations, as they move the nodes around with some small
910  // increments and can produce negative determinants. Thus we put a small
911  // value in the denominator. Note that here I3b < 0.
912  d = - I3b * 0.1;
913  }
914 
915  const double c = std::pow(d, -2.0/3.0);
916 
917  return ie.Get_I1() * c / 3.0;
918 }
919 
920 void TMOP_Metric_313::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
921 {
922  MFEM_ABORT("Metric not implemented yet.");
923 }
924 
925 void TMOP_Metric_313::AssembleH(const DenseMatrix &Jpt,
926  const DenseMatrix &DS,
927  const double weight,
928  DenseMatrix &A) const
929 {
930  MFEM_ABORT("Metric not implemented yet.");
931 }
932 
933 double TMOP_Metric_315::EvalW(const DenseMatrix &Jpt) const
934 {
935  // mu_315 = mu_15_3D = (det(J) - 1)^2
936  ie.SetJacobian(Jpt.GetData());
937  const double c1 = ie.Get_I3b() - 1.0;
938  return c1*c1;
939 }
940 
941 void TMOP_Metric_315::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
942 {
943  // mu_315 = (I3b - 1)^2
944  // P = 2*(I3b - 1)*dI3b
945  ie.SetJacobian(Jpt.GetData());
946  P.Set(2*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
947 }
948 
949 void TMOP_Metric_315::AssembleH(const DenseMatrix &Jpt,
950  const DenseMatrix &DS,
951  const double weight,
952  DenseMatrix &A) const
953 {
954  // P = 2*(I3b - 1)*dI3b
955  // dP = 2*(dI3b x dI3b) + 2*(I3b - 1)*ddI3b
956  ie.SetJacobian(Jpt.GetData());
957  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
958  ie.Assemble_TProd(2*weight, ie.Get_dI3b(), A.GetData());
959  ie.Assemble_ddI3b(2*weight*(ie.Get_I3b() - 1.0), A.GetData());
960 }
961 
962 double TMOP_Metric_316::EvalWMatrixForm(const DenseMatrix &Jpt) const
963 {
964  // mu_316 = 0.5 (det(J) + 1/det(J)) - 1.
965  return 0.5 * (Jpt.Det() + 1.0 / Jpt.Det()) - 1.0;
966 }
967 
968 double TMOP_Metric_316::EvalW(const DenseMatrix &Jpt) const
969 {
970  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
971  ie.SetJacobian(Jpt.GetData());
972  const double I3b = ie.Get_I3b();
973  return 0.5*(I3b + 1.0/I3b) - 1.0;
974 }
975 
976 void TMOP_Metric_316::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
977 {
978  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
979  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
980  ie.SetJacobian(Jpt.GetData());
981  P.Set(0.5 - 0.5/ie.Get_I3(), ie.Get_dI3b());
982 }
983 
984 void TMOP_Metric_316::AssembleH(const DenseMatrix &Jpt,
985  const DenseMatrix &DS,
986  const double weight,
987  DenseMatrix &A) const
988 {
989  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
990  // dP = (1/I3b^3)*(dI3b x dI3b) + (0.5 - 0.5/I3)*ddI3b
991  ie.SetJacobian(Jpt.GetData());
992  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
993  ie.Assemble_TProd(weight/(ie.Get_I3()*ie.Get_I3b()),
994  ie.Get_dI3b(), A.GetData());
995  ie.Assemble_ddI3b(weight*(0.5 - 0.5/ie.Get_I3()), A.GetData());
996 }
997 
998 double TMOP_Metric_321::EvalWMatrixForm(const DenseMatrix &Jpt) const
999 {
1000  // mu_321 = |J - J^-t|^2.
1001  ie.SetJacobian(Jpt.GetData());
1002  DenseMatrix invt(3);
1003  CalcInverseTranspose(Jpt, invt);
1004  invt.Add(-1.0, Jpt);
1005  return invt.FNorm2();
1006 }
1007 
1008 double TMOP_Metric_321::EvalW(const DenseMatrix &Jpt) const
1009 {
1010  // mu_321 = mu_21_3D = |J - J^{-t}|^2
1011  // = |J|^2 + |J^{-1}|^2 - 6
1012  // = |J|^2 + (l1^{-2} + l2^{-2} + l3^{-2}) - 6
1013  // = |J|^2 + (l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/det(J)^2 - 6
1014  // = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
1015  ie.SetJacobian(Jpt.GetData());
1016  return ie.Get_I1() + ie.Get_I2()/ie.Get_I3() - 6.0;
1017 }
1018 
1019 void TMOP_Metric_321::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1020 {
1021  // mu_321 = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
1022  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
1023  ie.SetJacobian(Jpt.GetData());
1024  const double I3 = ie.Get_I3();
1025  Add(1.0/I3, ie.Get_dI2(),
1026  -2*ie.Get_I2()/(I3*ie.Get_I3b()), ie.Get_dI3b(), P);
1027  P += ie.Get_dI1();
1028 }
1029 
1030 void TMOP_Metric_321::AssembleH(const DenseMatrix &Jpt,
1031  const DenseMatrix &DS,
1032  const double weight,
1033  DenseMatrix &A) const
1034 {
1035  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
1036  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b) + (1/I3)*ddI2 + (dI3b x dz) + z*ddI3b
1037  //
1038  // z = -2*I2/I3b^3
1039  // dz = (-2/I3b^3)*dI2 + (2*I2)*(3/I3b^4)*dI3b
1040  //
1041  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b + dI3b x dI2) + (1/I3)*ddI2
1042  // + (6*I2/I3b^4)*(dI3b x dI3b) + (-2*I2/I3b^3)*ddI3b
1043  ie.SetJacobian(Jpt.GetData());
1044  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1045  const double c0 = 1.0/ie.Get_I3b();
1046  const double c1 = weight*c0*c0;
1047  const double c2 = -2*c0*c1;
1048  const double c3 = c2*ie.Get_I2();
1049  ie.Assemble_ddI1(weight, A.GetData());
1050  ie.Assemble_ddI2(c1, A.GetData());
1051  ie.Assemble_ddI3b(c3, A.GetData());
1052  ie.Assemble_TProd(c2, ie.Get_dI2(), ie.Get_dI3b(), A.GetData());
1053  ie.Assemble_TProd(-3*c0*c3, ie.Get_dI3b(), A.GetData());
1054 }
1055 
1056 double TMOP_Metric_322::EvalWMatrixForm(const DenseMatrix &Jpt) const
1057 {
1058  // mu_322 = 1 / (6 det(J)) |J - adj(J)^t|^2
1059  DenseMatrix adj_J_t(3);
1060  CalcAdjugateTranspose(Jpt, adj_J_t);
1061  adj_J_t *= -1.0;
1062  adj_J_t.Add(1.0, Jpt);
1063  return 1.0 / 6.0 / Jpt.Det() * adj_J_t.FNorm2();
1064 }
1065 
1066 double TMOP_Metric_322::EvalW(const DenseMatrix &Jpt) const
1067 {
1068  // mu_322 = 1 / (6 det(J)) |J - adj(J)^t|^2
1069  // = 1 / (6 det(J)) |J|^2 + 1/6 det(J) |J^{-1}|^2 - 1
1070  // = I1b / (I3b^-1/3) / 6 + I2b (I3b^1/3) / 6 - 1
1071  ie.SetJacobian(Jpt.GetData());
1072 
1073  return ie.Get_I1b() / pow(ie.Get_I3b(), 1.0/3.0) / 6.0 +
1074  ie.Get_I2b() * pow(ie.Get_I3b(), 1.0/3.0) / 6.0 - 1.0;
1075 }
1076 
1077 void TMOP_Metric_322::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1078 {
1079  // mu_322 = I1b (I3b^-1/3) / 6 + I2b (I3b^1/3) / 6 - 1
1080  // P = 1/6 (I3b^-1/3) dI1b - 1/18 I1b (I3b^-4/3) dI3b
1081  // + 1/6 (I3b^1/3) dI2b + 1/18 I2b (I3b^-2/3) dI3b
1082  ie.SetJacobian(Jpt.GetData());
1083  P.Set(1.0/6.0 * pow(ie.Get_I3b(), -1.0/3.0),
1084  ie.Get_dI1b());
1085  P.Add(-1.0/18.0 * ie.Get_I1b() * pow(ie.Get_I3b(), -4.0/3.0),
1086  ie.Get_dI3b());
1087  P.Add(1.0/6.0 * pow(ie.Get_I3b(), 1.0/3.0),
1088  ie.Get_dI2b());
1089  P.Add(1.0/18.0 * ie.Get_I2b() * pow(ie.Get_I3b(), -2.0/3.0),
1090  ie.Get_dI3b());
1091 }
1092 
1093 void TMOP_Metric_322::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
1094  const double weight, DenseMatrix &A) const
1095 {
1096  // P = 1/6 (I3b^-1/3) dI1b - 1/18 I1b (I3b^-4/3) dI3b
1097  // + 1/6 (I3b^1/3) dI2b + 1/18 I2b (I3b^-2/3) dI3b
1098  // dP = 1/6 (I3b^-1/3) ddI1b - 1/18 (I3b^-4/3) (dI1b x dI3b)
1099  // - 1/18 I1b (I3b^-4/3) ddI3b
1100  // - 1/18 (I3b^-4/3) (dI3b x dI1b)
1101  // + 2/27 I1b (I3b^-7/3) (dI3b x dI3b)
1102  // + 1/6 (I3b^1/3) ddI2b + 1/18 (I3b^-2/3) (dI2b x dI3b)
1103  // + 1/18 I2b (I3b^-2/3) ddI3b
1104  // + 1/18 (I3b^-2/3) (dI3b x dI2b)
1105  // - 1/27 I2b (I3b^-5/3) (dI3b x dI3b)
1106  ie.SetJacobian(Jpt.GetData());
1107  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1108  const double p13 = weight * pow(ie.Get_I3b(), 1.0/3.0),
1109  m13 = weight * pow(ie.Get_I3b(), -1.0/3.0),
1110  m23 = weight * pow(ie.Get_I3b(), -2.0/3.0),
1111  m43 = weight * pow(ie.Get_I3b(), -4.0/3.0),
1112  m53 = weight * pow(ie.Get_I3b(), -5.0/3.0),
1113  m73 = weight * pow(ie.Get_I3b(), -7.0/3.0);
1114  ie.Assemble_ddI1b(1.0/6.0 * m13, A.GetData());
1115  // Combines - 1/18 (I3b^-4/3) (dI1b x dI3b) - 1/18 (I3b^-4/3) (dI3b x dI1b).
1116  ie.Assemble_TProd(-1.0/18.0 * m43,
1117  ie.Get_dI1b(), ie.Get_dI3b(), A.GetData());
1118  ie.Assemble_ddI3b(-1.0/18.0 * ie.Get_I1b() * m43, A.GetData());
1119  ie.Assemble_TProd(2.0/27.0 * ie.Get_I1b() * m73,
1120  ie.Get_dI3b(), A.GetData());
1121  ie.Assemble_ddI2b(1.0/6.0 * p13, A.GetData());
1122  // Combines + 1/18 (I3b^-2/3) (dI2b x dI3b) + 1/18 (I3b^-2/3) (dI3b x dI2b).
1123  ie.Assemble_TProd(1.0/18.0 * m23,
1124  ie.Get_dI2b(), ie.Get_dI3b(), A.GetData());
1125  ie.Assemble_ddI3b(1.0/18.0 * ie.Get_I2b() * m23, A.GetData());
1126  ie.Assemble_TProd(-1.0/27.0 * ie.Get_I2b() * m53,
1127  ie.Get_dI3b(), A.GetData());
1128 }
1129 
1130 double TMOP_Metric_323::EvalWMatrixForm(const DenseMatrix &Jpt) const
1131 {
1132  // mu_323 = |J|^3 - 3 sqrt(3) ln(det(J)) - 3 sqrt(3).
1133  double fnorm = Jpt.FNorm();
1134  return fnorm * fnorm * fnorm - 3.0 * sqrt(3.0) * (log(Jpt.Det()) + 1.0);
1135 }
1136 
1137 double TMOP_Metric_323::EvalW(const DenseMatrix &Jpt) const
1138 {
1139  // mu_323 = I1^3/2 - 3 sqrt(3) ln(I3b) - 3 sqrt(3).
1140  ie.SetJacobian(Jpt.GetData());
1141  return pow(ie.Get_I1(), 1.5) - 3.0 * sqrt(3.0) * (log(ie.Get_I3b()) + 1.0);
1142 }
1143 
1144 void TMOP_Metric_323::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1145 {
1146  // mu_323 = I1^3/2 - 3 sqrt(3) ln(I3b) - 3 sqrt(3).
1147  // P = 3/2 (I1^1/2) dI1 - 3 sqrt(3) (I3b^-1) dI3b.
1148  ie.SetJacobian(Jpt.GetData());
1149  P.Set(1.5 * sqrt(ie.Get_I1()), ie.Get_dI1());
1150  P.Add(- 3.0 * sqrt(3.0) / ie.Get_I3b(), ie.Get_dI3b());
1151 }
1152 
1153 void TMOP_Metric_323::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
1154  const double weight, DenseMatrix &A) const
1155 {
1156  // P = 3/2 (I1^1/2) dI1 - 3 sqrt(3) (I3b^-1) dI3b
1157  // dP = 3/2 (I1^1/2) ddI1 + 3/4 (I1^-1/2) (dI1 x dI1)
1158  // - 3 sqrt(3) (I3b^-1) ddI3b + 3 sqrt(3) (I3b^-2) (dI3b x dI3b)
1159  ie.SetJacobian(Jpt.GetData());
1160  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1161  ie.Assemble_ddI1(weight * 1.5 * sqrt(ie.Get_I1()), A.GetData());
1162  ie.Assemble_TProd(weight * 0.75 / sqrt(ie.Get_I1()),
1163  ie.Get_dI1(), A.GetData());
1164  ie.Assemble_ddI3b(- weight * 3.0 * sqrt(3.0) / ie.Get_I3b(), A.GetData());
1165  ie.Assemble_TProd(weight * 3.0 * sqrt(3.0) / ie.Get_I3b() / ie.Get_I3b(),
1166  ie.Get_dI3b(), A.GetData());
1167 }
1168 
1169 double TMOP_Metric_352::EvalW(const DenseMatrix &Jpt) const
1170 {
1171  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
1172  ie.SetJacobian(Jpt.GetData());
1173  const double I3b = ie.Get_I3b();
1174  return 0.5*(I3b - 1.0)*(I3b - 1.0)/(I3b - tau0);
1175 }
1176 
1177 void TMOP_Metric_352::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1178 {
1179  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
1180  // P = (I3b - 1)/(I3b - tau0)*dI3b + 0.5*(I3b - 1)^2*(-1/(I3b - tau0)^2)*dI3b
1181  // = [ (I3b - 1)/(I3b - tau0) - 0.5*(I3b - 1)^2/(I3b - tau0)^2 ] * dI3b
1182  // = (c - 0.5*c*c) * dI3b
1183  ie.SetJacobian(Jpt.GetData());
1184  const double I3b = ie.Get_I3b();
1185  const double c = (I3b - 1.0)/(I3b - tau0);
1186  P.Set(c - 0.5*c*c, ie.Get_dI3b());
1187 }
1188 
1189 void TMOP_Metric_352::AssembleH(const DenseMatrix &Jpt,
1190  const DenseMatrix &DS,
1191  const double weight,
1192  DenseMatrix &A) const
1193 {
1194  // c = (I3b - 1)/(I3b - tau0)
1195  //
1196  // P = (c - 0.5*c*c) * dI3b
1197  // dP = (1 - c)*(dI3b x dc) + (c - 0.5*c*c)*ddI3b
1198  //
1199  // dc = 1/(I3b - tau0)*dI3b - (I3b - 1)/(I3b - tau)^2*dI3b =
1200  // = (1 - c)/(I3b - tau0)*dI3b
1201  //
1202  // dP = (1 - c)^2/(I3b - tau0)*(dI3b x dI3b) + (c - 0.5*c*c)*ddI3b
1203  ie.SetJacobian(Jpt.GetData());
1204  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1205  const double I3b = ie.Get_I3b();
1206  const double c0 = 1.0/(I3b - tau0);
1207  const double c = c0*(I3b - 1.0);
1208  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI3b(), A.GetData());
1209  ie.Assemble_ddI3b(weight*(c - 0.5*c*c), A.GetData());
1210 }
1211 
1212 
1213 
1214 
1215 double TMOP_Metric_360::EvalWMatrixForm(const DenseMatrix &Jpt) const
1216 {
1217  // mu_360 = |J|^3 / 3^(3/2) - det(J)
1218  const double fnorm = Jpt.FNorm();
1219  return fnorm * fnorm * fnorm / pow(3.0, 1.5) - Jpt.Det();
1220 }
1221 
1222 double TMOP_Metric_360::EvalW(const DenseMatrix &Jpt) const
1223 {
1224  // mu_360 = (I1/3)^(3/2) - I3b.
1225  ie.SetJacobian(Jpt.GetData());
1226  return pow(ie.Get_I1()/3.0, 1.5) - ie.Get_I3b();
1227 }
1228 
1229 void TMOP_Metric_360::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1230 {
1231  // mu_360 = (I1/3)^(3/2) - I3b.
1232  // P = 3/2 * (I1/3)^1/2 * dI1 / 3 - dI3b
1233  // = 1/2 * (I1/3)^1/2 * dI1 - dI3b.
1234  ie.SetJacobian(Jpt.GetData());
1235  Add(0.5 * sqrt(ie.Get_I1()/3.0), ie.Get_dI1(), -1.0, ie.Get_dI3b(), P);
1236 }
1237 
1238 void TMOP_Metric_360::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
1239  const double weight, DenseMatrix &A) const
1240 {
1241  // P = 1/2 * (I1/3)^1/2 * dI1 - dI3b.
1242  // dP = 1/12 * (I1/3)^(-1/2) * (dI1 x dI1) + 1/2 * (I1/3)^1/2 * ddI1 - ddI3b
1243  ie.SetJacobian(Jpt.GetData());
1244  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1245  ie.Assemble_TProd(weight / 12.0 / sqrt(ie.Get_I1()/3.0),
1246  ie.Get_dI1(), A.GetData());
1247  ie.Assemble_ddI1(weight / 2.0 * sqrt(ie.Get_I1()/3.0), A.GetData());
1248  ie.Assemble_ddI3b(-weight, A.GetData());
1249 }
1250 
1251 double TMOP_AMetric_011::EvalW(const DenseMatrix &Jpt) const
1252 {
1253  MFEM_VERIFY(Jtr != NULL,
1254  "Requires a target Jacobian, use SetTargetJacobian().");
1255 
1256  int dim = Jpt.Size();
1257 
1258  DenseMatrix Jpr(dim, dim);
1259  Mult(Jpt, *Jtr, Jpr);
1260 
1261  double alpha = Jpr.Det(),
1262  omega = Jtr->Det();
1263 
1264  DenseMatrix AdjAt(dim), WtW(dim), WRK(dim), Jtrt(dim);
1265  CalcAdjugateTranspose(Jpr, AdjAt);
1266  Jtrt.Transpose(*Jtr);
1267  MultAAt(Jtrt, WtW);
1268  WtW *= 1./omega;
1269  Mult(AdjAt, WtW, WRK);
1270 
1271  WRK -= Jpr;
1272  WRK *= -1.;
1273 
1274  return (0.25/alpha)*WRK.FNorm2();
1275 }
1276 
1277 double TMOP_AMetric_014a::EvalW(const DenseMatrix &Jpt) const
1278 {
1279  MFEM_VERIFY(Jtr != NULL,
1280  "Requires a target Jacobian, use SetTargetJacobian().");
1281 
1282  int dim = Jpt.Size();
1283 
1284  DenseMatrix Jpr(dim, dim);
1285  Mult(Jpt, *Jtr, Jpr);
1286 
1287  double sqalpha = pow(Jpr.Det(), 0.5),
1288  sqomega = pow(Jtr->Det(), 0.5);
1289 
1290  return 0.5*pow(sqalpha/sqomega - sqomega/sqalpha, 2.);
1291 }
1292 
1293 double TMOP_AMetric_036::EvalW(const DenseMatrix &Jpt) const
1294 {
1295  MFEM_VERIFY(Jtr != NULL,
1296  "Requires a target Jacobian, use SetTargetJacobian().");
1297 
1298  int dim = Jpt.Size();
1299 
1300  DenseMatrix Jpr(dim, dim);
1301  Mult(Jpt, *Jtr, Jpr); // T*W = A
1302 
1303  double alpha = Jpr.Det(); // det(A)
1304  Jpr -= *Jtr; // A-W
1305 
1306  return (1./alpha)*(Jpr.FNorm2()); //(1/alpha)*(|A-W|^2)
1307 }
1308 
1309 double TMOP_AMetric_107a::EvalW(const DenseMatrix &Jpt) const
1310 {
1311  MFEM_VERIFY(Jtr != NULL,
1312  "Requires a target Jacobian, use SetTargetJacobian().");
1313 
1314  int dim = Jpt.Size();
1315 
1316  DenseMatrix Jpr(dim, dim);
1317  Mult(Jpt, *Jtr, Jpr);
1318 
1319  double alpha = Jpr.Det(),
1320  aw = Jpr.FNorm()/Jtr->FNorm();
1321 
1322  DenseMatrix W = *Jtr;
1323  W *= aw;
1324  Jpr -= W;
1325 
1326  return (0.5/alpha)*Jpr.FNorm2();
1327 }
1328 
1329 
1330 void TargetConstructor::ComputeAvgVolume() const
1331 {
1332  MFEM_VERIFY(nodes, "Nodes are not given!");
1333  MFEM_ASSERT(avg_volume == 0.0, "The average volume is already computed!");
1334 
1335  Mesh *mesh = nodes->FESpace()->GetMesh();
1336  const int NE = mesh->GetNE();
1337  IsoparametricTransformation Tr;
1338  double volume = 0.0;
1339 
1340  for (int i = 0; i < NE; i++)
1341  {
1342  mesh->GetElementTransformation(i, *nodes, &Tr);
1343  const IntegrationRule &ir =
1344  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr.OrderJ());
1345  for (int j = 0; j < ir.GetNPoints(); j++)
1346  {
1347  const IntegrationPoint &ip = ir.IntPoint(j);
1348  Tr.SetIntPoint(&ip);
1349  volume += ip.weight * Tr.Weight();
1350  }
1351  }
1352 
1353  NCMesh *ncmesh = mesh->ncmesh;
1354  if (Parallel() == false)
1355  {
1356  avg_volume = (ncmesh == NULL) ?
1357  volume / NE : volume / ncmesh->GetNumRootElements();
1358 
1359  }
1360 #ifdef MFEM_USE_MPI
1361  else
1362  {
1363  double area_NE[4];
1364  area_NE[0] = volume; area_NE[1] = NE;
1365  MPI_Allreduce(area_NE, area_NE + 2, 2, MPI_DOUBLE, MPI_SUM, comm);
1366  avg_volume = (ncmesh == NULL) ?
1367  area_NE[2] / area_NE[3] : area_NE[2] / ncmesh->GetNumRootElements();
1368  }
1369 #endif
1370 }
1371 
1372 void TargetConstructor::ComputeAllElementTargets_Fallback(
1373  const FiniteElementSpace &fes,
1374  const IntegrationRule &ir,
1375  const Vector &xe,
1376  DenseTensor &Jtr) const
1377 {
1378  // Fallback to the 1-element method, ComputeElementTargets()
1379 
1380  // When UsesPhysicalCoordinates() == true, we assume 'xe' uses
1381  // ElementDofOrdering::LEXICOGRAPHIC iff 'fe' is a TensorFiniteElement.
1382 
1383  const Mesh *mesh = fes.GetMesh();
1384  const int NE = mesh->GetNE();
1385  // Quick return for empty processors:
1386  if (NE == 0) { return; }
1387  const int dim = mesh->Dimension();
1388  MFEM_VERIFY(mesh->GetNumGeometries(dim) <= 1,
1389  "mixed meshes are not supported");
1390  MFEM_VERIFY(!fes.IsVariableOrder(), "variable orders are not supported");
1391  const FiniteElement &fe = *fes.GetFE(0);
1392  const int sdim = fes.GetVDim();
1393  const int nvdofs = sdim*fe.GetDof();
1394  MFEM_VERIFY(!UsesPhysicalCoordinates() ||
1395  xe.Size() == NE*nvdofs, "invalid input Vector 'xe'!");
1396  const int NQ = ir.GetNPoints();
1397  const Array<int> *dof_map = nullptr;
1398  if (UsesPhysicalCoordinates())
1399  {
1400  const TensorBasisElement *tfe =
1401  dynamic_cast<const TensorBasisElement *>(&fe);
1402  if (tfe)
1403  {
1404  dof_map = &tfe->GetDofMap();
1405  if (dof_map->Size() == 0) { dof_map = nullptr; }
1406  }
1407  }
1408 
1409  Vector elfun_lex, elfun_nat;
1410  DenseTensor J;
1411  xe.HostRead();
1412  Jtr.HostWrite();
1413  if (UsesPhysicalCoordinates() && dof_map != nullptr)
1414  {
1415  elfun_nat.SetSize(nvdofs);
1416  }
1417  for (int e = 0; e < NE; e++)
1418  {
1419  if (UsesPhysicalCoordinates())
1420  {
1421  if (!dof_map)
1422  {
1423  elfun_nat.SetDataAndSize(xe.GetData()+e*nvdofs, nvdofs);
1424  }
1425  else
1426  {
1427  elfun_lex.SetDataAndSize(xe.GetData()+e*nvdofs, nvdofs);
1428  const int ndofs = fe.GetDof();
1429  for (int d = 0; d < sdim; d++)
1430  {
1431  for (int i_lex = 0; i_lex < ndofs; i_lex++)
1432  {
1433  elfun_nat[(*dof_map)[i_lex]+d*ndofs] =
1434  elfun_lex[i_lex+d*ndofs];
1435  }
1436  }
1437  }
1438  }
1439  J.UseExternalData(Jtr(e*NQ).Data(), sdim, dim, NQ);
1440  ComputeElementTargets(e, fe, ir, elfun_nat, J);
1441  }
1442 }
1443 
1444 bool TargetConstructor::ContainsVolumeInfo() const
1445 {
1446  switch (target_type)
1447  {
1448  case IDEAL_SHAPE_UNIT_SIZE: return false;
1449  case IDEAL_SHAPE_EQUAL_SIZE:
1450  case IDEAL_SHAPE_GIVEN_SIZE:
1451  case GIVEN_SHAPE_AND_SIZE:
1452  case GIVEN_FULL: return true;
1453  default: MFEM_ABORT("TargetType not added to ContainsVolumeInfo.");
1454  }
1455  return false;
1456 }
1457 
1458 void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
1459  const IntegrationRule &ir,
1460  const Vector &elfun,
1461  DenseTensor &Jtr) const
1462 {
1463  MFEM_CONTRACT_VAR(elfun);
1464  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1465 
1466  const FiniteElement *nfe = (target_type != IDEAL_SHAPE_UNIT_SIZE) ?
1467  nodes->FESpace()->GetFE(e_id) : NULL;
1468  const DenseMatrix &Wideal =
1469  Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
1470  MFEM_ASSERT(Wideal.Height() == Jtr.SizeI(), "");
1471  MFEM_ASSERT(Wideal.Width() == Jtr.SizeJ(), "");
1472 
1473  switch (target_type)
1474  {
1475  case IDEAL_SHAPE_UNIT_SIZE:
1476  {
1477  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = Wideal; }
1478  break;
1479  }
1480  case IDEAL_SHAPE_EQUAL_SIZE:
1481  {
1482  if (avg_volume == 0.0) { ComputeAvgVolume(); }
1483  DenseMatrix W(Wideal.Height());
1484 
1485  NCMesh *ncmesh = nodes->FESpace()->GetMesh()->ncmesh;
1486  double el_volume = avg_volume;
1487  if (ncmesh)
1488  {
1489  el_volume = avg_volume / ncmesh->GetElementSizeReduction(e_id);
1490  }
1491 
1492  W.Set(std::pow(volume_scale * el_volume / Wideal.Det(),
1493  1./W.Height()), Wideal);
1494  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = W; }
1495  break;
1496  }
1497  case IDEAL_SHAPE_GIVEN_SIZE:
1498  case GIVEN_SHAPE_AND_SIZE:
1499  {
1500  const int dim = nfe->GetDim(), dof = nfe->GetDof();
1501  MFEM_ASSERT(dim == nodes->FESpace()->GetVDim(), "");
1502  DenseMatrix dshape(dof, dim), pos(dof, dim);
1503  Array<int> xdofs(dof * dim);
1504  Vector posV(pos.Data(), dof * dim);
1505  double detW;
1506 
1507  // always initialize detW to suppress a warning:
1508  detW = (target_type == IDEAL_SHAPE_GIVEN_SIZE) ? Wideal.Det() : 0.0;
1509  nodes->FESpace()->GetElementVDofs(e_id, xdofs);
1510  nodes->GetSubVector(xdofs, posV);
1511  for (int i = 0; i < ir.GetNPoints(); i++)
1512  {
1513  nfe->CalcDShape(ir.IntPoint(i), dshape);
1514  MultAtB(pos, dshape, Jtr(i));
1515  if (target_type == IDEAL_SHAPE_GIVEN_SIZE)
1516  {
1517  const double det = Jtr(i).Det();
1518  MFEM_VERIFY(det > 0.0, "The given mesh is inverted!");
1519  Jtr(i).Set(std::pow(det / detW, 1./dim), Wideal);
1520  }
1521  }
1522  break;
1523  }
1524  default:
1525  MFEM_ABORT("invalid target type!");
1526  }
1527 }
1528 
1529 void TargetConstructor::ComputeElementTargetsGradient(
1530  const IntegrationRule &ir,
1531  const Vector &elfun,
1532  IsoparametricTransformation &Tpr,
1533  DenseTensor &dJtr) const
1534 {
1535  MFEM_CONTRACT_VAR(elfun);
1536  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1537 
1538  // TODO: Compute derivative for targets with GIVEN_SHAPE or/and GIVEN_SIZE
1539  for (int i = 0; i < Tpr.GetFE()->GetDim()*ir.GetNPoints(); i++)
1540  { dJtr(i) = 0.; }
1541 }
1542 
1543 void AnalyticAdaptTC::SetAnalyticTargetSpec(Coefficient *sspec,
1544  VectorCoefficient *vspec,
1545  TMOPMatrixCoefficient *mspec)
1546 {
1547  scalar_tspec = sspec;
1548  vector_tspec = vspec;
1549  matrix_tspec = mspec;
1550 }
1551 
1552 void AnalyticAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
1553  const IntegrationRule &ir,
1554  const Vector &elfun,
1555  DenseTensor &Jtr) const
1556 {
1557  DenseMatrix point_mat;
1558  point_mat.UseExternalData(elfun.GetData(), fe.GetDof(), fe.GetDim());
1559 
1560  switch (target_type)
1561  {
1562  case GIVEN_FULL:
1563  {
1564  MFEM_VERIFY(matrix_tspec != NULL,
1565  "Target type GIVEN_FULL requires a MatrixCoefficient.");
1566 
1567  IsoparametricTransformation Tpr;
1568  Tpr.SetFE(&fe);
1569  Tpr.ElementNo = e_id;
1570  Tpr.ElementType = ElementTransformation::ELEMENT;
1571  Tpr.GetPointMat().Transpose(point_mat);
1572 
1573  for (int i = 0; i < ir.GetNPoints(); i++)
1574  {
1575  const IntegrationPoint &ip = ir.IntPoint(i);
1576  Tpr.SetIntPoint(&ip);
1577  matrix_tspec->Eval(Jtr(i), Tpr, ip);
1578  }
1579  break;
1580  }
1581  default:
1582  MFEM_ABORT("Incompatible target type for analytic adaptation!");
1583  }
1584 }
1585 
1586 void AnalyticAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
1587  const Vector &elfun,
1588  IsoparametricTransformation &Tpr,
1589  DenseTensor &dJtr) const
1590 {
1591  const FiniteElement *fe = Tpr.GetFE();
1592  DenseMatrix point_mat;
1593  point_mat.UseExternalData(elfun.GetData(), fe->GetDof(), fe->GetDim());
1594 
1595  switch (target_type)
1596  {
1597  case GIVEN_FULL:
1598  {
1599  MFEM_VERIFY(matrix_tspec != NULL,
1600  "Target type GIVEN_FULL requires a TMOPMatrixCoefficient.");
1601 
1602  for (int d = 0; d < fe->GetDim(); d++)
1603  {
1604  for (int i = 0; i < ir.GetNPoints(); i++)
1605  {
1606  const IntegrationPoint &ip = ir.IntPoint(i);
1607  Tpr.SetIntPoint(&ip);
1608  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1609  matrix_tspec->EvalGrad(dJtr_i, Tpr, ip, d);
1610  }
1611  }
1612  break;
1613  }
1614  default:
1615  MFEM_ABORT("Incompatible target type for analytic adaptation!");
1616  }
1617 }
1618 
1619 namespace internal
1620 {
1621 
1622 // MFEM_FORALL-based copy kernel -- used by protected methods below.
1623 // Needed as a workaround for the nvcc restriction that methods with MFEM_FORALL
1624 // in them must to be public.
1625 static inline void device_copy(double *d_dest, const double *d_src, int size)
1626 {
1627  MFEM_FORALL(i, size, d_dest[i] = d_src[i];);
1628 }
1629 
1630 } // namespace internal
1631 
1632 #ifdef MFEM_USE_MPI
1633 void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction &t)
1634 {
1635  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1636  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1637 
1638  ParFiniteElementSpace *ptspec_fes = t.ParFESpace();
1639 
1640  tspec_sav = tspec;
1641 
1642  delete tspec_fesv;
1643  tspec_fesv = new FiniteElementSpace(ptspec_fes->GetMesh(),
1644  ptspec_fes->FEColl(), ncomp);
1645 
1646  delete ptspec_fesv;
1647  ptspec_fesv = new ParFiniteElementSpace(ptspec_fes->GetParMesh(),
1648  ptspec_fes->FEColl(), ncomp);
1649 
1650  delete tspec_pgf;
1651  tspec_pgf = new ParGridFunction(ptspec_fesv, tspec);
1652  tspec_gf = tspec_pgf;
1653 
1654  adapt_eval->SetParMetaInfo(*ptspec_fes->GetParMesh(), *ptspec_fesv);
1655  adapt_eval->SetInitialField(*ptspec_fes->GetMesh()->GetNodes(), tspec);
1656 }
1657 
1658 void DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange()
1659 {
1660  ptspec_fesv->Update();
1661  if (tspec_fesv)
1662  {
1663  delete tspec_fesv;
1664  tspec_fesv = new FiniteElementSpace(ptspec_fesv->GetMesh(),
1665  ptspec_fesv->FEColl(), ncomp);
1666  }
1667  tspec_pgf->Update();
1668  tspec_gf = tspec_pgf;
1669  tspec.SetDataAndSize(tspec_pgf->GetData(), tspec_pgf->Size());
1670  tspec_sav = tspec;
1671 
1672  adapt_eval->SetParMetaInfo(*ptspec_fesv->GetParMesh(), *ptspec_fesv);
1673  adapt_eval->SetInitialField(*ptspec_fesv->GetMesh()->GetNodes(), tspec);
1674 }
1675 
1676 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const ParGridFunction &tspec_)
1677 {
1678  const int vdim = tspec_.FESpace()->GetVDim(),
1679  ndof = tspec_.FESpace()->GetNDofs();
1680  MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
1681 
1682  const auto tspec__d = tspec_.Read();
1683  auto tspec_d = tspec.ReadWrite();
1684  const int offset = idx*ndof;
1685  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1686  FinalizeParDiscreteTargetSpec(tspec_);
1687 }
1688 
1689 void DiscreteAdaptTC::SetParDiscreteTargetSize(const ParGridFunction &tspec_)
1690 {
1691  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1692  "Discrete target size should be ordered byNodes.");
1693  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1694  sizeidx = ncomp;
1695  SetDiscreteTargetBase(tspec_);
1696  FinalizeParDiscreteTargetSpec(tspec_);
1697 }
1698 
1699 void DiscreteAdaptTC::SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
1700 {
1701  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1702  "Discrete target skewness should be ordered byNodes.");
1703  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1704  skewidx = ncomp;
1705  SetDiscreteTargetBase(tspec_);
1706  FinalizeParDiscreteTargetSpec(tspec_);
1707 }
1708 
1709 void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction &ar)
1710 {
1711  MFEM_VERIFY(ar.FESpace()->GetOrdering() == Ordering::byNODES,
1712  "Discrete target aspect ratio should be ordered byNodes.");
1713  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
1714  aspectratioidx = ncomp;
1715  SetDiscreteTargetBase(ar);
1716  FinalizeParDiscreteTargetSpec(ar);
1717 }
1718 
1719 void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction &o)
1720 {
1721  MFEM_VERIFY(o.FESpace()->GetOrdering() == Ordering::byNODES,
1722  "Discrete target orientation should be ordered byNodes.");
1723  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
1724  orientationidx = ncomp;
1725  SetDiscreteTargetBase(o);
1726  FinalizeParDiscreteTargetSpec(o);
1727 }
1728 
1729 void DiscreteAdaptTC::SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
1730 {
1731  SetParDiscreteTargetSize(tspec_);
1732 }
1733 #endif // MFEM_USE_MPI
1734 
1735 void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
1736 {
1737  const int vdim = tspec_.FESpace()->GetVDim(),
1738  ndof = tspec_.FESpace()->GetNDofs();
1739 
1740  ncomp += vdim;
1741 
1742  // need to append data to tspec
1743  // make a copy of tspec->tspec_temp, increase its size, and
1744  // copy data from tspec_temp -> tspec, then add new entries
1745  Vector tspec_temp = tspec;
1746  tspec.UseDevice(true);
1747  tspec_sav.UseDevice(true);
1748  tspec.SetSize(ncomp*ndof);
1749 
1750  const auto tspec_temp_d = tspec_temp.Read();
1751  auto tspec_d = tspec.ReadWrite();
1752  internal::device_copy(tspec_d, tspec_temp_d, tspec_temp.Size());
1753 
1754  const auto tspec__d = tspec_.Read();
1755  const int offset = (ncomp-vdim)*ndof;
1756  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1757 }
1758 
1759 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
1760 {
1761  const int vdim = tspec_.FESpace()->GetVDim(),
1762  ndof = tspec_.FESpace()->GetNDofs();
1763  MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTargetSpec.");
1764 
1765  const auto tspec__d = tspec_.Read();
1766  auto tspec_d = tspec.ReadWrite();
1767  const int offset = idx*ndof;
1768  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1769  FinalizeSerialDiscreteTargetSpec(tspec_);
1770 }
1771 
1772 void DiscreteAdaptTC::SetSerialDiscreteTargetSize(const GridFunction &tspec_)
1773 {
1774  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1775  "Discrete target size should be ordered byNodes.");
1776  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1777  sizeidx = ncomp;
1778  SetDiscreteTargetBase(tspec_);
1779  FinalizeSerialDiscreteTargetSpec(tspec_);
1780 }
1781 
1782 void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
1783 {
1784  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1785  "Discrete target skewness should be ordered byNodes.");
1786  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1787  skewidx = ncomp;
1788  SetDiscreteTargetBase(tspec_);
1789  FinalizeSerialDiscreteTargetSpec(tspec_);
1790 }
1791 
1792 void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(const GridFunction &ar)
1793 {
1794  MFEM_VERIFY(ar.FESpace()->GetOrdering() == Ordering::byNODES,
1795  "Discrete target aspect ratio should be ordered byNodes.");
1796  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
1797  aspectratioidx = ncomp;
1798  SetDiscreteTargetBase(ar);
1799  FinalizeSerialDiscreteTargetSpec(ar);
1800 }
1801 
1802 void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(const GridFunction &o)
1803 {
1804  MFEM_VERIFY(o.FESpace()->GetOrdering() == Ordering::byNODES,
1805  "Discrete target orientation should be ordered byNodes.");
1806  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
1807  orientationidx = ncomp;
1808  SetDiscreteTargetBase(o);
1809  FinalizeSerialDiscreteTargetSpec(o);
1810 }
1811 
1812 void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec(const GridFunction &t)
1813 {
1814  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1815  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1816 
1817  const FiniteElementSpace *tspec_fes = t.FESpace();
1818 
1819  tspec_sav = tspec;
1820 
1821  delete tspec_fesv;
1822  tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
1823  tspec_fes->FEColl(), ncomp,
1824  Ordering::byNODES);
1825 
1826  delete tspec_gf;
1827  tspec_gf = new GridFunction(tspec_fesv, tspec);
1828 
1829  adapt_eval->SetSerialMetaInfo(*tspec_fes->GetMesh(), *tspec_fesv);
1830  adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
1831 }
1832 
1833 void DiscreteAdaptTC::GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
1834 {
1835  if (idx < 0) { return; }
1836  const int ndof = tspec_.FESpace()->GetNDofs(),
1837  vdim = tspec_.FESpace()->GetVDim();
1838  MFEM_VERIFY(ndof == tspec.Size()/ncomp,
1839  "Inconsistency in GetSerialDiscreteTargetSpec.");
1840 
1841  for (int i = 0; i < ndof*vdim; i++)
1842  {
1843  tspec_(i) = tspec(i + idx*ndof);
1844  }
1845 }
1846 
1847 void DiscreteAdaptTC::UpdateAfterMeshTopologyChange()
1848 {
1849  tspec_fesv->Update();
1850  tspec_gf->Update();
1851  tspec.SetDataAndSize(tspec_gf->GetData(), tspec_gf->Size());
1852  tspec_sav = tspec;
1853 
1854  adapt_eval->SetSerialMetaInfo(*tspec_fesv->GetMesh(), *tspec_fesv);
1855  adapt_eval->SetInitialField(*tspec_fesv->GetMesh()->GetNodes(), tspec);
1856 }
1857 
1858 void DiscreteAdaptTC::SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
1859 {
1860  SetSerialDiscreteTargetSize(tspec_);
1861 }
1862 
1863 
1864 void DiscreteAdaptTC::UpdateTargetSpecification(const Vector &new_x,
1865  bool use_flag,
1866  int new_x_ordering)
1867 {
1868  if (use_flag && good_tspec) { return; }
1869 
1870  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1871  adapt_eval->ComputeAtNewPosition(new_x, tspec, new_x_ordering);
1872  tspec_sav = tspec;
1873 
1874  good_tspec = use_flag;
1875 }
1876 
1877 void DiscreteAdaptTC::UpdateTargetSpecification(Vector &new_x,
1878  Vector &IntData,
1879  int new_x_ordering)
1880 {
1881  adapt_eval->ComputeAtNewPosition(new_x, IntData, new_x_ordering);
1882 }
1883 
1884 void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
1885  ElementTransformation &T,
1886  int dofidx, int dir,
1887  const Vector &IntData)
1888 {
1889  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1890 
1891  Array<int> dofs;
1892  tspec_fesv->GetElementDofs(T.ElementNo, dofs);
1893  const int cnt = tspec.Size()/ncomp; // dofs per scalar-field
1894 
1895  for (int i = 0; i < ncomp; i++)
1896  {
1897  tspec(dofs[dofidx]+i*cnt) = IntData(dofs[dofidx] + i*cnt + dir*cnt*ncomp);
1898  }
1899 }
1900 
1901 void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
1902  int dofidx)
1903 {
1904  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1905 
1906  Array<int> dofs;
1907  tspec_fesv->GetElementDofs(T.ElementNo, dofs);
1908  const int cnt = tspec.Size()/ncomp;
1909  for (int i = 0; i < ncomp; i++)
1910  {
1911  tspec(dofs[dofidx] + i*cnt) = tspec_sav(dofs[dofidx] + i*cnt);
1912  }
1913 }
1914 
1915 void DiscreteAdaptTC::SetTspecFromIntRule(int e_id,
1916  const IntegrationRule &intrule)
1917 {
1918  switch (target_type)
1919  {
1920  case IDEAL_SHAPE_GIVEN_SIZE:
1921  case GIVEN_SHAPE_AND_SIZE:
1922  {
1923  const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
1924  ntspec_dofs = ndofs*ncomp;
1925 
1926  Vector tspec_vals(ntspec_dofs);
1927 
1928  Array<int> dofs;
1929  tspec_fesv->GetElementVDofs(e_id, dofs);
1930  tspec.GetSubVector(dofs, tspec_vals);
1931  DenseMatrix tr;
1932  tspec_gf->GetVectorValues(e_id, intrule, tspec_refine, tr);
1933  tspec_refine.Transpose();
1934  break;
1935  }
1936  default:
1937  MFEM_ABORT("Incompatible target type for discrete adaptation!");
1938  }
1939 }
1940 
1941 void DiscreteAdaptTC::SetTspecDataForDerefinement(FiniteElementSpace *fes)
1942 {
1943  coarse_tspec_fesv = fes;
1944  const Operator *c_op = fes->GetUpdateOperator();
1945  tspec_derefine.SetSize(c_op->Height());
1946  c_op->Mult(tspec, tspec_derefine);
1947 }
1948 
1949 void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
1950  const IntegrationRule &ir,
1951  const Vector &elfun,
1952  DenseTensor &Jtr) const
1953 {
1954  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
1955  const int dim = fe.GetDim(),
1956  nqp = ir.GetNPoints();
1957  Jtrcomp.SetSize(dim, dim, 4*nqp);
1958 
1959  FiniteElementSpace *src_fes = tspec_fesv;
1960 
1961  switch (target_type)
1962  {
1963  case IDEAL_SHAPE_GIVEN_SIZE:
1964  case GIVEN_SHAPE_AND_SIZE:
1965  {
1966  const DenseMatrix &Wideal =
1967  Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
1968  const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
1969  ntspec_dofs = ndofs*ncomp;
1970 
1971  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
1972  par_vals_c1, par_vals_c2, par_vals_c3;
1973 
1974  Array<int> dofs;
1975  DenseMatrix D_rho(dim), Q_phi(dim), R_theta(dim);
1976  tspec_fesv->GetElementVDofs(e_id, dofs);
1977  tspec.UseDevice(true);
1978  tspec.GetSubVector(dofs, tspec_vals);
1979  if (tspec_refine.NumCols() > 0) // Refinement
1980  {
1981  MFEM_VERIFY(amr_el >= 0, " Target being constructed for an AMR element.");
1982  for (int i = 0; i < ncomp; i++)
1983  {
1984  for (int j = 0; j < ndofs; j++)
1985  {
1986  tspec_vals(j + i*ndofs) = tspec_refine(j + amr_el*ndofs, i);
1987  }
1988  }
1989  }
1990  else if (tspec_derefine.Size() > 0) // Derefinement
1991  {
1992  dofs.SetSize(0);
1993  coarse_tspec_fesv->GetElementVDofs(e_id, dofs);
1994  tspec_derefine.GetSubVector(dofs, tspec_vals);
1995  src_fes = coarse_tspec_fesv;
1996  }
1997 
1998  for (int q = 0; q < nqp; q++)
1999  {
2000  const IntegrationPoint &ip = ir.IntPoint(q);
2001  src_fes->GetFE(e_id)->CalcShape(ip, shape);
2002  Jtr(q) = Wideal; // Initialize to identity
2003  for (int d = 0; d < 4; d++)
2004  {
2005  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(d + 4*q), dim, dim);
2006  Jtrcomp_q = Wideal; // Initialize to identity
2007  }
2008 
2009  if (sizeidx != -1) // Set size
2010  {
2011  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
2012  double min_size = par_vals.Min();//0.001; //
2013  if (lim_min_size > 0.)
2014  {
2015  min_size = lim_min_size;
2016  }
2017  else
2018  {
2019  MFEM_VERIFY(min_size > 0.0,
2020  "Non-positive size propagated in the target definition.");
2021  }
2022  const double size = std::max(shape * par_vals, min_size);
2023  Jtr(q).Set(std::pow(size, 1.0/dim), Jtr(q));
2024  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(0 + 4*q), dim, dim);
2025  Jtrcomp_q = Jtr(q);
2026  } // Done size
2027 
2028  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
2029 
2030  if (aspectratioidx != -1) // Set aspect ratio
2031  {
2032  if (dim == 2)
2033  {
2034  par_vals.SetDataAndSize(tspec_vals.GetData()+
2035  aspectratioidx*ndofs, ndofs);
2036  const double min_size = par_vals.Min();
2037  MFEM_VERIFY(min_size > 0.0,
2038  "Non-positive aspect-ratio propagated in the target definition.");
2039 
2040  const double aspectratio = shape * par_vals;
2041  D_rho = 0.;
2042  D_rho(0,0) = 1./pow(aspectratio,0.5);
2043  D_rho(1,1) = pow(aspectratio,0.5);
2044  }
2045  else
2046  {
2047  par_vals.SetDataAndSize(tspec_vals.GetData()+
2048  aspectratioidx*ndofs, ndofs*3);
2049  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
2050  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
2051  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
2052 
2053  const double rho1 = shape * par_vals_c1;
2054  const double rho2 = shape * par_vals_c2;
2055  const double rho3 = shape * par_vals_c3;
2056  D_rho = 0.;
2057  D_rho(0,0) = pow(rho1,2./3.);
2058  D_rho(1,1) = pow(rho2,2./3.);
2059  D_rho(2,2) = pow(rho3,2./3.);
2060  }
2061  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(1 + 4*q), dim, dim);
2062  Jtrcomp_q = D_rho;
2063  DenseMatrix Temp = Jtr(q);
2064  Mult(D_rho, Temp, Jtr(q));
2065  } // Done aspect ratio
2066 
2067  if (skewidx != -1) // Set skew
2068  {
2069  if (dim == 2)
2070  {
2071  par_vals.SetDataAndSize(tspec_vals.GetData()+
2072  skewidx*ndofs, ndofs);
2073 
2074  const double skew = shape * par_vals;
2075 
2076  Q_phi = 0.;
2077  Q_phi(0,0) = 1.;
2078  Q_phi(0,1) = cos(skew);
2079  Q_phi(1,1) = sin(skew);
2080  }
2081  else
2082  {
2083  par_vals.SetDataAndSize(tspec_vals.GetData()+
2084  skewidx*ndofs, ndofs*3);
2085  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
2086  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
2087  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
2088 
2089  const double phi12 = shape * par_vals_c1;
2090  const double phi13 = shape * par_vals_c2;
2091  const double chi = shape * par_vals_c3;
2092 
2093  Q_phi = 0.;
2094  Q_phi(0,0) = 1.;
2095  Q_phi(0,1) = cos(phi12);
2096  Q_phi(0,2) = cos(phi13);
2097 
2098  Q_phi(1,1) = sin(phi12);
2099  Q_phi(1,2) = sin(phi13)*cos(chi);
2100 
2101  Q_phi(2,2) = sin(phi13)*sin(chi);
2102  }
2103  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*q), dim, dim);
2104  Jtrcomp_q = Q_phi;
2105  DenseMatrix Temp = Jtr(q);
2106  Mult(Q_phi, Temp, Jtr(q));
2107  } // Done skew
2108 
2109  if (orientationidx != -1) // Set orientation
2110  {
2111  if (dim == 2)
2112  {
2113  par_vals.SetDataAndSize(tspec_vals.GetData()+
2114  orientationidx*ndofs, ndofs);
2115 
2116  const double theta = shape * par_vals;
2117  R_theta(0,0) = cos(theta);
2118  R_theta(0,1) = -sin(theta);
2119  R_theta(1,0) = sin(theta);
2120  R_theta(1,1) = cos(theta);
2121  }
2122  else
2123  {
2124  par_vals.SetDataAndSize(tspec_vals.GetData()+
2125  orientationidx*ndofs, ndofs*3);
2126  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
2127  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
2128  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
2129 
2130  const double theta = shape * par_vals_c1;
2131  const double psi = shape * par_vals_c2;
2132  const double beta = shape * par_vals_c3;
2133 
2134  double ct = cos(theta), st = sin(theta),
2135  cp = cos(psi), sp = sin(psi),
2136  cb = cos(beta), sb = sin(beta);
2137 
2138  R_theta = 0.;
2139  R_theta(0,0) = ct*sp;
2140  R_theta(1,0) = st*sp;
2141  R_theta(2,0) = cp;
2142 
2143  R_theta(0,1) = -st*cb + ct*cp*sb;
2144  R_theta(1,1) = ct*cb + st*cp*sb;
2145  R_theta(2,1) = -sp*sb;
2146 
2147  R_theta(0,0) = -st*sb - ct*cp*cb;
2148  R_theta(1,0) = ct*sb - st*cp*cb;
2149  R_theta(2,0) = sp*cb;
2150  }
2151  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(3 + 4*q), dim, dim);
2152  Jtrcomp_q = R_theta;
2153  DenseMatrix Temp = Jtr(q);
2154  Mult(R_theta, Temp, Jtr(q));
2155  } // Done orientation
2156  }
2157  break;
2158  }
2159  default:
2160  MFEM_ABORT("Incompatible target type for discrete adaptation!");
2161  }
2162 }
2163 
2164 void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
2165  const Vector &elfun,
2166  IsoparametricTransformation &Tpr,
2167  DenseTensor &dJtr) const
2168 {
2169  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
2170 
2171  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
2172 
2173  dJtr = 0.;
2174  const int e_id = Tpr.ElementNo;
2175  const FiniteElement *fe = Tpr.GetFE();
2176 
2177  switch (target_type)
2178  {
2179  case IDEAL_SHAPE_GIVEN_SIZE:
2180  case GIVEN_SHAPE_AND_SIZE:
2181  {
2182  const DenseMatrix &Wideal =
2183  Geometries.GetGeomToPerfGeomJac(fe->GetGeomType());
2184  const int dim = Wideal.Height(),
2185  ndofs = fe->GetDof(),
2186  ntspec_dofs = ndofs*ncomp;
2187 
2188  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
2189  par_vals_c1(ndofs), par_vals_c2(ndofs), par_vals_c3(ndofs);
2190 
2191  Array<int> dofs;
2192  DenseMatrix dD_rho(dim), dQ_phi(dim), dR_theta(dim);
2193  DenseMatrix dQ_phi13(dim), dQ_phichi(dim); // dQ_phi is used for dQ/dphi12 in 3D
2194  DenseMatrix dR_psi(dim), dR_beta(dim);
2195  tspec_fesv->GetElementVDofs(e_id, dofs);
2196  tspec.GetSubVector(dofs, tspec_vals);
2197 
2198  DenseMatrix grad_e_c1(ndofs, dim),
2199  grad_e_c2(ndofs, dim),
2200  grad_e_c3(ndofs, dim);
2201  Vector grad_ptr_c1(grad_e_c1.GetData(), ndofs*dim),
2202  grad_ptr_c2(grad_e_c2.GetData(), ndofs*dim),
2203  grad_ptr_c3(grad_e_c3.GetData(), ndofs*dim);
2204 
2205  DenseMatrix grad_phys; // This will be (dof x dim, dof).
2206  fe->ProjectGrad(*fe, Tpr, grad_phys);
2207 
2208  for (int i = 0; i < ir.GetNPoints(); i++)
2209  {
2210  const IntegrationPoint &ip = ir.IntPoint(i);
2211  DenseMatrix Jtrcomp_s(Jtrcomp.GetData(0 + 4*i), dim, dim); // size
2212  DenseMatrix Jtrcomp_d(Jtrcomp.GetData(1 + 4*i), dim, dim); // aspect-ratio
2213  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*i), dim, dim); // skew
2214  DenseMatrix Jtrcomp_r(Jtrcomp.GetData(3 + 4*i), dim, dim); // orientation
2215  DenseMatrix work1(dim), work2(dim), work3(dim);
2216 
2217  if (sizeidx != -1) // Set size
2218  {
2219  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
2220 
2221  grad_phys.Mult(par_vals, grad_ptr_c1);
2222  Vector grad_q(dim);
2223  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2224  grad_e_c1.MultTranspose(shape, grad_q);
2225 
2226  const double min_size = par_vals.Min();
2227  MFEM_VERIFY(min_size > 0.0,
2228  "Non-positive size propagated in the target definition.");
2229  const double size = std::max(shape * par_vals, min_size);
2230  double dz_dsize = (1./dim)*pow(size, 1./dim - 1.);
2231 
2232  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2233  Mult(Jtrcomp_r, work1, work2); // R*Q*D
2234 
2235  for (int d = 0; d < dim; d++)
2236  {
2237  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2238  work1 = Wideal;
2239  work1.Set(dz_dsize, work1); // dz/dsize
2240  work1 *= grad_q(d); // dz/dsize*dsize/dx
2241  AddMult(work1, work2, dJtr_i); // dz/dx*R*Q*D
2242  }
2243  } // Done size
2244 
2245  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
2246 
2247  if (aspectratioidx != -1) // Set aspect ratio
2248  {
2249  if (dim == 2)
2250  {
2251  par_vals.SetDataAndSize(tspec_vals.GetData()+
2252  aspectratioidx*ndofs, ndofs);
2253 
2254  grad_phys.Mult(par_vals, grad_ptr_c1);
2255  Vector grad_q(dim);
2256  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2257  grad_e_c1.MultTranspose(shape, grad_q);
2258 
2259  const double aspectratio = shape * par_vals;
2260  dD_rho = 0.;
2261  dD_rho(0,0) = -0.5*pow(aspectratio,-1.5);
2262  dD_rho(1,1) = 0.5*pow(aspectratio,-0.5);
2263 
2264  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
2265  Mult(work1, Jtrcomp_q, work2); // z*R*Q
2266 
2267  for (int d = 0; d < dim; d++)
2268  {
2269  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2270  work1 = dD_rho;
2271  work1 *= grad_q(d); // work1 = dD/drho*drho/dx
2272  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
2273  }
2274  }
2275  else // 3D
2276  {
2277  par_vals.SetDataAndSize(tspec_vals.GetData()+
2278  aspectratioidx*ndofs, ndofs*3);
2279  par_vals_c1.SetData(par_vals.GetData());
2280  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2281  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2282 
2283  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2284  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2285  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2286  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2287  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2288  grad_e_c1.MultTranspose(shape, grad_q1);
2289  grad_e_c2.MultTranspose(shape, grad_q2);
2290  grad_e_c3.MultTranspose(shape, grad_q3);
2291 
2292  const double rho1 = shape * par_vals_c1;
2293  const double rho2 = shape * par_vals_c2;
2294  const double rho3 = shape * par_vals_c3;
2295  dD_rho = 0.;
2296  dD_rho(0,0) = (2./3.)*pow(rho1,-1./3.);
2297  dD_rho(1,1) = (2./3.)*pow(rho2,-1./3.);
2298  dD_rho(2,2) = (2./3.)*pow(rho3,-1./3.);
2299 
2300  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
2301  Mult(work1, Jtrcomp_q, work2); // z*R*Q
2302 
2303 
2304  for (int d = 0; d < dim; d++)
2305  {
2306  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2307  work1 = dD_rho;
2308  work1(0,0) *= grad_q1(d);
2309  work1(1,2) *= grad_q2(d);
2310  work1(2,2) *= grad_q3(d);
2311  // work1 = dD/dx = dD/drho1*drho1/dx + dD/drho2*drho2/dx
2312  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
2313  }
2314  }
2315  } // Done aspect ratio
2316 
2317  if (skewidx != -1) // Set skew
2318  {
2319  if (dim == 2)
2320  {
2321  par_vals.SetDataAndSize(tspec_vals.GetData()+
2322  skewidx*ndofs, ndofs);
2323 
2324  grad_phys.Mult(par_vals, grad_ptr_c1);
2325  Vector grad_q(dim);
2326  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2327  grad_e_c1.MultTranspose(shape, grad_q);
2328 
2329  const double skew = shape * par_vals;
2330 
2331  dQ_phi = 0.;
2332  dQ_phi(0,0) = 1.;
2333  dQ_phi(0,1) = -sin(skew);
2334  dQ_phi(1,1) = cos(skew);
2335 
2336  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
2337 
2338  for (int d = 0; d < dim; d++)
2339  {
2340  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2341  work1 = dQ_phi;
2342  work1 *= grad_q(d); // work1 = dQ/dphi*dphi/dx
2343  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
2344  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
2345  }
2346  }
2347  else
2348  {
2349  par_vals.SetDataAndSize(tspec_vals.GetData()+
2350  skewidx*ndofs, ndofs*3);
2351  par_vals_c1.SetData(par_vals.GetData());
2352  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2353  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2354 
2355  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2356  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2357  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2358  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2359  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2360  grad_e_c1.MultTranspose(shape, grad_q1);
2361  grad_e_c2.MultTranspose(shape, grad_q2);
2362  grad_e_c3.MultTranspose(shape, grad_q3);
2363 
2364  const double phi12 = shape * par_vals_c1;
2365  const double phi13 = shape * par_vals_c2;
2366  const double chi = shape * par_vals_c3;
2367 
2368  dQ_phi = 0.;
2369  dQ_phi(0,0) = 1.;
2370  dQ_phi(0,1) = -sin(phi12);
2371  dQ_phi(1,1) = cos(phi12);
2372 
2373  dQ_phi13 = 0.;
2374  dQ_phi13(0,2) = -sin(phi13);
2375  dQ_phi13(1,2) = cos(phi13)*cos(chi);
2376  dQ_phi13(2,2) = cos(phi13)*sin(chi);
2377 
2378  dQ_phichi = 0.;
2379  dQ_phichi(1,2) = -sin(phi13)*sin(chi);
2380  dQ_phichi(2,2) = sin(phi13)*cos(chi);
2381 
2382  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
2383 
2384  for (int d = 0; d < dim; d++)
2385  {
2386  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2387  work1 = dQ_phi;
2388  work1 *= grad_q1(d); // work1 = dQ/dphi12*dphi12/dx
2389  work1.Add(grad_q2(d), dQ_phi13); // + dQ/dphi13*dphi13/dx
2390  work1.Add(grad_q3(d), dQ_phichi); // + dQ/dchi*dchi/dx
2391  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
2392  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
2393  }
2394  }
2395  } // Done skew
2396 
2397  if (orientationidx != -1) // Set orientation
2398  {
2399  if (dim == 2)
2400  {
2401  par_vals.SetDataAndSize(tspec_vals.GetData()+
2402  orientationidx*ndofs, ndofs);
2403 
2404  grad_phys.Mult(par_vals, grad_ptr_c1);
2405  Vector grad_q(dim);
2406  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2407  grad_e_c1.MultTranspose(shape, grad_q);
2408 
2409  const double theta = shape * par_vals;
2410  dR_theta(0,0) = -sin(theta);
2411  dR_theta(0,1) = -cos(theta);
2412  dR_theta(1,0) = cos(theta);
2413  dR_theta(1,1) = -sin(theta);
2414 
2415  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2416  Mult(Jtrcomp_s, work1, work2); // z*Q*D
2417  for (int d = 0; d < dim; d++)
2418  {
2419  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2420  work1 = dR_theta;
2421  work1 *= grad_q(d); // work1 = dR/dtheta*dtheta/dx
2422  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
2423  }
2424  }
2425  else
2426  {
2427  par_vals.SetDataAndSize(tspec_vals.GetData()+
2428  orientationidx*ndofs, ndofs*3);
2429  par_vals_c1.SetData(par_vals.GetData());
2430  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2431  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2432 
2433  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2434  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2435  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2436  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2437  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2438  grad_e_c1.MultTranspose(shape, grad_q1);
2439  grad_e_c2.MultTranspose(shape, grad_q2);
2440  grad_e_c3.MultTranspose(shape, grad_q3);
2441 
2442  const double theta = shape * par_vals_c1;
2443  const double psi = shape * par_vals_c2;
2444  const double beta = shape * par_vals_c3;
2445 
2446  const double ct = cos(theta), st = sin(theta),
2447  cp = cos(psi), sp = sin(psi),
2448  cb = cos(beta), sb = sin(beta);
2449 
2450  dR_theta = 0.;
2451  dR_theta(0,0) = -st*sp;
2452  dR_theta(1,0) = ct*sp;
2453  dR_theta(2,0) = 0;
2454 
2455  dR_theta(0,1) = -ct*cb - st*cp*sb;
2456  dR_theta(1,1) = -st*cb + ct*cp*sb;
2457  dR_theta(2,1) = 0.;
2458 
2459  dR_theta(0,0) = -ct*sb + st*cp*cb;
2460  dR_theta(1,0) = -st*sb - ct*cp*cb;
2461  dR_theta(2,0) = 0.;
2462 
2463  dR_beta = 0.;
2464  dR_beta(0,0) = 0.;
2465  dR_beta(1,0) = 0.;
2466  dR_beta(2,0) = 0.;
2467 
2468  dR_beta(0,1) = st*sb + ct*cp*cb;
2469  dR_beta(1,1) = -ct*sb + st*cp*cb;
2470  dR_beta(2,1) = -sp*cb;
2471 
2472  dR_beta(0,0) = -st*cb + ct*cp*sb;
2473  dR_beta(1,0) = ct*cb + st*cp*sb;
2474  dR_beta(2,0) = 0.;
2475 
2476  dR_psi = 0.;
2477  dR_psi(0,0) = ct*cp;
2478  dR_psi(1,0) = st*cp;
2479  dR_psi(2,0) = -sp;
2480 
2481  dR_psi(0,1) = 0. - ct*sp*sb;
2482  dR_psi(1,1) = 0. + st*sp*sb;
2483  dR_psi(2,1) = -cp*sb;
2484 
2485  dR_psi(0,0) = 0. + ct*sp*cb;
2486  dR_psi(1,0) = 0. + st*sp*cb;
2487  dR_psi(2,0) = cp*cb;
2488 
2489  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2490  Mult(Jtrcomp_s, work1, work2); // z*Q*D
2491  for (int d = 0; d < dim; d++)
2492  {
2493  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2494  work1 = dR_theta;
2495  work1 *= grad_q1(d); // work1 = dR/dtheta*dtheta/dx
2496  work1.Add(grad_q2(d), dR_psi); // +dR/dpsi*dpsi/dx
2497  work1.Add(grad_q3(d), dR_beta); // +dR/dbeta*dbeta/dx
2498  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
2499  }
2500  }
2501  } // Done orientation
2502  }
2503  break;
2504  }
2505  default:
2506  MFEM_ABORT("Incompatible target type for discrete adaptation!");
2507  }
2508  Jtrcomp.Clear();
2509 }
2510 
2511 void DiscreteAdaptTC:: UpdateGradientTargetSpecification(const Vector &x,
2512  const double dx,
2513  bool use_flag,
2514  int x_ordering)
2515 {
2516  if (use_flag && good_tspec_grad) { return; }
2517 
2518  const int dim = tspec_fesv->GetFE(0)->GetDim(),
2519  cnt = x.Size()/dim;
2520 
2521  tspec_pert1h.SetSize(x.Size()*ncomp);
2522 
2523  Vector TSpecTemp;
2524  Vector xtemp = x;
2525  for (int j = 0; j < dim; j++)
2526  {
2527  for (int i = 0; i < cnt; i++)
2528  {
2529  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2530  xtemp(idx) += dx;
2531  }
2532 
2533  TSpecTemp.NewDataAndSize(tspec_pert1h.GetData() + j*cnt*ncomp, cnt*ncomp);
2534  UpdateTargetSpecification(xtemp, TSpecTemp, x_ordering);
2535 
2536  for (int i = 0; i < cnt; i++)
2537  {
2538  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2539  xtemp(idx) -= dx;
2540  }
2541  }
2542 
2543  good_tspec_grad = use_flag;
2544 }
2545 
2546 void DiscreteAdaptTC::UpdateHessianTargetSpecification(const Vector &x,
2547  double dx, bool use_flag,
2548  int x_ordering)
2549 {
2550 
2551  if (use_flag && good_tspec_hess) { return; }
2552 
2553  const int dim = tspec_fesv->GetFE(0)->GetDim(),
2554  cnt = x.Size()/dim,
2555  totmix = 1+2*(dim-2);
2556 
2557  tspec_pert2h.SetSize(cnt*dim*ncomp);
2558  tspec_pertmix.SetSize(cnt*totmix*ncomp);
2559 
2560  Vector TSpecTemp;
2561  Vector xtemp = x;
2562 
2563  // T(x+2h)
2564  for (int j = 0; j < dim; j++)
2565  {
2566  for (int i = 0; i < cnt; i++)
2567  {
2568  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2569  xtemp(idx) += 2*dx;
2570  }
2571 
2572  TSpecTemp.NewDataAndSize(tspec_pert2h.GetData() + j*cnt*ncomp, cnt*ncomp);
2573  UpdateTargetSpecification(xtemp, TSpecTemp, x_ordering);
2574 
2575  for (int i = 0; i < cnt; i++)
2576  {
2577  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2578  xtemp(idx) -= 2*dx;
2579  }
2580  }
2581 
2582  // T(x+h,y+h)
2583  int j = 0;
2584  for (int k1 = 0; k1 < dim; k1++)
2585  {
2586  for (int k2 = 0; (k1 != k2) && (k2 < dim); k2++)
2587  {
2588  for (int i = 0; i < cnt; i++)
2589  {
2590  int idx1 = x_ordering == Ordering::byNODES ? k1*cnt+i : i*dim + k1;
2591  int idx2 = x_ordering == Ordering::byNODES ? k2*cnt+i : i*dim + k2;
2592  xtemp(idx1) += dx;
2593  xtemp(idx2) += dx;
2594  }
2595 
2596  TSpecTemp.NewDataAndSize(tspec_pertmix.GetData() + j*cnt*ncomp, cnt*ncomp);
2597  UpdateTargetSpecification(xtemp, TSpecTemp, x_ordering);
2598 
2599  for (int i = 0; i < cnt; i++)
2600  {
2601  int idx1 = x_ordering == Ordering::byNODES ? k1*cnt+i : i*dim + k1;
2602  int idx2 = x_ordering == Ordering::byNODES ? k2*cnt+i : i*dim + k2;
2603  xtemp(idx1) -= dx;
2604  xtemp(idx2) -= dx;
2605  }
2606  j++;
2607  }
2608  }
2609 
2610  good_tspec_hess = use_flag;
2611 }
2612 
2613 DiscreteAdaptTC::~DiscreteAdaptTC()
2614 {
2615  delete tspec_gf;
2616  delete adapt_eval;
2617  delete tspec_fesv;
2618 #ifdef MFEM_USE_MPI
2619  delete ptspec_fesv;
2620 #endif
2621 }
2622 
2623 void AdaptivityEvaluator::SetSerialMetaInfo(const Mesh &m,
2624  const FiniteElementSpace &f)
2625 {
2626  delete fes;
2627  delete mesh;
2628  mesh = new Mesh(m, true);
2629  fes = new FiniteElementSpace(mesh, f.FEColl(),
2630  f.GetVDim(), f.GetOrdering());
2631 }
2632 
2633 #ifdef MFEM_USE_MPI
2634 void AdaptivityEvaluator::SetParMetaInfo(const ParMesh &m,
2635  const ParFiniteElementSpace &f)
2636 {
2637  delete pfes;
2638  delete pmesh;
2639  pmesh = new ParMesh(m, true);
2640  pfes = new ParFiniteElementSpace(pmesh, f.FEColl(),
2641  f.GetVDim(), f.GetOrdering());
2642 }
2643 #endif
2644 
2645 void AdaptivityEvaluator::ClearGeometricFactors()
2646 {
2647 #ifdef MFEM_USE_MPI
2648  if (pmesh) { pmesh->DeleteGeometricFactors(); }
2649 #else
2650  if (mesh) { mesh->DeleteGeometricFactors(); }
2651 #endif
2652 }
2653 
2654 AdaptivityEvaluator::~AdaptivityEvaluator()
2655 {
2656  delete fes;
2657  delete mesh;
2658 #ifdef MFEM_USE_MPI
2659  delete pfes;
2660  delete pmesh;
2661 #endif
2662 }
2663 
2664 void TMOP_Integrator::ReleasePADeviceMemory(bool copy_to_host)
2665 {
2666  if (PA.enabled)
2667  {
2668  PA.H.GetMemory().DeleteDevice(copy_to_host);
2669  PA.H0.GetMemory().DeleteDevice(copy_to_host);
2670  if (!copy_to_host && !PA.Jtr.GetMemory().HostIsValid())
2671  {
2672  PA.Jtr_needs_update = true;
2673  }
2674  PA.Jtr.GetMemory().DeleteDevice(copy_to_host);
2675  }
2676 }
2677 
2678 TMOP_Integrator::~TMOP_Integrator()
2679 {
2680  delete lim_func;
2681  delete adapt_lim_gf;
2682  delete surf_fit_gf;
2683  for (int i = 0; i < ElemDer.Size(); i++)
2684  {
2685  delete ElemDer[i];
2686  delete ElemPertEnergy[i];
2687  }
2688 }
2689 
2690 void TMOP_Integrator::EnableLimiting(const GridFunction &n0,
2691  const GridFunction &dist, Coefficient &w0,
2692  TMOP_LimiterFunction *lfunc)
2693 {
2694  lim_nodes0 = &n0;
2695  lim_coeff = &w0;
2696  lim_dist = &dist;
2697  MFEM_VERIFY(lim_dist->FESpace()->GetVDim() == 1,
2698  "'dist' must be a scalar GridFunction!");
2699 
2700  delete lim_func;
2701  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2702 }
2703 
2704 void TMOP_Integrator::EnableLimiting(const GridFunction &n0, Coefficient &w0,
2705  TMOP_LimiterFunction *lfunc)
2706 {
2707  lim_nodes0 = &n0;
2708  lim_coeff = &w0;
2709  lim_dist = NULL;
2710 
2711  delete lim_func;
2712  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2713 }
2714 
2715 void TMOP_Integrator::EnableAdaptiveLimiting(const GridFunction &z0,
2716  Coefficient &coeff,
2717  AdaptivityEvaluator &ae)
2718 {
2719  adapt_lim_gf0 = &z0;
2720  delete adapt_lim_gf;
2721  adapt_lim_gf = new GridFunction(z0);
2722  adapt_lim_coeff = &coeff;
2723  adapt_lim_eval = &ae;
2724 
2725  adapt_lim_eval->SetSerialMetaInfo(*z0.FESpace()->GetMesh(),
2726  *z0.FESpace());
2727  adapt_lim_eval->SetInitialField
2728  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2729 }
2730 
2731 #ifdef MFEM_USE_MPI
2732 void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
2733  Coefficient &coeff,
2734  AdaptivityEvaluator &ae)
2735 {
2736  adapt_lim_gf0 = &z0;
2737  adapt_lim_pgf0 = &z0;
2738  delete adapt_lim_gf;
2739  adapt_lim_gf = new GridFunction(z0);
2740  adapt_lim_coeff = &coeff;
2741  adapt_lim_eval = &ae;
2742 
2743  adapt_lim_eval->SetParMetaInfo(*z0.ParFESpace()->GetParMesh(),
2744  *z0.ParFESpace());
2745  adapt_lim_eval->SetInitialField
2746  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2747 }
2748 #endif
2749 
2750 void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &s0,
2751  const Array<bool> &smarker,
2752  Coefficient &coeff,
2753  AdaptivityEvaluator &ae)
2754 {
2755  delete surf_fit_gf;
2756  surf_fit_gf = new GridFunction(s0);
2757  surf_fit_marker = &smarker;
2758  surf_fit_coeff = &coeff;
2759  surf_fit_eval = &ae;
2760 
2761  surf_fit_eval->SetSerialMetaInfo(*s0.FESpace()->GetMesh(),
2762  *s0.FESpace());
2763  surf_fit_eval->SetInitialField
2764  (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
2765 }
2766 
2767 #ifdef MFEM_USE_MPI
2768 void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
2769  const Array<bool> &smarker,
2770  Coefficient &coeff,
2771  AdaptivityEvaluator &ae)
2772 {
2773  delete surf_fit_gf;
2774  surf_fit_gf = new GridFunction(s0);
2775  surf_fit_marker = &smarker;
2776  surf_fit_coeff = &coeff;
2777  surf_fit_eval = &ae;
2778 
2779  surf_fit_eval->SetParMetaInfo(*s0.ParFESpace()->GetParMesh(),
2780  *s0.ParFESpace());
2781  surf_fit_eval->SetInitialField
2782  (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
2783 }
2784 #endif
2785 
2786 void TMOP_Integrator::GetSurfaceFittingErrors(double &err_avg, double &err_max)
2787 {
2788  MFEM_VERIFY(surf_fit_gf, "Surface fitting has not been enabled.");
2789 
2790  int loc_cnt = 0;
2791  double loc_max = 0.0, loc_sum = 0.0;
2792  for (int i = 0; i < surf_fit_marker->Size(); i++)
2793  {
2794  if ((*surf_fit_marker)[i] == true)
2795  {
2796  loc_cnt++;
2797  loc_max = std::max(loc_max, std::abs((*surf_fit_gf)(i)));
2798  loc_sum += std::abs((*surf_fit_gf)(i));
2799  }
2800  }
2801  err_avg = loc_sum / loc_cnt;
2802  err_max = loc_max;
2803 
2804 #ifdef MFEM_USE_MPI
2805  if (targetC->Parallel() == false) { return; }
2806  int glob_cnt;
2807  MPI_Comm comm = targetC->GetComm();
2808  MPI_Allreduce(&loc_max, &err_max, 1, MPI_DOUBLE, MPI_MAX, comm);
2809  MPI_Allreduce(&loc_cnt, &glob_cnt, 1, MPI_INT, MPI_SUM, comm);
2810  MPI_Allreduce(&loc_sum, &err_avg, 1, MPI_DOUBLE, MPI_SUM, comm);
2811  err_avg = err_avg / glob_cnt;
2812 #endif
2813 }
2814 
2815 void TMOP_Integrator::UpdateAfterMeshTopologyChange()
2816 {
2817  if (adapt_lim_gf)
2818  {
2819  adapt_lim_gf->Update();
2820  adapt_lim_eval->SetSerialMetaInfo(*adapt_lim_gf->FESpace()->GetMesh(),
2821  *adapt_lim_gf->FESpace());
2822  adapt_lim_eval->SetInitialField
2823  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2824  }
2825 }
2826 
2827 #ifdef MFEM_USE_MPI
2828 void TMOP_Integrator::ParUpdateAfterMeshTopologyChange()
2829 {
2830  if (adapt_lim_gf)
2831  {
2832  adapt_lim_gf->Update();
2833  adapt_lim_eval->SetParMetaInfo(*adapt_lim_pgf0->ParFESpace()->GetParMesh(),
2834  *adapt_lim_pgf0->ParFESpace());
2835  adapt_lim_eval->SetInitialField
2836  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2837  }
2838 }
2839 #endif
2840 
2841 double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
2842  ElementTransformation &T,
2843  const Vector &elfun)
2844 {
2845  const int dof = el.GetDof(), dim = el.GetDim();
2846  const int el_id = T.ElementNo;
2847  double energy;
2848 
2849  // No adaptive limiting / surface fitting terms if the function is called
2850  // as part of a FD derivative computation (because we include the exact
2851  // derivatives of these terms in FD computations).
2852  const bool adaptive_limiting = (adapt_lim_gf && fd_call_flag == false);
2853  const bool surface_fit = (surf_fit_gf && fd_call_flag == false);
2854 
2855  DSh.SetSize(dof, dim);
2856  Jrt.SetSize(dim);
2857  Jpr.SetSize(dim);
2858  Jpt.SetSize(dim);
2859  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2860 
2861  const IntegrationRule &ir = EnergyIntegrationRule(el);
2862 
2863  energy = 0.0;
2864  DenseTensor Jtr(dim, dim, ir.GetNPoints());
2865  targetC->ComputeElementTargets(el_id, el, ir, elfun, Jtr);
2866 
2867  // Limited case.
2868  Vector shape, p, p0, d_vals;
2869  DenseMatrix pos0;
2870  if (lim_coeff)
2871  {
2872  shape.SetSize(dof);
2873  p.SetSize(dim);
2874  p0.SetSize(dim);
2875  pos0.SetSize(dof, dim);
2876  Vector pos0V(pos0.Data(), dof * dim);
2877  Array<int> pos_dofs;
2878  lim_nodes0->FESpace()->GetElementVDofs(el_id, pos_dofs);
2879  lim_nodes0->GetSubVector(pos_dofs, pos0V);
2880  if (lim_dist)
2881  {
2882  lim_dist->GetValues(el_id, ir, d_vals);
2883  }
2884  else
2885  {
2886  d_vals.SetSize(ir.GetNPoints()); d_vals = 1.0;
2887  }
2888  }
2889 
2890  // Define ref->physical transformation, when a Coefficient is specified.
2891  IsoparametricTransformation *Tpr = NULL;
2892  if (metric_coeff || lim_coeff || adaptive_limiting || surface_fit)
2893  {
2894  Tpr = new IsoparametricTransformation;
2895  Tpr->SetFE(&el);
2896  Tpr->ElementNo = el_id;
2897  Tpr->ElementType = ElementTransformation::ELEMENT;
2898  Tpr->Attribute = T.Attribute;
2899  Tpr->mesh = T.mesh;
2900  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2901  }
2902  // TODO: computing the coefficients 'metric_coeff' and 'lim_coeff' in physical
2903  // coordinates means that, generally, the gradient and Hessian of the
2904  // TMOP_Integrator will depend on the derivatives of the coefficients.
2905  //
2906  // In some cases the coefficients are independent of any movement of
2907  // the physical coordinates (i.e. changes in 'elfun'), e.g. when the
2908  // coefficient is a ConstantCoefficient or a GridFunctionCoefficient.
2909 
2910  Vector adapt_lim_gf_q, adapt_lim_gf0_q;
2911  if (adaptive_limiting)
2912  {
2913  adapt_lim_gf->GetValues(el_id, ir, adapt_lim_gf_q);
2914  adapt_lim_gf0->GetValues(el_id, ir, adapt_lim_gf0_q);
2915  }
2916 
2917  for (int i = 0; i < ir.GetNPoints(); i++)
2918  {
2919  const IntegrationPoint &ip = ir.IntPoint(i);
2920 
2921  const DenseMatrix &Jtr_i = Jtr(i);
2922  metric->SetTargetJacobian(Jtr_i);
2923  CalcInverse(Jtr_i, Jrt);
2924  const double weight = ip.weight * Jtr_i.Det();
2925 
2926  el.CalcDShape(ip, DSh);
2927  MultAtB(PMatI, DSh, Jpr);
2928  Mult(Jpr, Jrt, Jpt);
2929 
2930  double val = metric_normal * metric->EvalW(Jpt);
2931  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
2932 
2933  if (lim_coeff)
2934  {
2935  el.CalcShape(ip, shape);
2936  PMatI.MultTranspose(shape, p);
2937  pos0.MultTranspose(shape, p0);
2938  val += lim_normal *
2939  lim_func->Eval(p, p0, d_vals(i)) *
2940  lim_coeff->Eval(*Tpr, ip);
2941  }
2942 
2943  // Contribution from the adaptive limiting term.
2944  if (adaptive_limiting)
2945  {
2946  const double diff = adapt_lim_gf_q(i) - adapt_lim_gf0_q(i);
2947  val += adapt_lim_coeff->Eval(*Tpr, ip) * lim_normal * diff * diff;
2948  }
2949 
2950  energy += weight * val;
2951  }
2952 
2953  // Contribution from the surface fitting term.
2954  if (surface_fit)
2955  {
2956  const IntegrationRule &ir_s =
2957  surf_fit_gf->FESpace()->GetFE(el_id)->GetNodes();
2958  Array<int> dofs;
2959  Vector sigma_e;
2960  surf_fit_gf->FESpace()->GetElementDofs(el_id, dofs);
2961  surf_fit_gf->GetSubVector(dofs, sigma_e);
2962  for (int s = 0; s < dofs.Size(); s++)
2963  {
2964  if ((*surf_fit_marker)[dofs[s]] == true)
2965  {
2966  const IntegrationPoint &ip_s = ir_s.IntPoint(s);
2967  Tpr->SetIntPoint(&ip_s);
2968  energy += surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
2969  sigma_e(s) * sigma_e(s);
2970  }
2971  }
2972  }
2973 
2974  delete Tpr;
2975  return energy;
2976 }
2977 
2978 double TMOP_Integrator::GetRefinementElementEnergy(const FiniteElement &el,
2979  ElementTransformation &T,
2980  const Vector &elfun,
2981  const IntegrationRule &irule)
2982 {
2983  int dof = el.GetDof(), dim = el.GetDim(),
2984  NEsplit = elfun.Size() / (dof*dim), el_id = T.ElementNo;
2985  double energy = 0.;
2986 
2987  TargetConstructor *tc = const_cast<TargetConstructor *>(targetC);
2988  DiscreteAdaptTC *dtc = dynamic_cast<DiscreteAdaptTC *>(tc);
2989  // For DiscreteAdaptTC the GridFunctions used to set the targets must be
2990  // mapped onto the fine elements.
2991  if (dtc) { dtc->SetTspecFromIntRule(el_id, irule); }
2992 
2993  for (int e = 0; e < NEsplit; e++)
2994  {
2995  DSh.SetSize(dof, dim);
2996  Jrt.SetSize(dim);
2997  Jpr.SetSize(dim);
2998  Jpt.SetSize(dim);
2999  Vector elfun_child(dof*dim);
3000  for (int i = 0; i < dof; i++)
3001  {
3002  for (int d = 0; d < dim; d++)
3003  {
3004  // elfun is (xe1,xe2,...xen,ye1,ye2...yen) and has nodal coordinates
3005  // for all the children element of the parent element being considered.
3006  // So we must index and get (xek, yek) i.e. nodal coordinates for
3007  // the fine element being considered.
3008  elfun_child(i + d*dof) = elfun(i + e*dof + d*dof*NEsplit);
3009  }
3010  }
3011  PMatI.UseExternalData(elfun_child.GetData(), dof, dim);
3012 
3013  const IntegrationRule &ir = EnergyIntegrationRule(el);
3014 
3015  double el_energy = 0;
3016  DenseTensor Jtr(dim, dim, ir.GetNPoints());
3017  if (dtc)
3018  {
3019  // This is used to index into the tspec vector inside DiscreteAdaptTC.
3020  dtc->SetRefinementSubElement(e);
3021  }
3022  targetC->ComputeElementTargets(el_id, el, ir, elfun_child, Jtr);
3023 
3024  // Define ref->physical transformation, wn a Coefficient is specified.
3025  IsoparametricTransformation *Tpr = NULL;
3026  if (metric_coeff || lim_coeff)
3027  {
3028  Tpr = new IsoparametricTransformation;
3029  Tpr->SetFE(&el);
3030  Tpr->ElementNo = T.ElementNo;
3031  Tpr->ElementType = ElementTransformation::ELEMENT;
3032  Tpr->Attribute = T.Attribute;
3033  Tpr->mesh = T.mesh;
3034  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3035  }
3036 
3037  for (int i = 0; i < ir.GetNPoints(); i++)
3038  {
3039  const IntegrationPoint &ip = ir.IntPoint(i);
3040  const DenseMatrix &Jtr_i = Jtr(i);
3041  h_metric->SetTargetJacobian(Jtr_i);
3042  CalcInverse(Jtr_i, Jrt);
3043  const double weight = ip.weight * Jtr_i.Det();
3044 
3045  el.CalcDShape(ip, DSh);
3046  MultAtB(PMatI, DSh, Jpr);
3047  Mult(Jpr, Jrt, Jpt);
3048 
3049  double val = metric_normal * h_metric->EvalW(Jpt);
3050  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
3051 
3052  el_energy += weight * val;
3053  delete Tpr;
3054  }
3055  energy += el_energy;
3056  }
3057  energy /= NEsplit;
3058 
3059  if (dtc) { dtc->ResetRefinementTspecData(); }
3060 
3061  return energy;
3062 }
3063 
3064 double TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
3065  ElementTransformation &T,
3066  const Vector &elfun)
3067 {
3068  int dof = el.GetDof(), dim = el.GetDim();
3069  double energy = 0.;
3070 
3071  DSh.SetSize(dof, dim);
3072  Jrt.SetSize(dim);
3073  Jpr.SetSize(dim);
3074  Jpt.SetSize(dim);
3075  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3076 
3077  const IntegrationRule &ir = EnergyIntegrationRule(el);
3078 
3079  energy = 0.0;
3080  DenseTensor Jtr(dim, dim, ir.GetNPoints());
3081  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3082 
3083  // Define ref->physical transformation, wn a Coefficient is specified.
3084  IsoparametricTransformation *Tpr = NULL;
3085  if (metric_coeff)
3086  {
3087  Tpr = new IsoparametricTransformation;
3088  Tpr->SetFE(&el);
3089  Tpr->ElementNo = T.ElementNo;
3090  Tpr->ElementType = ElementTransformation::ELEMENT;
3091  Tpr->Attribute = T.Attribute;
3092  Tpr->mesh = T.mesh;
3093  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3094  }
3095 
3096  for (int i = 0; i < ir.GetNPoints(); i++)
3097  {
3098  const IntegrationPoint &ip = ir.IntPoint(i);
3099  const DenseMatrix &Jtr_i = Jtr(i);
3100  h_metric->SetTargetJacobian(Jtr_i);
3101  CalcInverse(Jtr_i, Jrt);
3102  const double weight = ip.weight * Jtr_i.Det();
3103 
3104  el.CalcDShape(ip, DSh);
3105  MultAtB(PMatI, DSh, Jpr);
3106  Mult(Jpr, Jrt, Jpt);
3107 
3108  double val = metric_normal * h_metric->EvalW(Jpt);
3109  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
3110 
3111  energy += weight * val;
3112  }
3113 
3114  delete Tpr;
3115  return energy;
3116 }
3117 
3118 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
3119  ElementTransformation &T,
3120  const Vector &elfun, Vector &elvect)
3121 {
3122  if (!fdflag)
3123  {
3124  AssembleElementVectorExact(el, T, elfun, elvect);
3125  }
3126  else
3127  {
3128  AssembleElementVectorFD(el, T, elfun, elvect);
3129  }
3130 }
3131 
3132 void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,
3133  ElementTransformation &T,
3134  const Vector &elfun,
3135  DenseMatrix &elmat)
3136 {
3137  if (!fdflag)
3138  {
3139  AssembleElementGradExact(el, T, elfun, elmat);
3140  }
3141  else
3142  {
3143  AssembleElementGradFD(el, T, elfun, elmat);
3144  }
3145 }
3146 
3147 void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
3148  ElementTransformation &T,
3149  const Vector &elfun,
3150  Vector &elvect)
3151 {
3152  const int dof = el.GetDof(), dim = el.GetDim();
3153 
3154  DenseMatrix Amat(dim), work1(dim), work2(dim);
3155  DSh.SetSize(dof, dim);
3156  DS.SetSize(dof, dim);
3157  Jrt.SetSize(dim);
3158  Jpt.SetSize(dim);
3159  P.SetSize(dim);
3160  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3161  elvect.SetSize(dof*dim);
3162  PMatO.UseExternalData(elvect.GetData(), dof, dim);
3163 
3164  const IntegrationRule &ir = ActionIntegrationRule(el);
3165  const int nqp = ir.GetNPoints();
3166 
3167  elvect = 0.0;
3168  Vector weights(nqp);
3169  DenseTensor Jtr(dim, dim, nqp);
3170  DenseTensor dJtr(dim, dim, dim*nqp);
3171  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3172 
3173  // Limited case.
3174  DenseMatrix pos0;
3175  Vector shape, p, p0, d_vals, grad;
3176  shape.SetSize(dof);
3177  if (lim_coeff)
3178  {
3179  p.SetSize(dim);
3180  p0.SetSize(dim);
3181  pos0.SetSize(dof, dim);
3182  Vector pos0V(pos0.Data(), dof * dim);
3183  Array<int> pos_dofs;
3184  lim_nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
3185  lim_nodes0->GetSubVector(pos_dofs, pos0V);
3186  if (lim_dist)
3187  {
3188  lim_dist->GetValues(T.ElementNo, ir, d_vals);
3189  }
3190  else
3191  {
3192  d_vals.SetSize(nqp); d_vals = 1.0;
3193  }
3194  }
3195 
3196  // Define ref->physical transformation, when a Coefficient is specified.
3197  IsoparametricTransformation *Tpr = NULL;
3198  if (metric_coeff || lim_coeff || adapt_lim_gf || surf_fit_gf || exact_action)
3199  {
3200  Tpr = new IsoparametricTransformation;
3201  Tpr->SetFE(&el);
3202  Tpr->ElementNo = T.ElementNo;
3203  Tpr->ElementType = ElementTransformation::ELEMENT;
3204  Tpr->Attribute = T.Attribute;
3205  Tpr->mesh = T.mesh;
3206  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3207  if (exact_action)
3208  {
3209  targetC->ComputeElementTargetsGradient(ir, elfun, *Tpr, dJtr);
3210  }
3211  }
3212 
3213  Vector d_detW_dx(dim);
3214  Vector d_Winv_dx(dim);
3215 
3216  for (int q = 0; q < nqp; q++)
3217  {
3218  const IntegrationPoint &ip = ir.IntPoint(q);
3219  const DenseMatrix &Jtr_q = Jtr(q);
3220  metric->SetTargetJacobian(Jtr_q);
3221  CalcInverse(Jtr_q, Jrt);
3222  weights(q) = ip.weight * Jtr_q.Det();
3223  double weight_m = weights(q) * metric_normal;
3224 
3225  el.CalcDShape(ip, DSh);
3226  Mult(DSh, Jrt, DS);
3227  MultAtB(PMatI, DS, Jpt);
3228 
3229  metric->EvalP(Jpt, P);
3230 
3231  if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
3232 
3233  P *= weight_m;
3234  AddMultABt(DS, P, PMatO); // w_q det(W) dmu/dx : dA/dx Winv
3235 
3236  if (exact_action)
3237  {
3238  el.CalcShape(ip, shape);
3239  // Derivatives of adaptivity-based targets.
3240  // First term: w_q d*(Det W)/dx * mu(T)
3241  // d(Det W)/dx = det(W)*Tr[Winv*dW/dx]
3242  DenseMatrix dwdx(dim);
3243  for (int d = 0; d < dim; d++)
3244  {
3245  const DenseMatrix &dJtr_q = dJtr(q + d * nqp);
3246  Mult(Jrt, dJtr_q, dwdx );
3247  d_detW_dx(d) = dwdx.Trace();
3248  }
3249  d_detW_dx *= weight_m*metric->EvalW(Jpt); // *[w_q*det(W)]*mu(T)
3250 
3251  // Second term: w_q det(W) dmu/dx : AdWinv/dx
3252  // dWinv/dx = -Winv*dW/dx*Winv
3253  MultAtB(PMatI, DSh, Amat);
3254  for (int d = 0; d < dim; d++)
3255  {
3256  const DenseMatrix &dJtr_q = dJtr(q + d*nqp);
3257  Mult(Jrt, dJtr_q, work1); // Winv*dw/dx
3258  Mult(work1, Jrt, work2); // Winv*dw/dx*Winv
3259  Mult(Amat, work2, work1); // A*Winv*dw/dx*Winv
3260  MultAtB(P, work1, work2); // dmu/dT^T*A*Winv*dw/dx*Winv
3261  d_Winv_dx(d) = work2.Trace(); // Tr[dmu/dT : AWinv*dw/dx*Winv]
3262  }
3263  d_Winv_dx *= -weight_m; // Include (-) factor as well
3264 
3265  d_detW_dx += d_Winv_dx;
3266  AddMultVWt(shape, d_detW_dx, PMatO);
3267  }
3268 
3269  if (lim_coeff)
3270  {
3271  if (!exact_action) { el.CalcShape(ip, shape); }
3272  PMatI.MultTranspose(shape, p);
3273  pos0.MultTranspose(shape, p0);
3274  lim_func->Eval_d1(p, p0, d_vals(q), grad);
3275  grad *= weights(q) * lim_normal * lim_coeff->Eval(*Tpr, ip);
3276  AddMultVWt(shape, grad, PMatO);
3277  }
3278  }
3279 
3280  if (adapt_lim_gf) { AssembleElemVecAdaptLim(el, *Tpr, ir, weights, PMatO); }
3281  if (surf_fit_gf) { AssembleElemVecSurfFit(el, *Tpr, PMatO); }
3282 
3283  delete Tpr;
3284 }
3285 
3286 void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
3287  ElementTransformation &T,
3288  const Vector &elfun,
3289  DenseMatrix &elmat)
3290 {
3291  const int dof = el.GetDof(), dim = el.GetDim();
3292 
3293  DSh.SetSize(dof, dim);
3294  DS.SetSize(dof, dim);
3295  Jrt.SetSize(dim);
3296  Jpt.SetSize(dim);
3297  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3298  elmat.SetSize(dof*dim);
3299 
3300  const IntegrationRule &ir = GradientIntegrationRule(el);
3301  const int nqp = ir.GetNPoints();
3302 
3303  elmat = 0.0;
3304  Vector weights(nqp);
3305  DenseTensor Jtr(dim, dim, nqp);
3306  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3307 
3308  // Limited case.
3309  DenseMatrix pos0, hess;
3310  Vector shape, p, p0, d_vals;
3311  if (lim_coeff)
3312  {
3313  shape.SetSize(dof);
3314  p.SetSize(dim);
3315  p0.SetSize(dim);
3316  pos0.SetSize(dof, dim);
3317  Vector pos0V(pos0.Data(), dof * dim);
3318  Array<int> pos_dofs;
3319  lim_nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
3320  lim_nodes0->GetSubVector(pos_dofs, pos0V);
3321  if (lim_dist)
3322  {
3323  lim_dist->GetValues(T.ElementNo, ir, d_vals);
3324  }
3325  else
3326  {
3327  d_vals.SetSize(nqp); d_vals = 1.0;
3328  }
3329  }
3330 
3331  // Define ref->physical transformation, when a Coefficient is specified.
3332  IsoparametricTransformation *Tpr = NULL;
3333  if (metric_coeff || lim_coeff || adapt_lim_gf || surf_fit_gf)
3334  {
3335  Tpr = new IsoparametricTransformation;
3336  Tpr->SetFE(&el);
3337  Tpr->ElementNo = T.ElementNo;
3338  Tpr->ElementType = ElementTransformation::ELEMENT;
3339  Tpr->Attribute = T.Attribute;
3340  Tpr->mesh = T.mesh;
3341  Tpr->GetPointMat().Transpose(PMatI);
3342  }
3343 
3344  for (int q = 0; q < nqp; q++)
3345  {
3346  const IntegrationPoint &ip = ir.IntPoint(q);
3347  const DenseMatrix &Jtr_q = Jtr(q);
3348  metric->SetTargetJacobian(Jtr_q);
3349  CalcInverse(Jtr_q, Jrt);
3350  weights(q) = ip.weight * Jtr_q.Det();
3351  double weight_m = weights(q) * metric_normal;
3352 
3353  el.CalcDShape(ip, DSh);
3354  Mult(DSh, Jrt, DS);
3355  MultAtB(PMatI, DS, Jpt);
3356 
3357  if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
3358 
3359  metric->AssembleH(Jpt, DS, weight_m, elmat);
3360 
3361  // TODO: derivatives of adaptivity-based targets.
3362 
3363  if (lim_coeff)
3364  {
3365  el.CalcShape(ip, shape);
3366  PMatI.MultTranspose(shape, p);
3367  pos0.MultTranspose(shape, p0);
3368  weight_m = weights(q) * lim_normal * lim_coeff->Eval(*Tpr, ip);
3369  lim_func->Eval_d2(p, p0, d_vals(q), hess);
3370  for (int i = 0; i < dof; i++)
3371  {
3372  const double w_shape_i = weight_m * shape(i);
3373  for (int j = 0; j < dof; j++)
3374  {
3375  const double w = w_shape_i * shape(j);
3376  for (int d1 = 0; d1 < dim; d1++)
3377  {
3378  for (int d2 = 0; d2 < dim; d2++)
3379  {
3380  elmat(d1*dof + i, d2*dof + j) += w * hess(d1, d2);
3381  }
3382  }
3383  }
3384  }
3385  }
3386  }
3387 
3388  if (adapt_lim_gf) { AssembleElemGradAdaptLim(el, *Tpr, ir, weights, elmat); }
3389  if (surf_fit_gf) { AssembleElemGradSurfFit(el, *Tpr, elmat); }
3390 
3391  delete Tpr;
3392 }
3393 
3394 void TMOP_Integrator::AssembleElemVecAdaptLim(const FiniteElement &el,
3395  IsoparametricTransformation &Tpr,
3396  const IntegrationRule &ir,
3397  const Vector &weights,
3398  DenseMatrix &mat)
3399 {
3400  const int dof = el.GetDof(), dim = el.GetDim(), nqp = weights.Size();
3401  Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3402 
3403  Array<int> dofs;
3404  adapt_lim_gf->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
3405  adapt_lim_gf->GetSubVector(dofs, adapt_lim_gf_e);
3406  adapt_lim_gf->GetValues(Tpr.ElementNo, ir, adapt_lim_gf_q);
3407  adapt_lim_gf0->GetValues(Tpr.ElementNo, ir, adapt_lim_gf0_q);
3408 
3409  // Project the gradient of adapt_lim_gf in the same space.
3410  // The FE coefficients of the gradient go in adapt_lim_gf_grad_e.
3411  DenseMatrix adapt_lim_gf_grad_e(dof, dim);
3412  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3413  el.ProjectGrad(el, Tpr, grad_phys);
3414  Vector grad_ptr(adapt_lim_gf_grad_e.GetData(), dof*dim);
3415  grad_phys.Mult(adapt_lim_gf_e, grad_ptr);
3416 
3417  Vector adapt_lim_gf_grad_q(dim);
3418 
3419  for (int q = 0; q < nqp; q++)
3420  {
3421  const IntegrationPoint &ip = ir.IntPoint(q);
3422  el.CalcShape(ip, shape);
3423  adapt_lim_gf_grad_e.MultTranspose(shape, adapt_lim_gf_grad_q);
3424  adapt_lim_gf_grad_q *= 2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q));
3425  adapt_lim_gf_grad_q *= weights(q) * lim_normal * adapt_lim_coeff->Eval(Tpr, ip);
3426  AddMultVWt(shape, adapt_lim_gf_grad_q, mat);
3427  }
3428 }
3429 
3430 void TMOP_Integrator::AssembleElemGradAdaptLim(const FiniteElement &el,
3431  IsoparametricTransformation &Tpr,
3432  const IntegrationRule &ir,
3433  const Vector &weights,
3434  DenseMatrix &mat)
3435 {
3436  const int dof = el.GetDof(), dim = el.GetDim(), nqp = weights.Size();
3437  Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3438 
3439  Array<int> dofs;
3440  adapt_lim_gf->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
3441  adapt_lim_gf->GetSubVector(dofs, adapt_lim_gf_e);
3442  adapt_lim_gf->GetValues(Tpr.ElementNo, ir, adapt_lim_gf_q);
3443  adapt_lim_gf0->GetValues(Tpr.ElementNo, ir, adapt_lim_gf0_q);
3444 
3445  // Project the gradient of adapt_lim_gf in the same space.
3446  // The FE coefficients of the gradient go in adapt_lim_gf_grad_e.
3447  DenseMatrix adapt_lim_gf_grad_e(dof, dim);
3448  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3449  el.ProjectGrad(el, Tpr, grad_phys);
3450  Vector grad_ptr(adapt_lim_gf_grad_e.GetData(), dof*dim);
3451  grad_phys.Mult(adapt_lim_gf_e, grad_ptr);
3452 
3453  // Project the gradient of each gradient of adapt_lim_gf in the same space.
3454  // The FE coefficients of the second derivatives go in adapt_lim_gf_hess_e.
3455  DenseMatrix adapt_lim_gf_hess_e(dof*dim, dim);
3456  Mult(grad_phys, adapt_lim_gf_grad_e, adapt_lim_gf_hess_e);
3457  // Reshape to be more convenient later (no change in the data).
3458  adapt_lim_gf_hess_e.SetSize(dof, dim*dim);
3459 
3460  Vector adapt_lim_gf_grad_q(dim);
3461  DenseMatrix adapt_lim_gf_hess_q(dim, dim);
3462 
3463  for (int q = 0; q < nqp; q++)
3464  {
3465  const IntegrationPoint &ip = ir.IntPoint(q);
3466  el.CalcShape(ip, shape);
3467 
3468  adapt_lim_gf_grad_e.MultTranspose(shape, adapt_lim_gf_grad_q);
3469  Vector gg_ptr(adapt_lim_gf_hess_q.GetData(), dim*dim);
3470  adapt_lim_gf_hess_e.MultTranspose(shape, gg_ptr);
3471 
3472  const double w = weights(q) * lim_normal * adapt_lim_coeff->Eval(Tpr, ip);
3473  for (int i = 0; i < dof * dim; i++)
3474  {
3475  const int idof = i % dof, idim = i / dof;
3476  for (int j = 0; j <= i; j++)
3477  {
3478  const int jdof = j % dof, jdim = j / dof;
3479  const double entry =
3480  w * ( 2.0 * adapt_lim_gf_grad_q(idim) * shape(idof) *
3481  /* */ adapt_lim_gf_grad_q(jdim) * shape(jdof) +
3482  2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q)) *
3483  adapt_lim_gf_hess_q(idim, jdim) * shape(idof) * shape(jdof));
3484  mat(i, j) += entry;
3485  if (i != j) { mat(j, i) += entry; }
3486  }
3487  }
3488  }
3489 }
3490 
3491 void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x,
3492  IsoparametricTransformation &Tpr,
3493  DenseMatrix &mat)
3494 {
3495  const int el_id = Tpr.ElementNo;
3496  // Check if the element has any DOFs marked for surface fitting.
3497  Array<int> dofs;
3498  surf_fit_gf->FESpace()->GetElementDofs(el_id, dofs);
3499  int count = 0;
3500  for (int s = 0; s < dofs.Size(); s++)
3501  {
3502  count += ((*surf_fit_marker)[dofs[s]]) ? 1 : 0;
3503  }
3504  if (count == 0) { return; }
3505 
3506  const FiniteElement &el_s = *surf_fit_gf->FESpace()->GetFE(el_id);
3507 
3508  const int dof_x = el_x.GetDof(), dim = el_x.GetDim(),
3509  dof_s = el_s.GetDof();
3510 
3511  Vector sigma_e;
3512  surf_fit_gf->GetSubVector(dofs, sigma_e);
3513 
3514  // Project the gradient of sigma in the same space.
3515  // The FE coefficients of the gradient go in surf_fit_grad_e.
3516  DenseMatrix surf_fit_grad_e(dof_s, dim);
3517  Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
3518  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3519  el_s.ProjectGrad(el_s, Tpr, grad_phys);
3520  grad_phys.Mult(sigma_e, grad_ptr);
3521 
3522  Vector shape_x(dof_x), shape_s(dof_s);
3523  const IntegrationRule &ir = el_s.GetNodes();
3524  Vector surf_fit_grad_s(dim);
3525  surf_fit_grad_s = 0.0;
3526 
3527  for (int s = 0; s < dof_s; s++)
3528  {
3529  if ((*surf_fit_marker)[dofs[s]] == false) { continue; }
3530 
3531  const IntegrationPoint &ip = ir.IntPoint(s);
3532  Tpr.SetIntPoint(&ip);
3533  el_x.CalcShape(ip, shape_x);
3534  el_s.CalcShape(ip, shape_s);
3535 
3536  // Note that this gradient is already in physical space.
3537  surf_fit_grad_e.MultTranspose(shape_s, surf_fit_grad_s);
3538  surf_fit_grad_s *= 2.0 * surf_fit_normal *
3539  surf_fit_coeff->Eval(Tpr, ip) * sigma_e(s);
3540 
3541  AddMultVWt(shape_x, surf_fit_grad_s, mat);
3542  }
3543 }
3544 
3545 void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
3546  IsoparametricTransformation &Tpr,
3547  DenseMatrix &mat)
3548 {
3549  const int el_id = Tpr.ElementNo;
3550  // Check if the element has any DOFs marked for surface fitting.
3551  Array<int> dofs;
3552  surf_fit_gf->FESpace()->GetElementDofs(el_id, dofs);
3553  int count = 0;
3554  for (int s = 0; s < dofs.Size(); s++)
3555  {
3556  count += ((*surf_fit_marker)[dofs[s]]) ? 1 : 0;
3557  }
3558  if (count == 0) { return; }
3559 
3560  const FiniteElement &el_s = *surf_fit_gf->FESpace()->GetFE(el_id);
3561 
3562  const int dof_x = el_x.GetDof(), dim = el_x.GetDim(),
3563  dof_s = el_s.GetDof();
3564 
3565  Vector sigma_e;
3566  surf_fit_gf->GetSubVector(dofs, sigma_e);
3567 
3568  DenseMatrix surf_fit_grad_e(dof_s, dim);
3569  Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
3570  DenseMatrix grad_phys;
3571  el_s.ProjectGrad(el_s, Tpr, grad_phys);
3572  grad_phys.Mult(sigma_e, grad_ptr);
3573 
3574  DenseMatrix surf_fit_hess_e(dof_s, dim*dim);
3575  Vector hess_ptr(surf_fit_hess_e.GetData(), dof_s*dim*dim);
3576  surf_fit_hess_e.SetSize(dof_s*dim, dim);
3577  Mult(grad_phys, surf_fit_grad_e, surf_fit_hess_e);
3578  surf_fit_hess_e.SetSize(dof_s, dim * dim);
3579 
3580  const IntegrationRule &ir = el_s.GetNodes();
3581  Vector shape_x(dof_x), shape_s(dof_s);
3582 
3583  Vector surf_fit_grad_s(dim);
3584  DenseMatrix surf_fit_hess_s(dim, dim);
3585 
3586 
3587  for (int s = 0; s < dof_s; s++)
3588  {
3589  if ((*surf_fit_marker)[dofs[s]] == false) { continue; }
3590 
3591  const IntegrationPoint &ip = ir.IntPoint(s);
3592  Tpr.SetIntPoint(&ip);
3593  el_x.CalcShape(ip, shape_x);
3594  el_s.CalcShape(ip, shape_s);
3595 
3596  // These are the sums over k at the dof s (looking at the notes).
3597  surf_fit_grad_e.MultTranspose(shape_s, surf_fit_grad_s);
3598  Vector gg_ptr(surf_fit_hess_s.GetData(), dim * dim);
3599  surf_fit_hess_e.MultTranspose(shape_s, gg_ptr);
3600 
3601  // Loops over the local matrix.
3602  const double w = surf_fit_normal * surf_fit_coeff->Eval(Tpr, ip);
3603  for (int i = 0; i < dof_x * dim; i++)
3604  {
3605  const int idof = i % dof_x, idim = i / dof_x;
3606  for (int j = 0; j <= i; j++)
3607  {
3608  const int jdof = j % dof_x, jdim = j / dof_x;
3609  const double entry =
3610  w * ( 2.0 * surf_fit_grad_s(idim) * shape_x(idof) *
3611  /* */ surf_fit_grad_s(jdim) * shape_x(jdof) +
3612  2.0 * sigma_e(s) * surf_fit_hess_s(idim, jdim) *
3613  /* */ shape_x(idof) * shape_x(jdof));
3614  mat(i, j) += entry;
3615  if (i != j) { mat(j, i) += entry; }
3616  }
3617  }
3618  }
3619 }
3620 
3621 double TMOP_Integrator::GetFDDerivative(const FiniteElement &el,
3622  ElementTransformation &T,
3623  Vector &elfun, const int dofidx,
3624  const int dir, const double e_fx,
3625  bool update_stored)
3626 {
3627  int dof = el.GetDof();
3628  int idx = dir*dof+dofidx;
3629  elfun[idx] += dx;
3630  double e_fxph = GetElementEnergy(el, T, elfun);
3631  elfun[idx] -= dx;
3632  double dfdx = (e_fxph-e_fx)/dx;
3633 
3634  if (update_stored)
3635  {
3636  (*(ElemPertEnergy[T.ElementNo]))(idx) = e_fxph;
3637  (*(ElemDer[T.ElementNo]))(idx) = dfdx;
3638  }
3639 
3640  return dfdx;
3641 }
3642 
3643 void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
3644  ElementTransformation &T,
3645  const Vector &elfun,
3646  Vector &elvect)
3647 {
3648  const int dof = el.GetDof(), dim = el.GetDim(), elnum = T.ElementNo;
3649  if (elnum>=ElemDer.Size())
3650  {
3651  ElemDer.Append(new Vector);
3652  ElemPertEnergy.Append(new Vector);
3653  ElemDer[elnum]->SetSize(dof*dim);
3654  ElemPertEnergy[elnum]->SetSize(dof*dim);
3655  }
3656 
3657  elvect.SetSize(dof*dim);
3658  Vector elfunmod(elfun);
3659 
3660  // In GetElementEnergy(), skip terms that have exact derivative calculations.
3661  fd_call_flag = true;
3662 
3663  // Energy for unperturbed configuration.
3664  const double e_fx = GetElementEnergy(el, T, elfun);
3665 
3666  for (int j = 0; j < dim; j++)
3667  {
3668  for (int i = 0; i < dof; i++)
3669  {
3670  if (discr_tc)
3671  {
3672  discr_tc->UpdateTargetSpecificationAtNode(
3673  el, T, i, j, discr_tc->GetTspecPert1H());
3674  }
3675  elvect(j*dof+i) = GetFDDerivative(el, T, elfunmod, i, j, e_fx, true);
3676  if (discr_tc) { discr_tc->RestoreTargetSpecificationAtNode(T, i); }
3677  }
3678  }
3679  fd_call_flag = false;
3680 
3681  // Contributions from adaptive limiting, surface fitting (exact derivatives).
3682  if (adapt_lim_gf || surf_fit_gf)
3683  {
3684  const IntegrationRule &ir = ActionIntegrationRule(el);
3685  const int nqp = ir.GetNPoints();
3686  DenseTensor Jtr(dim, dim, nqp);
3687  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3688 
3689  IsoparametricTransformation Tpr;
3690  Tpr.SetFE(&el);
3691  Tpr.ElementNo = T.ElementNo;
3692  Tpr.Attribute = T.Attribute;
3693  Tpr.mesh = T.mesh;
3694  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3695  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3696 
3697  Vector weights(nqp);
3698  for (int q = 0; q < nqp; q++)
3699  {
3700  weights(q) = ir.IntPoint(q).weight * Jtr(q).Det();
3701  }
3702 
3703  PMatO.UseExternalData(elvect.GetData(), dof, dim);
3704  if (adapt_lim_gf) { AssembleElemVecAdaptLim(el, Tpr, ir, weights, PMatO); }
3705  if (surf_fit_gf) { AssembleElemVecSurfFit(el, Tpr, PMatO); }
3706  }
3707 }
3708 
3709 void TMOP_Integrator::AssembleElementGradFD(const FiniteElement &el,
3710  ElementTransformation &T,
3711  const Vector &elfun,
3712  DenseMatrix &elmat)
3713 {
3714  const int dof = el.GetDof(), dim = el.GetDim();
3715 
3716  elmat.SetSize(dof*dim);
3717  Vector elfunmod(elfun);
3718 
3719  const Vector &ElemDerLoc = *(ElemDer[T.ElementNo]);
3720  const Vector &ElemPertLoc = *(ElemPertEnergy[T.ElementNo]);
3721 
3722  // In GetElementEnergy(), skip terms that have exact derivative calculations.
3723  fd_call_flag = true;
3724  for (int i = 0; i < dof; i++)
3725  {
3726  for (int j = 0; j < i+1; j++)
3727  {
3728  for (int k1 = 0; k1 < dim; k1++)
3729  {
3730  for (int k2 = 0; k2 < dim; k2++)
3731  {
3732  elfunmod(k2*dof+j) += dx;
3733 
3734  if (discr_tc)
3735  {
3736  discr_tc->UpdateTargetSpecificationAtNode(
3737  el, T, j, k2, discr_tc->GetTspecPert1H());
3738  if (j != i)
3739  {
3740  discr_tc->UpdateTargetSpecificationAtNode(
3741  el, T, i, k1, discr_tc->GetTspecPert1H());
3742  }
3743  else // j==i
3744  {
3745  if (k1 != k2)
3746  {
3747  int idx = k1+k2-1;
3748  discr_tc->UpdateTargetSpecificationAtNode(
3749  el, T, i, idx, discr_tc->GetTspecPertMixH());
3750  }
3751  else // j==i && k1==k2
3752  {
3753  discr_tc->UpdateTargetSpecificationAtNode(
3754  el, T, i, k1, discr_tc->GetTspecPert2H());
3755  }
3756  }
3757  }
3758 
3759  double e_fx = ElemPertLoc(k2*dof+j);
3760  double e_fpxph = GetFDDerivative(el, T, elfunmod, i, k1, e_fx,
3761  false);
3762  elfunmod(k2*dof+j) -= dx;
3763  double e_fpx = ElemDerLoc(k1*dof+i);
3764 
3765  elmat(k1*dof+i, k2*dof+j) = (e_fpxph - e_fpx) / dx;
3766  elmat(k2*dof+j, k1*dof+i) = (e_fpxph - e_fpx) / dx;
3767 
3768  if (discr_tc)
3769  {
3770  discr_tc->RestoreTargetSpecificationAtNode(T, i);
3771  discr_tc->RestoreTargetSpecificationAtNode(T, j);
3772  }
3773  }
3774  }
3775  }
3776  }
3777  fd_call_flag = false;
3778 
3779  // Contributions from adaptive limiting.
3780  if (adapt_lim_gf || surf_fit_gf)
3781  {
3782  const IntegrationRule &ir = GradientIntegrationRule(el);
3783  const int nqp = ir.GetNPoints();
3784  DenseTensor Jtr(dim, dim, nqp);
3785  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3786 
3787  IsoparametricTransformation Tpr;
3788  Tpr.SetFE(&el);
3789  Tpr.ElementNo = T.ElementNo;
3790  Tpr.Attribute = T.Attribute;
3791  Tpr.mesh = T.mesh;
3792  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3793  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3794 
3795  Vector weights(nqp);
3796  for (int q = 0; q < nqp; q++)
3797  {
3798  weights(q) = ir.IntPoint(q).weight * Jtr(q).Det();
3799  }
3800 
3801  if (adapt_lim_gf) { AssembleElemGradAdaptLim(el, Tpr, ir, weights, elmat); }
3802  if (surf_fit_gf) { AssembleElemGradSurfFit(el, Tpr, elmat); }
3803  }
3804 }
3805 
3806 void TMOP_Integrator::UpdateSurfaceFittingWeight(double factor)
3807 {
3808  if (!surf_fit_coeff) { return; }
3809 
3810  if (surf_fit_coeff)
3811  {
3812  auto cf = dynamic_cast<ConstantCoefficient *>(surf_fit_coeff);
3813  MFEM_VERIFY(cf, "Dynamic weight works only with a ConstantCoefficient.");
3814  cf->constant *= factor;
3815  }
3816 }
3817 
3818 double TMOP_Integrator::GetSurfaceFittingWeight()
3819 {
3820  if (surf_fit_coeff)
3821  {
3822  auto cf = dynamic_cast<ConstantCoefficient *>(surf_fit_coeff);
3823  MFEM_VERIFY(cf, "Dynamic weight works only with a ConstantCoefficient.");
3824  return cf->constant;
3825  }
3826  return 0.0;
3827 }
3828 
3829 void TMOP_Integrator::EnableNormalization(const GridFunction &x)
3830 {
3831  ComputeNormalizationEnergies(x, metric_normal, lim_normal, surf_fit_normal);
3832  metric_normal = 1.0 / metric_normal;
3833  lim_normal = 1.0 / lim_normal;
3834  //if (surf_fit_gf) { surf_fit_normal = 1.0 / surf_fit_normal; }
3835  if (surf_fit_gf) { surf_fit_normal = lim_normal; }
3836 }
3837 
3838 #ifdef MFEM_USE_MPI
3839 void TMOP_Integrator::ParEnableNormalization(const ParGridFunction &x)
3840 {
3841  double loc[3];
3842  ComputeNormalizationEnergies(x, loc[0], loc[1], loc[2]);
3843  double rdc[3];
3844  MPI_Allreduce(loc, rdc, 3, MPI_DOUBLE, MPI_SUM, x.ParFESpace()->GetComm());
3845  metric_normal = 1.0 / rdc[0];
3846  lim_normal = 1.0 / rdc[1];
3847  // if (surf_fit_gf) { surf_fit_normal = 1.0 / rdc[2]; }
3848  if (surf_fit_gf) { surf_fit_normal = lim_normal; }
3849 }
3850 #endif
3851 
3852 void TMOP_Integrator::ComputeNormalizationEnergies(const GridFunction &x,
3853  double &metric_energy,
3854  double &lim_energy,
3855  double &surf_fit_gf_energy)
3856 {
3857  Array<int> vdofs;
3858  Vector x_vals;
3859  const FiniteElementSpace* const fes = x.FESpace();
3860 
3861  const int dim = fes->GetMesh()->Dimension();
3862  Jrt.SetSize(dim);
3863  Jpr.SetSize(dim);
3864  Jpt.SetSize(dim);
3865 
3866  metric_energy = 0.0;
3867  lim_energy = 0.0;
3868  surf_fit_gf_energy = 0.0;
3869  for (int i = 0; i < fes->GetNE(); i++)
3870  {
3871  const FiniteElement *fe = fes->GetFE(i);
3872  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
3873  const int nqp = ir.GetNPoints();
3874  DenseTensor Jtr(dim, dim, nqp);
3875  const int dof = fe->GetDof();
3876  DSh.SetSize(dof, dim);
3877 
3878  fes->GetElementVDofs(i, vdofs);
3879  x.GetSubVector(vdofs, x_vals);
3880  PMatI.UseExternalData(x_vals.GetData(), dof, dim);
3881 
3882  targetC->ComputeElementTargets(i, *fe, ir, x_vals, Jtr);
3883 
3884  for (int q = 0; q < nqp; q++)
3885  {
3886  const IntegrationPoint &ip = ir.IntPoint(q);
3887  metric->SetTargetJacobian(Jtr(q));
3888  CalcInverse(Jtr(q), Jrt);
3889  const double weight = ip.weight * Jtr(q).Det();
3890 
3891  fe->CalcDShape(ip, DSh);
3892  MultAtB(PMatI, DSh, Jpr);
3893  Mult(Jpr, Jrt, Jpt);
3894 
3895  metric_energy += weight * metric->EvalW(Jpt);
3896  lim_energy += weight;
3897  }
3898 
3899  // Normalization of the surface fitting term.
3900  if (surf_fit_gf)
3901  {
3902  Array<int> dofs;
3903  Vector sigma_e;
3904  surf_fit_gf->FESpace()->GetElementDofs(i, dofs);
3905  surf_fit_gf->GetSubVector(dofs, sigma_e);
3906  for (int s = 0; s < dofs.Size(); s++)
3907  {
3908  if ((*surf_fit_marker)[dofs[s]] == true)
3909  {
3910  surf_fit_gf_energy += sigma_e(s) * sigma_e(s);
3911  }
3912  }
3913  }
3914  }
3915 
3916  if (targetC->ContainsVolumeInfo() == false)
3917  {
3918  // Special case when the targets don't contain volumetric information.
3919  lim_energy = fes->GetNE();
3920  }
3921 }
3922 
3923 void TMOP_Integrator::ComputeMinJac(const Vector &x,
3924  const FiniteElementSpace &fes)
3925 {
3926  const FiniteElement *fe = fes.GetFE(0);
3927  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
3928  const int NE = fes.GetMesh()->GetNE(), dim = fe->GetDim(),
3929  dof = fe->GetDof(), nsp = ir.GetNPoints();
3930 
3931  Array<int> xdofs(dof * dim);
3932  DenseMatrix dshape(dof, dim), pos(dof, dim);
3933  Vector posV(pos.Data(), dof * dim);
3934  Jpr.SetSize(dim);
3935 
3936  dx = std::numeric_limits<float>::max();
3937 
3938  double detv_sum;
3939  double detv_avg_min = std::numeric_limits<float>::max();
3940  for (int i = 0; i < NE; i++)
3941  {
3942  fes.GetElementVDofs(i, xdofs);
3943  x.GetSubVector(xdofs, posV);
3944  detv_sum = 0.;
3945  for (int j = 0; j < nsp; j++)
3946  {
3947  fes.GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
3948  MultAtB(pos, dshape, Jpr);
3949  detv_sum += std::fabs(Jpr.Det());
3950  }
3951  double detv_avg = pow(detv_sum/nsp, 1./dim);
3952  detv_avg_min = std::min(detv_avg, detv_avg_min);
3953  }
3954  dx = detv_avg_min / dxscale;
3955 }
3956 
3957 void TMOP_Integrator::UpdateAfterMeshPositionChange(const Vector &new_x,
3958  int new_x_ordering)
3959 {
3960  if (discr_tc)
3961  {
3962  PA.Jtr_needs_update = true;
3963  }
3964  // Update adapt_lim_gf if adaptive limiting is enabled.
3965  if (adapt_lim_gf)
3966  {
3967  adapt_lim_eval->ComputeAtNewPosition(new_x, *adapt_lim_gf, new_x_ordering);
3968  }
3969 
3970  // Update surf_fit_gf if surface fitting is enabled.
3971  if (surf_fit_gf)
3972  {
3973  surf_fit_eval->ComputeAtNewPosition(new_x, *surf_fit_gf, new_x_ordering);
3974  }
3975 }
3976 
3977 void TMOP_Integrator::ComputeFDh(const Vector &x, const FiniteElementSpace &fes)
3978 {
3979  if (!fdflag) { return; }
3980  ComputeMinJac(x, fes);
3981 #ifdef MFEM_USE_MPI
3982  const ParFiniteElementSpace *pfes =
3983  dynamic_cast<const ParFiniteElementSpace *>(&fes);
3984  if (pfes)
3985  {
3986  double min_jac_all;
3987  MPI_Allreduce(&dx, &min_jac_all, 1, MPI_DOUBLE, MPI_MIN, pfes->GetComm());
3988  dx = min_jac_all;
3989  }
3990 #endif
3991 }
3992 
3993 void TMOP_Integrator::EnableFiniteDifferences(const GridFunction &x)
3994 {
3995  fdflag = true;
3996  const FiniteElementSpace *fes = x.FESpace();
3997  ComputeFDh(x,*fes);
3998  if (discr_tc)
3999  {
4000  discr_tc->UpdateTargetSpecification(x, false, fes->GetOrdering());
4001  discr_tc->UpdateGradientTargetSpecification(x, dx, false, fes->GetOrdering());
4002  discr_tc->UpdateHessianTargetSpecification(x, dx, false, fes->GetOrdering());
4003  }
4004 }
4005 
4006 #ifdef MFEM_USE_MPI
4007 void TMOP_Integrator::EnableFiniteDifferences(const ParGridFunction &x)
4008 {
4009  fdflag = true;
4010  const ParFiniteElementSpace *pfes = x.ParFESpace();
4011  ComputeFDh(x,*pfes);
4012  if (discr_tc)
4013  {
4014  discr_tc->UpdateTargetSpecification(x, false, pfes->GetOrdering());
4015  discr_tc->UpdateGradientTargetSpecification(x, dx, false, pfes->GetOrdering());
4016  discr_tc->UpdateHessianTargetSpecification(x, dx, false, pfes->GetOrdering());
4017  }
4018 }
4019 #endif
4020 
4021 double TMOP_Integrator::ComputeMinDetT(const Vector &x,
4022  const FiniteElementSpace &fes)
4023 {
4024  double min_detT = std::numeric_limits<double>::infinity();
4025  const int NE = fes.GetMesh()->GetNE();
4026  const int dim = fes.GetMesh()->Dimension();
4027  Array<int> xdofs;
4028  Jpr.SetSize(dim);
4029  Jpt.SetSize(dim);
4030  Jrt.SetSize(dim);
4031 
4032  for (int i = 0; i < NE; i++)
4033  {
4034  const FiniteElement *fe = fes.GetFE(i);
4035  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
4036  const int dof = fe->GetDof(), nsp = ir.GetNPoints();
4037 
4038  DSh.SetSize(dof, dim);
4039  Vector posV(dof * dim);
4040  PMatI.UseExternalData(posV.GetData(), dof, dim);
4041 
4042  fes.GetElementVDofs(i, xdofs);
4043  x.GetSubVector(xdofs, posV);
4044 
4045  DenseTensor Jtr(dim, dim, ir.GetNPoints());
4046  targetC->ComputeElementTargets(i, *fe, ir, posV, Jtr);
4047 
4048  for (int q = 0; q < nsp; q++)
4049  {
4050  const IntegrationPoint &ip = ir.IntPoint(q);
4051  const DenseMatrix &Jtr_q = Jtr(q);
4052  CalcInverse(Jtr_q, Jrt);
4053  fe->CalcDShape(ip, DSh);
4054  MultAtB(PMatI, DSh, Jpr);
4055  Mult(Jpr, Jrt, Jpt);
4056  double detT = Jpt.Det();
4057  min_detT = std::min(min_detT, detT);
4058  }
4059  }
4060  return min_detT;
4061 }
4062 
4063 double TMOP_Integrator::ComputeUntanglerMaxMuBarrier(const Vector &x,
4064  const FiniteElementSpace &fes)
4065 {
4066  double max_muT = -std::numeric_limits<double>::infinity();
4067  const int NE = fes.GetMesh()->GetNE();
4068  const int dim = fes.GetMesh()->Dimension();
4069  Array<int> xdofs;
4070  Jpr.SetSize(dim);
4071  Jpt.SetSize(dim);
4072  Jrt.SetSize(dim);
4073 
4074  TMOP_WorstCaseUntangleOptimizer_Metric *wcuo =
4075  dynamic_cast<TMOP_WorstCaseUntangleOptimizer_Metric *>(metric);
4076 
4077  if (!wcuo || wcuo->GetWorstCaseType() !=
4078  TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::Beta)
4079  {
4080  return 0.0;
4081  }
4082 
4083  for (int i = 0; i < NE; i++)
4084  {
4085  const FiniteElement *fe = fes.GetFE(i);
4086  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
4087  const int dof = fe->GetDof(), nsp = ir.GetNPoints();
4088  Jpr.SetSize(dim);
4089  Jrt.SetSize(dim);
4090  Jpt.SetSize(dim);
4091 
4092  DSh.SetSize(dof, dim);
4093  Vector posV(dof * dim);
4094  PMatI.UseExternalData(posV.GetData(), dof, dim);
4095 
4096  fes.GetElementVDofs(i, xdofs);
4097  x.GetSubVector(xdofs, posV);
4098 
4099  DenseTensor Jtr(dim, dim, ir.GetNPoints());
4100  targetC->ComputeElementTargets(i, *fe, ir, posV, Jtr);
4101 
4102  for (int q = 0; q < nsp; q++)
4103  {
4104  const IntegrationPoint &ip = ir.IntPoint(q);
4105  const DenseMatrix &Jtr_q = Jtr(q);
4106  CalcInverse(Jtr_q, Jrt);
4107 
4108  fe->CalcDShape(ip, DSh);
4109  MultAtB(PMatI, DSh, Jpr);
4110  Mult(Jpr, Jrt, Jpt);
4111 
4112  double metric_val = 0.0;
4113  if (wcuo)
4114  {
4115  wcuo->SetTargetJacobian(Jtr_q);
4116  metric_val = wcuo->EvalWBarrier(Jpt);
4117  }
4118 
4119  max_muT = std::max(max_muT, metric_val);
4120  }
4121  }
4122  return max_muT;
4123 }
4124 
4125 void TMOP_Integrator::ComputeUntangleMetricQuantiles(const Vector &x,
4126  const FiniteElementSpace &fes)
4127 {
4128  TMOP_WorstCaseUntangleOptimizer_Metric *wcuo =
4129  dynamic_cast<TMOP_WorstCaseUntangleOptimizer_Metric *>(metric);
4130 
4131  if (!wcuo) { return; }
4132 
4133 #ifdef MFEM_USE_MPI
4134  const ParFiniteElementSpace *pfes =
4135  dynamic_cast<const ParFiniteElementSpace *>(&fes);
4136 #endif
4137 
4138  if (wcuo && wcuo->GetBarrierType() ==
4139  TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Shifted)
4140  {
4141  double min_detT = ComputeMinDetT(x, fes);
4142  double min_detT_all = min_detT;
4143 #ifdef MFEM_USE_MPI
4144  if (pfes)
4145  {
4146  MPI_Allreduce(&min_detT, &min_detT_all, 1, MPI_DOUBLE, MPI_MIN,
4147  pfes->GetComm());
4148  }
4149 #endif
4150  if (wcuo) { wcuo->SetMinDetT(min_detT_all); }
4151  }
4152 
4153  double max_muT = ComputeUntanglerMaxMuBarrier(x, fes);
4154  double max_muT_all = max_muT;
4155 #ifdef MFEM_USE_MPI
4156  if (pfes)
4157  {
4158  MPI_Allreduce(&max_muT, &max_muT_all, 1, MPI_DOUBLE, MPI_MAX,
4159  pfes->GetComm());
4160  }
4161 #endif
4162  wcuo->SetMaxMuT(max_muT_all);
4163 }
4164 
4165 void TMOPComboIntegrator::EnableLimiting(const GridFunction &n0,
4166  const GridFunction &dist,
4167  Coefficient &w0,
4168  TMOP_LimiterFunction *lfunc)
4169 {
4170  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
4171 
4172  tmopi[0]->EnableLimiting(n0, dist, w0, lfunc);
4173  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
4174 }
4175 
4176 void TMOPComboIntegrator::EnableLimiting(const GridFunction &n0,
4177  Coefficient &w0,
4178  TMOP_LimiterFunction *lfunc)
4179 {
4180  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
4181 
4182  tmopi[0]->EnableLimiting(n0, w0, lfunc);
4183  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
4184 }
4185 
4186 void TMOPComboIntegrator::SetLimitingNodes(const GridFunction &n0)
4187 {
4188  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
4189 
4190  tmopi[0]->SetLimitingNodes(n0);
4191  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
4192 }
4193 
4194 double TMOPComboIntegrator::GetElementEnergy(const FiniteElement &el,
4195  ElementTransformation &T,
4196  const Vector &elfun)
4197 {
4198  double energy= 0.0;
4199  for (int i = 0; i < tmopi.Size(); i++)
4200  {
4201  energy += tmopi[i]->GetElementEnergy(el, T, elfun);
4202  }
4203  return energy;
4204 }
4205 
4206 void TMOPComboIntegrator::AssembleElementVector(const FiniteElement &el,
4207  ElementTransformation &T,
4208  const Vector &elfun,
4209  Vector &elvect)
4210 {
4211  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
4212 
4213  tmopi[0]->AssembleElementVector(el, T, elfun, elvect);
4214  for (int i = 1; i < tmopi.Size(); i++)
4215  {
4216  Vector elvect_i;
4217  tmopi[i]->AssembleElementVector(el, T, elfun, elvect_i);
4218  elvect += elvect_i;
4219  }
4220 }
4221 
4222 void TMOPComboIntegrator::AssembleElementGrad(const FiniteElement &el,
4223  ElementTransformation &T,
4224  const Vector &elfun,
4225  DenseMatrix &elmat)
4226 {
4227  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
4228 
4229  tmopi[0]->AssembleElementGrad(el, T, elfun, elmat);
4230  for (int i = 1; i < tmopi.Size(); i++)
4231  {
4232  DenseMatrix elmat_i;
4233  tmopi[i]->AssembleElementGrad(el, T, elfun, elmat_i);
4234  elmat += elmat_i;
4235  }
4236 }
4237 
4238 double TMOPComboIntegrator::GetRefinementElementEnergy(const FiniteElement &el,
4239  ElementTransformation &T,
4240  const Vector &elfun,
4241  const IntegrationRule &irule)
4242 {
4243  double energy= 0.0;
4244  for (int i = 0; i < tmopi.Size(); i++)
4245  {
4246  energy += tmopi[i]->GetRefinementElementEnergy(el, T, elfun, irule);
4247  }
4248  return energy;
4249 }
4250 
4251 double TMOPComboIntegrator::GetDerefinementElementEnergy(
4252  const FiniteElement &el,
4253  ElementTransformation &T,
4254  const Vector &elfun)
4255 {
4256  double energy= 0.0;
4257  for (int i = 0; i < tmopi.Size(); i++)
4258  {
4259  energy += tmopi[i]->GetDerefinementElementEnergy(el, T, elfun);
4260  }
4261  return energy;
4262 }
4263 
4264 void TMOPComboIntegrator::EnableNormalization(const GridFunction &x)
4265 {
4266  const int cnt = tmopi.Size();
4267  double total_integral = 0.0;
4268  for (int i = 0; i < cnt; i++)
4269  {
4270  tmopi[i]->EnableNormalization(x);
4271  total_integral += 1.0 / tmopi[i]->metric_normal;
4272  }
4273  for (int i = 0; i < cnt; i++)
4274  {
4275  tmopi[i]->metric_normal = 1.0 / total_integral;
4276  }
4277 }
4278 
4279 #ifdef MFEM_USE_MPI
4280 void TMOPComboIntegrator::ParEnableNormalization(const ParGridFunction &x)
4281 {
4282  const int cnt = tmopi.Size();
4283  double total_integral = 0.0;
4284  for (int i = 0; i < cnt; i++)
4285  {
4286  tmopi[i]->ParEnableNormalization(x);
4287  total_integral += 1.0 / tmopi[i]->metric_normal;
4288  }
4289  for (int i = 0; i < cnt; i++)
4290  {
4291  tmopi[i]->metric_normal = 1.0 / total_integral;
4292  }
4293 }
4294 #endif
4295 
4296 void TMOPComboIntegrator::AssemblePA(const FiniteElementSpace &fes)
4297 {
4298  for (int i = 0; i < tmopi.Size(); i++)
4299  {
4300  tmopi[i]->AssemblePA(fes);
4301  }
4302 }
4303 
4304 void TMOPComboIntegrator::AssembleGradPA(const Vector &xe,
4305  const FiniteElementSpace &fes)
4306 {
4307  for (int i = 0; i < tmopi.Size(); i++)
4308  {
4309  tmopi[i]->AssembleGradPA(xe,fes);
4310  }
4311 }
4312 
4313 void TMOPComboIntegrator::AssembleGradDiagonalPA(Vector &de) const
4314 {
4315  for (int i = 0; i < tmopi.Size(); i++)
4316  {
4317  tmopi[i]->AssembleGradDiagonalPA(de);
4318  }
4319 }
4320 
4321 void TMOPComboIntegrator::AddMultPA(const Vector &xe, Vector &ye) const
4322 {
4323  for (int i = 0; i < tmopi.Size(); i++)
4324  {
4325  tmopi[i]->AddMultPA(xe, ye);
4326  }
4327 }
4328 
4329 void TMOPComboIntegrator::AddMultGradPA(const Vector &re, Vector &ce) const
4330 {
4331  for (int i = 0; i < tmopi.Size(); i++)
4332  {
4333  tmopi[i]->AddMultGradPA(re, ce);
4334  }
4335 }
4336 
4337 double TMOPComboIntegrator::GetLocalStateEnergyPA(const Vector &xe) const
4338 {
4339  double energy = 0.0;
4340  for (int i = 0; i < tmopi.Size(); i++)
4341  {
4342  energy += tmopi[i]->GetLocalStateEnergyPA(xe);
4343  }
4344  return energy;
4345 }
4346 
4347 void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric,
4348  const TargetConstructor &tc,
4349  const Mesh &mesh, GridFunction &metric_gf)
4350 {
4351  const int NE = mesh.GetNE();
4352  const GridFunction &nodes = *mesh.GetNodes();
4353  const int dim = mesh.Dimension();
4354  DenseMatrix Winv(dim), T(dim), A(dim), dshape, pos;
4355  Array<int> pos_dofs, gf_dofs;
4356  DenseTensor W;
4357  Vector posV;
4358 
4359  for (int i = 0; i < NE; i++)
4360  {
4361  const FiniteElement &fe_pos = *nodes.FESpace()->GetFE(i);
4362  const IntegrationRule &ir = metric_gf.FESpace()->GetFE(i)->GetNodes();
4363  const int nsp = ir.GetNPoints(), dof = fe_pos.GetDof();
4364 
4365  dshape.SetSize(dof, dim);
4366  pos.SetSize(dof, dim);
4367  posV.SetDataAndSize(pos.Data(), dof * dim);
4368 
4369  metric_gf.FESpace()->GetElementDofs(i, gf_dofs);
4370  nodes.FESpace()->GetElementVDofs(i, pos_dofs);
4371  nodes.GetSubVector(pos_dofs, posV);
4372 
4373  W.SetSize(dim, dim, nsp);
4374  tc.ComputeElementTargets(i, fe_pos, ir, posV, W);
4375 
4376  for (int j = 0; j < nsp; j++)
4377  {
4378  const DenseMatrix &Wj = W(j);
4379  metric.SetTargetJacobian(Wj);
4380  CalcInverse(Wj, Winv);
4381 
4382  const IntegrationPoint &ip = ir.IntPoint(j);
4383  fe_pos.CalcDShape(ip, dshape);
4384  MultAtB(pos, dshape, A);
4385  Mult(A, Winv, T);
4386 
4387  metric_gf(gf_dofs[j]) = metric.EvalW(T);
4388  }
4389  }
4390 }
4391 
4392 } // namespace mfem
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
mfem::FiniteElement
Abstract class for all finite elements.
Definition: fe_base.hpp:235
mfem::TMOP_Integrator::ComputeUntangleMetricQuantiles
void ComputeUntangleMetricQuantiles(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:4125
mfem::TMOP_Metric_077::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:478
mfem::DenseMatrix::Size
int Size() const
For backward compatibility define Size to be synonym of Width()
Definition: densemat.hpp:102
mfem::TMOP_Integrator::GetDerefinementElementEnergy
virtual double GetDerefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Definition: tmop.cpp:3064
mfem::TMOP_Metric_022::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:416
mfem::FiniteElementSpace::GetOrdering
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
mfem::TMOP_QualityMetric::Id
virtual int Id() const
Return the metric ID.
Definition: tmop.hpp:78
mfem::DiscreteAdaptTC::SetParDiscreteTargetSize
virtual void SetParDiscreteTargetSize(const ParGridFunction &tspec_)
Definition: tmop.cpp:1689
mfem::TMOP_Metric_001::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:116
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
mfem::TMOP_AMetric_014a::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1277
mfem::TMOP_Metric_022::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:386
mfem::TMOP_Metric_311::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:671
mfem::AdaptivityEvaluator::SetInitialField
virtual void SetInitialField(const Vector &init_nodes, const Vector &init_field)=0
mfem::FiniteElement::GetDim
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
mfem::TMOP_QualityMetric::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const =0
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
mfem::AddMultVWt
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3127
mfem::FiniteElementSpace::GetNDofs
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:584
mfem::TMOP_Metric_322::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:768
mfem::TMOP_Integrator::AssembleElemVecSurfFit
void AssembleElemVecSurfFit(const FiniteElement &el_x, IsoparametricTransformation &Tpr, DenseMatrix &mat)
Definition: tmop.cpp:3491
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
mfem::AnalyticAdaptTC::vector_tspec
VectorCoefficient * vector_tspec
Definition: tmop.hpp:1311
mfem::TMOP_Metric_002::EvalWMatrixForm
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.cpp:245
mfem::GridFunction::GetVectorValues
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:706
mfem::TMOP_Integrator::EnableFiniteDifferences
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:3993
mfem::Vector::NewDataAndSize
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:163
mfem::TargetConstructor::target_type
const TargetType target_type
Definition: tmop.hpp:1193
mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
mfem::TMOP_Integrator::PMatO
DenseMatrix PMatO
Definition: tmop.hpp:1638
mfem::TMOPComboIntegrator::AddMultPA
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
Definition: tmop.cpp:4321
mfem::TMOP_Integrator::Jtr
DenseTensor Jtr
Definition: tmop.hpp:1668
mfem::DenseTensor::UseExternalData
void UseExternalData(double *ext_data, int i, int j, int k)
Definition: densemat.hpp:1014
mfem::TMOP_Metric_022::min_detT
double & min_detT
Definition: tmop.hpp:366
mfem::FiniteElementSpace::Update
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3334
mfem::TMOP_Metric_313::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:690
mfem::TMOP_QualityMetric::Jtr
const DenseMatrix * Jtr
Definition: tmop.hpp:26
mfem::TMOP_Metric_321::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:747
mfem::FiniteElementSpace::IsVariableOrder
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:463
mfem::InvariantsEvaluator2D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:384
mfem::IntegrationRules::Get
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
mfem::TMOP_Metric_302::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:776
mfem::TMOP_Integrator::GetSurfaceFittingWeight
double GetSurfaceFittingWeight()
Get the surface fitting weight.
Definition: tmop.cpp:3818
mfem::ParFiniteElementSpace::GetParMesh
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::detT_ep
double detT_ep
Definition: tmop.hpp:149
mfem::VectorCoefficient
Base class for vector Coefficients that optionally depend on time and space.
Definition: coefficient.hpp:448
mfem::ConstantCoefficient
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
mfem::InvariantsEvaluator3D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:1090
mfem::TMOP_Integrator::Jrt
DenseMatrix Jrt
Definition: tmop.hpp:1638
mfem::InvariantsEvaluator2D::SetDerivativeMatrix
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 2, using column-major storage.
Definition: invariants.hpp:184
mfem::TMOP_Integrator::lim_dist
const GridFunction * lim_dist
Definition: tmop.hpp:1590
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetOrientation
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
Definition: tmop.cpp:1802
mfem::TMOPComboIntegrator::GetElementEnergy
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Compute the local energy.
Definition: tmop.cpp:4194
mfem::TMOP_Integrator::PA
struct mfem::TMOP_Integrator::@23 PA
mfem::AdaptivityEvaluator::ComputeAtNewPosition
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field, int new_nodes_ordering=Ordering::byNODES)=0
mfem::TMOPMatrixCoefficient::EvalGrad
virtual void EvalGrad(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip, int comp)=0
Evaluate the derivative of the matrix coefficient with respect to comp in the element described by T ...
mfem::Mesh::GetNumGeometries
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:5908
mfem::TMOP_Metric_055::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:401
mfem::TMOP_Metric_311::eps
const double eps
Definition: tmop.hpp:670
mfem::TMOP_Integrator::EnableNormalization
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
Definition: tmop.cpp:3829
mfem::InvariantsEvaluator3D::Get_dI3b
const scalar_t * Get_dI3b()
Definition: invariants.hpp:783
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
mfem::InterpolateTMOP_QualityMetric
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:4347
mfem::TMOP_Integrator::metric
TMOP_QualityMetric * metric
Definition: tmop.hpp:1572
mfem::TMOP_Metric_056::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:418
mfem::DenseMatrix::Det
double Det() const
Definition: densemat.cpp:449
mfem::TMOP_LimiterFunction
Base class for limiting functions to be used in class TMOP_Integrator.
Definition: tmop.hpp:1057
mfem::TMOP_QualityMetric::EvalW
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...
mfem::TMOP_AMetric_011::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1251
mfem::TMOP_LimiterFunction::Eval_d2
virtual void Eval_d2(const Vector &x, const Vector &x0, double dist, DenseMatrix &d2) const =0
Returns the Hessian of the limiting function f(x, x0, d) with respect to x.
mfem::DiscreteAdaptTC::ComputeElementTargetsGradient
virtual void ComputeElementTargetsGradient(const IntegrationRule &ir, const Vector &elfun, IsoparametricTransformation &Tpr, DenseTensor &dJtr) const
Definition: tmop.cpp:2164
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
mfem::TMOP_Metric_321::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:1030
mfem::Operator::Width
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
mfem::TMOP_Metric_315::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:941
mfem::TMOP_Integrator::lim_normal
double lim_normal
Definition: tmop.hpp:1594
mfem::ParFiniteElementSpace::Update
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3287
mfem::Vector::Norml2
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:812
mfem::DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange
void ParUpdateAfterMeshTopologyChange()
Definition: tmop.cpp:1658
mfem::TMOP_Metric_002::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:250
mfem::TMOP_Integrator::Jpr
DenseMatrix Jpr
Definition: tmop.hpp:1638
mfem::MatrixCoefficient::Eval
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
mfem::TMOP_Metric_252::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:667
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
mfem::Vector::GetSubVector
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:547
mfem::TMOP_AMetric_107a::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1309
mfem::TMOP_QualityMetric::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
mfem::TMOP_Metric_001::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:200
mfem::AnalyticAdaptTC::ComputeElementTargetsGradient
virtual void ComputeElementTargetsGradient(const IntegrationRule &ir, const Vector &elfun, IsoparametricTransformation &Tpr, DenseTensor &dJtr) const
Definition: tmop.cpp:1586
mfem::AdaptivityEvaluator::pfes
ParFiniteElementSpace * pfes
Definition: tmop.hpp:1122
mfem::TMOP_Metric_302::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:766
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
mfem::TMOP_Integrator::GetSurfaceFittingErrors
void GetSurfaceFittingErrors(double &err_avg, double &err_max)
Definition: tmop.cpp:2786
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
mfem::TMOP_Metric_322::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1066
mfem::TMOP_Integrator::metric_normal
double metric_normal
Definition: tmop.hpp:1582
mfem::petsc::Mat
struct::_p_Mat * Mat
Definition: petsc.hpp:70
mfem::TMOP_Integrator::EnableAdaptiveLimiting
void EnableAdaptiveLimiting(const GridFunction &z0, Coefficient &coeff, AdaptivityEvaluator &ae)
Restriction of the node positions to certain regions.
Definition: tmop.cpp:2715
mfem::TMOP_WorstCaseUntangleOptimizer_Metric
Definition: tmop.hpp:127
mfem::TMOP_Metric_311::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:885
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
mfem::TMOPComboIntegrator::EnableNormalization
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
Definition: tmop.cpp:4264
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::tmop_metric
TMOP_QualityMetric & tmop_metric
Definition: tmop.hpp:144
mfem::TMOP_Metric_252::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:679
mfem::TMOP_Metric_313::min_detT
double & min_detT
Definition: tmop.hpp:689
mfem::Operator::Mult
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
mfem::TMOP_Metric_085::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:603
mfem::ParFiniteElementSpace
Abstract parallel finite element space.
Definition: pfespace.hpp:28
mfem::TMOPComboIntegrator::AssemblePA
virtual void AssemblePA(const FiniteElementSpace &)
Method defining partial assembly.
Definition: tmop.cpp:4296
mfem::IsoparametricTransformation::GetFE
const FiniteElement * GetFE() const
Get the current element used to compute the transformations.
Definition: eltrans.hpp:389
mfem::TMOP_Metric_009::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:335
mfem::AdaptivityEvaluator::SetSerialMetaInfo
void SetSerialMetaInfo(const Mesh &m, const FiniteElementSpace &f)
Definition: tmop.cpp:2623
mfem::DiscreteAdaptTC::UpdateGradientTargetSpecification
void UpdateGradientTargetSpecification(const Vector &x, double dx, bool use_flag=false, int x_ordering=Ordering::byNODES)
Definition: tmop.cpp:2511
mfem::TMOP_Integrator::ir
const IntegrationRule * ir
Definition: tmop.hpp:1676
mfem::TMOP_Metric_058::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:436
mfem::TMOP_Integrator::AssembleElementVectorFD
void AssembleElementVectorFD(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Definition: tmop.cpp:3643
mfem::TargetConstructor::nodes
const GridFunction * nodes
Definition: tmop.hpp:1190
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSkew
virtual void SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
Definition: tmop.cpp:1782
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::EvalWBarrier
virtual double EvalWBarrier(const DenseMatrix &Jpt) const
Definition: tmop.cpp:88
mfem::TMOP_Integrator::UpdateSurfaceFittingWeight
void UpdateSurfaceFittingWeight(double factor)
Update the surface fitting weight as surf_fit_coeff *= factor;.
Definition: tmop.cpp:3806
mfem::DiscreteAdaptTC::SetRefinementSubElement
void SetRefinementSubElement(int amr_el_)
Definition: tmop.hpp:1553
mfem::InvariantsEvaluator2D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition: invariants.hpp:432
mfem::InvariantsEvaluator2D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:226
mfem::DenseMatrix::GetData
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
mfem::IsoparametricTransformation::OrderJ
virtual int OrderJ() const
Return the order of the elements of the Jacobian of the transformation.
Definition: eltrans.cpp:443
mfem::FiniteElement::ProjectGrad
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:160
mfem::Vector::GetData
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
mfem::TMOP_Metric_321::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:1019
mfem::TMOP_QuadraticLimiter
Default limiter function in TMOP_Integrator.
Definition: tmop.hpp:1078
mfem::Mesh::GetElementBaseGeometry
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1069
mfem::Vector::UseDevice
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:118
mfem::DenseMatrix::Set
void Set(double alpha, const double *A)
Set the matrix to alpha * A, assuming that A has the same dimensions as the matrix and uses column-ma...
Definition: densemat.cpp:533
mfem::DiscreteAdaptTC::ComputeElementTargets
virtual void ComputeElementTargets(int e_id, const FiniteElement &fe, const IntegrationRule &ir, const Vector &elfun, DenseTensor &Jtr) const
Given an element and quadrature rule, computes ref-&gt;target transformation Jacobians for each quadratur...
Definition: tmop.cpp:1949
mfem::TMOP_Integrator::AssembleElementGradExact
void AssembleElementGradExact(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Definition: tmop.cpp:3286
mfem::TMOP_Integrator::ElemDer
Array< Vector * > ElemDer
Definition: tmop.hpp:1625
mfem::DiscreteAdaptTC::GetTspecPertMixH
const Vector & GetTspecPertMixH()
Definition: tmop.hpp:1503
mfem::TMOPComboIntegrator::EnableLimiting
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Adds the limiting term to the first integrator. Disables it for the rest.
Definition: tmop.cpp:4165
mfem::TMOP_Metric_323::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:1153
mfem::InvariantsEvaluator2D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:335
mfem::TMOP_Metric_056::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:504
mfem::AnalyticAdaptTC::ComputeElementTargets
virtual void ComputeElementTargets(int e_id, const FiniteElement &fe, const IntegrationRule &ir, const Vector &elfun, DenseTensor &Jtr) const
Given an element and quadrature rule, computes ref-&gt;target transformation Jacobians for each quadratur...
Definition: tmop.cpp:1552
mfem::AnalyticAdaptTC::scalar_tspec
Coefficient * scalar_tspec
Definition: tmop.hpp:1310
mfem::TMOP_Metric_252::tau0
double & tau0
Definition: tmop.hpp:568
mfem::TMOP_Metric_301::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:713
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio
virtual void SetSerialDiscreteTargetAspectRatio(const GridFunction &tspec_)
Definition: tmop.cpp:1792
mfem::Geometry::GetGeomToPerfGeomJac
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
mfem::TMOP_Metric_311::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:868
mfem::TMOP_Metric_skew3D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:150
mfem::TMOP_Metric_303::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:822
mfem::TMOP_Metric_056::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:520
mfem::TMOP_Integrator::ComputeUntanglerMaxMuBarrier
double ComputeUntanglerMaxMuBarrier(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:4063
mfem::TMOP_Combo_QualityMetric::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:43
mfem::TMOPComboIntegrator::AddMultGradPA
virtual void AddMultGradPA(const Vector &, Vector &) const
Method for partially assembled gradient action.
Definition: tmop.cpp:4329
mfem::TargetConstructor::ComputeElementTargetsGradient
virtual void ComputeElementTargetsGradient(const IntegrationRule &ir, const Vector &elfun, IsoparametricTransformation &Tpr, DenseTensor &dJtr) const
Definition: tmop.cpp:1529
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::exponent
int exponent
Definition: tmop.hpp:147
mfem::FiniteElement::GetGeomType
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
mfem::Add
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2282
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
mfem::TMOP_QualityMetric
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
mfem::DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec
void FinalizeSerialDiscreteTargetSpec(const GridFunction &tspec_)
Definition: tmop.cpp:1812
mfem::Geometries
Geometry Geometries
Definition: fe.cpp:49
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
mfem::TMOP_Integrator::adapt_lim_coeff
Coefficient * adapt_lim_coeff
Definition: tmop.hpp:1602
mfem::TMOP_Metric_304::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:848
mfem::DenseMatrix::MultTranspose
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:201
mfem::AdaptivityEvaluator::~AdaptivityEvaluator
virtual ~AdaptivityEvaluator()
Definition: tmop.cpp:2654
mfem::InvariantsEvaluator2D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:181
mfem::CalcAdjugateTranspose
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2512
mfem::InvariantsEvaluator3D::Assemble_ddI3b
void Assemble_ddI3b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:1253
mfem::TMOP_Metric_315::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:710
mfem::FiniteElementSpace::GetElementVDofs
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
mfem::TMOPComboIntegrator::SetLimitingNodes
void SetLimitingNodes(const GridFunction &n0)
Update the original/reference nodes used for limiting.
Definition: tmop.cpp:4186
mfem::DiscreteAdaptTC::SetDiscreteTargetBase
void SetDiscreteTargetBase(const GridFunction &tspec_)
Definition: tmop.cpp:1735
mfem::TMOP_Integrator::adapt_lim_eval
AdaptivityEvaluator * adapt_lim_eval
Definition: tmop.hpp:1603
mfem::ParFiniteElementSpace::GetComm
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
mfem::f
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
mfem::TMOP_Metric_302::EvalWMatrixForm
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.cpp:757
mfem::TMOP_Metric_316::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:976
mfem::AddMult
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2408
mfem::TMOP_Integrator::discr_tc
DiscreteAdaptTC * discr_tc
Definition: tmop.hpp:1612
mfem::TMOP_Integrator::EnableLimiting
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Limiting of the mesh displacements (general version).
Definition: tmop.cpp:2690
mfem::TMOP_Metric_211::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:637
mfem::TMOP_Metric_352::tau0
double & tau0
Definition: tmop.hpp:923
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::GetWorstCaseType
virtual WorstCaseType GetWorstCaseType()
Definition: tmop.hpp:193
mfem::DiscreteAdaptTC::UpdateTargetSpecificationAtNode
void UpdateTargetSpecificationAtNode(const FiniteElement &el, ElementTransformation &T, int nodenum, int idir, const Vector &IntData)
Definition: tmop.cpp:1884
mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
mfem::DiscreteAdaptTC::SetTspecFromIntRule
void SetTspecFromIntRule(int e_id, const IntegrationRule &intrule)
Definition: tmop.cpp:1915
mfem::TMOP_Metric_050::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:462
mfem::ParGridFunction::Update
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:81
mfem::TMOP_Combo_QualityMetric::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:56
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
mfem::TMOP_Integrator::UpdateAfterMeshPositionChange
void UpdateAfterMeshPositionChange(const Vector &new_x, int new_x_ordering=Ordering::byNODES)
Definition: tmop.cpp:3957
mfem::TMOP_Metric_055::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:491
mfem::TMOP_Metric_009::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:336
mfem::DenseMatrix::Add
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:542
mfem::TMOP_Metric_301::EvalWMatrixForm
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.cpp:697
mfem::DiscreteAdaptTC::ResetRefinementTspecData
void ResetRefinementTspecData()
Definition: tmop.hpp:1538
mfem::TMOP_Metric_211::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:550
mfem::TMOP_Metric_360::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:1229
mfem::TMOP_Metric_303::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:807
mfem::TMOP_Integrator::AssembleElementVectorExact
void AssembleElementVectorExact(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Definition: tmop.cpp:3147
mfem::TMOPComboIntegrator::AssembleGradDiagonalPA
virtual void AssembleGradDiagonalPA(Vector &) const
Method for computing the diagonal of the gradient with partial assembly.
Definition: tmop.cpp:4313
mfem::TMOP_Integrator::targetC
const TargetConstructor * targetC
Definition: tmop.hpp:1573
mfem::TMOP_Metric_304::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:841
mfem::TMOP_Integrator::ElemPertEnergy
Array< Vector * > ElemPertEnergy
Definition: tmop.hpp:1626
mfem::TMOP_Integrator::lim_nodes0
const GridFunction * lim_nodes0
Definition: tmop.hpp:1587
mfem::TMOP_Metric_211::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:646
mfem::FiniteElement::CalcShape
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
mfem::TMOP_Integrator::adapt_lim_pgf0
const ParGridFunction * adapt_lim_pgf0
Definition: tmop.hpp:1599
mfem::TMOP_Metric_303::EvalWMatrixForm
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.cpp:800
mfem::InvariantsEvaluator2D::Get_dI2b
const scalar_t * Get_dI2b()
Definition: invariants.hpp:216
mfem::TMOP_Metric_058::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:561
mfem::TMOP_Metric_022::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:367
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::PMean
mfem::FiniteElement::GetNodes
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
mfem::InvariantsEvaluator3D::SetDerivativeMatrix
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 3, using column-major storage.
Definition: invariants.hpp:740
mfem::TMOPComboIntegrator::AssembleGradPA
virtual void AssembleGradPA(const Vector &, const FiniteElementSpace &)
Prepare the integrator for partial assembly (PA) gradient evaluations on the given FE space fes at th...
Definition: tmop.cpp:4304
mfem::TMOP_Metric_aspratio2D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:190
mfem::TMOP_Metric_323::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:789
mfem::DiscreteAdaptTC::SetParDiscreteTargetAspectRatio
virtual void SetParDiscreteTargetAspectRatio(const ParGridFunction &tspec_)
Definition: tmop.cpp:1709
mfem::DenseTensor::GetData
double * GetData(int k)
Definition: densemat.hpp:1058
mfem::TMOP_Metric_323::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1137
mfem::ElementTransformation::Weight
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
mfem::TMOP_Metric_077::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:574
mfem::DenseTensor::SetSize
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:1005
mfem::TMOPComboIntegrator::ParEnableNormalization
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:4280
mfem::TMOP_Metric_252::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:569
mfem::NCMesh::GetNumRootElements
int GetNumRootElements()
Return the number of root elements.
Definition: ncmesh.hpp:366
mfem::GridFunction::GetValues
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:513
mfem::TMOP_Integrator::AssembleElementVector
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: tmop.cpp:3118
mfem::TMOP_Metric_004::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:299
mfem::TMOPComboIntegrator::GetLocalStateEnergyPA
virtual double GetLocalStateEnergyPA(const Vector &) const
Compute the local (to the MPI rank) energy with partial assembly.
Definition: tmop.cpp:4337
mfem::TMOP_Metric_313::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:925
mfem::DiscreteAdaptTC::tspec_fesv
FiniteElementSpace * tspec_fesv
Definition: tmop.hpp:1375
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
mfem::TMOP_Metric_007::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:311
mfem::NCMesh
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
Definition: ncmesh.hpp:121
mfem::TMOP_Integrator::metric_coeff
Coefficient * metric_coeff
Definition: tmop.hpp:1580
mfem::TMOP_Metric_001::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:110
mfem::TMOP_Integrator::AssembleElemGradSurfFit
void AssembleElemGradSurfFit(const FiniteElement &el_x, IsoparametricTransformation &Tpr, DenseMatrix &mat)
Definition: tmop.cpp:3545
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::max_muT
double max_muT
Definition: tmop.hpp:146
mfem::Vector::SetData
void SetData(double *d)
Definition: vector.hpp:150
mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:668
mfem::TargetConstructor::ContainsVolumeInfo
virtual bool ContainsVolumeInfo() const
Checks if the target matrices contain non-trivial size specification.
Definition: tmop.cpp:1444
mfem::TMOP_Integrator::h_metric
TMOP_QualityMetric * h_metric
Definition: tmop.hpp:1571
mfem::TMOP_Combo_QualityMetric::tmop_q_arr
Array< TMOP_QualityMetric * > tmop_q_arr
Definition: tmop.hpp:85
mfem::Mesh::Dimension
int Dimension() const
Definition: mesh.hpp:1006
mfem::TMOP_Metric_321::EvalWMatrixForm
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.cpp:998
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::GetBarrierType
virtual BarrierType GetBarrierType()
Definition: tmop.hpp:191
mfem::DenseMatrix::GetColumn
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1292
mfem::InvariantsEvaluator2D::Assemble_ddI1b
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:257
mfem::Vector::HostRead
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:453
mfem::DiscreteAdaptTC::RestoreTargetSpecificationAtNode
void RestoreTargetSpecificationAtNode(ElementTransformation &T, int nodenum)
Definition: tmop.cpp:1901
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Pseudo
mfem::TMOP_Metric_007::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:303
mfem::TMOP_Metric_058::EvalWMatrixForm
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.cpp:534
mfem::FiniteElementSpace::GetElementDofs
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2678
mfem::DiscreteAdaptTC::SetTspecAtIndex
void SetTspecAtIndex(int idx, const GridFunction &tspec_)
Definition: tmop.cpp:1759
mfem::TMOP_Integrator::GetRefinementElementEnergy
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:2978
mfem::TargetConstructor::ComputeAvgVolume
void ComputeAvgVolume() const
Definition: tmop.cpp:1330
mfem::DenseTensor::SizeI
int SizeI() const
Definition: densemat.hpp:999
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2548
mfem::TMOP_Metric_321::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1008
mfem::TMOP_Integrator::~TMOP_Integrator
~TMOP_Integrator()
Definition: tmop.cpp:2678
omega
double omega
Definition: ex25.cpp:141
mfem::TMOP_Integrator::dxscale
double dxscale
Definition: tmop.hpp:1617
mfem::TMOP_Integrator::lim_coeff
Coefficient * lim_coeff
Definition: tmop.hpp:1588
tmop_tools.hpp
mfem::Vector::Min
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1202
mfem::TMOP_Combo_QualityMetric::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:33
mfem::TMOP_Metric_323::EvalWMatrixForm
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.cpp:1130
mfem::TMOP_Integrator::surf_fit_eval
AdaptivityEvaluator * surf_fit_eval
Definition: tmop.hpp:1609
mfem::DenseMatrix::Data
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::min_detT
double min_detT
Definition: tmop.hpp:145
p
double p(const Vector &x, double t)
Definition: navier_mms.cpp:53
mfem::TMOP_QualityMetric::SetTargetJacobian
virtual void SetTargetJacobian(const DenseMatrix &Jtr_)
Specify the reference-element -&gt; target-element Jacobian matrix for the point of interest.
Definition: tmop.hpp:43
mfem::DenseMatrix::Transpose
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1381
mfem::MultAtB
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
Definition: complex_densemat.cpp:328
mfem::TMOP_Integrator::AssembleElemGradAdaptLim
void AssembleElemGradAdaptLim(const FiniteElement &el, IsoparametricTransformation &Tpr, const IntegrationRule &ir, const Vector &weights, DenseMatrix &m)
Definition: tmop.cpp:3430
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::muT_ep
double muT_ep
Definition: tmop.hpp:150
mfem::DenseMatrix::Trace
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:425
mfem::AddMultABt
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2847
mfem::TMOP_Integrator::AssembleElementGrad
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: tmop.cpp:3132
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
mfem::FiniteElement::GetDof
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
mfem::InvariantsEvaluator3D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:951
mfem::TMOP_Metric_211::eps
const double eps
Definition: tmop.hpp:549
mfem::DiscreteAdaptTC::tspec_refine
DenseMatrix tspec_refine
Definition: tmop.hpp:1364
mfem::InvariantsEvaluator3D::Get_dI1
const scalar_t * Get_dI1()
Definition: invariants.hpp:763
mfem::TMOP_Metric_077::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:581
mfem::Coefficient
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
mfem::TMOP_Metric_352::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:924
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::WorstCaseType::Beta
mfem::AnalyticAdaptTC::matrix_tspec
TMOPMatrixCoefficient * matrix_tspec
Definition: tmop.hpp:1312
mfem::InvariantsEvaluator2D::Get_dI1
const scalar_t * Get_dI1()
Definition: invariants.hpp:204
mfem::TargetConstructor::ComputeElementTargets
virtual void ComputeElementTargets(int e_id, const FiniteElement &fe, const IntegrationRule &ir, const Vector &elfun, DenseTensor &Jtr) const
Given an element and quadrature rule, computes ref-&gt;target transformation Jacobians for each quadratur...
Definition: tmop.cpp:1458
mfem::TMOP_Integrator::ActionIntegrationRule
const IntegrationRule & ActionIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1746
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
mfem::TMOP_Metric_304::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:856
mfem::TMOP_Metric_050::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:385
mfem::MultAAt
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2680
mfem::Operator::NumCols
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:76
mfem::TMOP_Metric_skew2D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:126
mfem::TMOP_Integrator::ReleasePADeviceMemory
void ReleasePADeviceMemory(bool copy_to_host=true)
Definition: tmop.cpp:2664
mfem::TargetConstructor::UsesPhysicalCoordinates
bool UsesPhysicalCoordinates() const
Return true if the methods ComputeElementTargets(), ComputeAllElementTargets(), and ComputeElementTar...
Definition: tmop.hpp:1259
mfem::TMOP_Metric_322::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:1077
mfem::TMOP_Metric_002::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:262
mfem::DenseTensor::HostWrite
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:1099
mfem::TMOPComboIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: tmop.cpp:4222
mfem::TMOP_Metric_050::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:443
mfem::CalcInverseTranspose
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2615
mfem::TMOP_Metric_360::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1222
mfem::TMOP_Metric_055::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:475
mfem::InvariantsEvaluator2D::Get_dI1b
const scalar_t * Get_dI1b()
Definition: invariants.hpp:208
mfem::DiscreteAdaptTC::Jtrcomp
DenseTensor Jtrcomp
Definition: tmop.hpp:1371
mfem::Vector::SetDataAndSize
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:157
mfem::AdaptivityEvaluator::SetParMetaInfo
void SetParMetaInfo(const ParMesh &m, const ParFiniteElementSpace &f)
Parallel version of SetSerialMetaInfo.
Definition: tmop.cpp:2634
mfem::TMOP_Metric_252::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:659
mfem::GridFunction::Update
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:165
mfem::TMOP_Metric_323::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:1144
mfem::TMOP_Metric_058::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:545
mfem::DiscreteAdaptTC::tspec_gf
GridFunction * tspec_gf
Definition: tmop.hpp:1377
mfem::InvariantsEvaluator3D::Get_dI1b
const scalar_t * Get_dI1b()
Definition: invariants.hpp:767
mfem::TMOP_Metric_313::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:900
mfem::TMOP_Metric_014::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:370
mfem::TensorBasisElement::GetDofMap
const Array< int > & GetDofMap() const
Get an Array&lt;int&gt; that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:1182
mfem::TMOP_Metric_007::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:296
mfem::TMOP_Metric_009::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:353
mfem::TMOP_Metric_315::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:949
mfem::TMOP_Metric_058::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:553
mfem::TMOP_Integrator::UpdateAfterMeshTopologyChange
void UpdateAfterMeshTopologyChange()
Definition: tmop.cpp:2815
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::SetMaxMuT
virtual void SetMaxMuT(double max_muT_)
Definition: tmop.hpp:189
mfem::TMOP_Integrator::PMatI
DenseMatrix PMatI
Definition: tmop.hpp:1638
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::SetMinDetT
virtual void SetMinDetT(double min_detT_)
Definition: tmop.hpp:187
mfem::TMOP_Metric_316::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:728
mfem::TMOP_Metric_009::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:344
mfem::DenseTensor::SizeJ
int SizeJ() const
Definition: densemat.hpp:1000
a
double a
Definition: lissajous.cpp:41
mfem::TMOP_Integrator::ComputeFDh
void ComputeFDh(const Vector &x, const FiniteElementSpace &fes)
Determines the perturbation, h, for FD-based approximation.
Definition: tmop.cpp:3977
mfem::InvariantsEvaluator3D::Get_dI2b
const scalar_t * Get_dI2b()
Definition: invariants.hpp:775
mfem::TMOP_Metric_360::EvalWMatrixForm
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.cpp:1215
mfem::Vector::ReadWrite
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:465
mfem::TMOP_Integrator::ComputeNormalizationEnergies
void ComputeNormalizationEnergies(const GridFunction &x, double &metric_energy, double &lim_energy, double &surf_fit_gf_energy)
Definition: tmop.cpp:3852
mfem::TMOP_LimiterFunction::Eval
virtual double Eval(const Vector &x, const Vector &x0, double d) const =0
Returns the limiting function, f(x, x0, d).
mfem::TMOP_Integrator::surf_fit_coeff
Coefficient * surf_fit_coeff
Definition: tmop.hpp:1608
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::DiscreteAdaptTC::UpdateHessianTargetSpecification
void UpdateHessianTargetSpecification(const Vector &x, double dx, bool use_flag=false, int x_ordering=Ordering::byNODES)
Definition: tmop.cpp:2546
mfem::DiscreteAdaptTC::adapt_eval
AdaptivityEvaluator * adapt_eval
Definition: tmop.hpp:1394
mfem::TMOP_Metric_316::EvalWMatrixForm
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.cpp:962
mfem::DiscreteAdaptTC::ptspec_fesv
ParFiniteElementSpace * ptspec_fesv
Definition: tmop.hpp:1380
mfem::TMOP_Integrator::dx
double dx
Definition: tmop.hpp:1616
mfem::TMOP_Metric_301::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:588
mfem::InvariantsEvaluator3D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition: invariants.hpp:1297
mfem::TMOP_Integrator::fes
const FiniteElementSpace * fes
Definition: tmop.hpp:1675
mfem::TMOPComboIntegrator::tmopi
Array< TMOP_Integrator * > tmopi
Definition: tmop.hpp:1997
mfem::IsoparametricTransformation
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
mfem::InvariantsEvaluator3D::Get_dI2
const scalar_t * Get_dI2()
Definition: invariants.hpp:771
mfem::TMOP_Integrator::ComputeMinJac
void ComputeMinJac(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:3923
mfem::TMOP_Metric_302::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:784
mfem::Mesh::GetElementTransformation
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
mfem::TMOP_Metric_304::EvalWMatrixForm
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.cpp:834
mfem::TMOP_Integrator::EnergyIntegrationRule
const IntegrationRule & EnergyIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1737
mfem::TMOP_Integrator::GradientIntegrationRule
const IntegrationRule & GradientIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1751
mfem::TMOP_Metric_004::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:278
mfem::TMOP_Integrator::GetFDDerivative
double GetFDDerivative(const FiniteElement &el, ElementTransformation &T, Vector &elfun, const int nodenum, const int idir, const double baseenergy, bool update_stored)
Definition: tmop.cpp:3621
mfem::TMOP_Integrator::fdflag
bool fdflag
Definition: tmop.hpp:1615
mfem::TMOP_Metric_056::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:512
mfem::TMOP_Integrator::surf_fit_normal
double surf_fit_normal
Definition: tmop.hpp:1610
dim
int dim
Definition: ex24.cpp:53
mfem::TMOP_Integrator::GetElementEnergy
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Computes the integral of W(Jacobian(Trt)) over a target zone.
Definition: tmop.cpp:2841
mfem::TMOP_Integrator::AssembleElemVecAdaptLim
void AssembleElemVecAdaptLim(const FiniteElement &el, IsoparametricTransformation &Tpr, const IntegrationRule &ir, const Vector &weights, DenseMatrix &mat)
Definition: tmop.cpp:3394
mfem::TMOP_Metric_001::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:104
mfem::IsoparametricTransformation::SetFE
void SetFE(const FiniteElement *FE)
Set the element that will be used to compute the transformations.
Definition: eltrans.hpp:381
mfem::TMOP_Integrator::ParUpdateAfterMeshTopologyChange
void ParUpdateAfterMeshTopologyChange()
Definition: tmop.cpp:2828
mfem::TMOPComboIntegrator::GetDerefinementElementEnergy
virtual double GetDerefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Definition: tmop.cpp:4251
mfem::InvariantsEvaluator3D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:793
mfem::TMOP_LimiterFunction::Eval_d1
virtual void Eval_d1(const Vector &x, const Vector &x0, double dist, Vector &d1) const =0
Returns the gradient of the limiting function f(x, x0, d) with respect to x.
mfem::FiniteElementSpace::GetFE
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
mfem::FiniteElementSpace::GetUpdateOperator
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
Definition: fespace.hpp:897
mfem::TargetConstructor::ComputeAllElementTargets_Fallback
void ComputeAllElementTargets_Fallback(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const
Definition: tmop.cpp:1372
mfem::DiscreteAdaptTC::GetTspecPert2H
const Vector & GetTspecPert2H()
Definition: tmop.hpp:1502
mfem::DiscreteAdaptTC::GetTspecPert1H
const Vector & GetTspecPert1H()
Definition: tmop.hpp:1501
mfem::TMOP_Integrator::AssembleElementGradFD
void AssembleElementGradFD(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Definition: tmop.cpp:3709
mfem::TMOP_Metric_316::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:968
mfem::TMOP_Metric_002::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:278
mfem::infinity
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
mfem::TMOP_Metric_301::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:722
mfem::DiscreteAdaptTC::SetParDiscreteTargetSpec
virtual void SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
Definition: tmop.cpp:1729
mfem::TMOP_Metric_007::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:317
mfem::Mesh::ncmesh
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition: mesh.hpp:273
mfem::TMOP_Metric_aspratio3D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:210
mfem::TMOP_Metric_022::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:406
t
RefCoord t[3]
Definition: ncmesh_tables.hpp:158
mfem::Coefficient::Eval
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
mfem::TMOP_Integrator::adapt_lim_gf
GridFunction * adapt_lim_gf
Definition: tmop.hpp:1601
mfem::InvariantsEvaluator2D::Get_dI2
const scalar_t * Get_dI2()
Definition: invariants.hpp:212
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::BarrierType::Shifted
mfem::TMOP_Metric_316::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:984
mfem::TMOP_Metric_004::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:272
mfem::TMOP_Integrator::surf_fit_gf
GridFunction * surf_fit_gf
Definition: tmop.hpp:1606
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::btype
BarrierType btype
Definition: tmop.hpp:151
alpha
const double alpha
Definition: ex15.cpp:369
mfem::TMOP_Metric_301::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:706
mfem::DiscreteAdaptTC::FinalizeParDiscreteTargetSpec
void FinalizeParDiscreteTargetSpec(const ParGridFunction &tspec_)
Definition: tmop.cpp:1633
mfem::DiscreteAdaptTC::tspec_pgf
ParGridFunction * tspec_pgf
Definition: tmop.hpp:1382
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::FiniteElementSpace::FEColl
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
mfem::DiscreteAdaptTC::GetDiscreteTargetSpec
void GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
Get one of the discrete fields from tspec.
Definition: tmop.cpp:1833
mfem::DenseMatrix::Mult
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
mfem::TMOP_Integrator::DSh
DenseMatrix DSh
Definition: tmop.hpp:1638
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::alpha
double alpha
Definition: tmop.hpp:148
mfem::Mesh::GetNodes
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7908
mfem::TMOP_Integrator::surf_fit_marker
const Array< bool > * surf_fit_marker
Definition: tmop.hpp:1607
mfem::TMOP_Integrator::DS
DenseMatrix DS
Definition: tmop.hpp:1638
mfem::TMOP_Metric_322::EvalWMatrixForm
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.cpp:1056
mfem::TMOP_Metric_303::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:628
mfem::DiscreteAdaptTC::SetTspecDataForDerefinement
void SetTspecDataForDerefinement(FiniteElementSpace *fes)
Computes target specification data with respect to the coarse FE space.
Definition: tmop.cpp:1941
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:71
mfem::TMOP_Metric_098::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:621
mfem::DiscreteAdaptTC::UpdateTargetSpecification
void UpdateTargetSpecification(const Vector &new_x, bool use_flag=false, int new_x_ordering=Ordering::byNODES)
Definition: tmop.cpp:1864
mfem::TMOP_Metric_055::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:483
mfem::TMOP_Integrator::EnableSurfaceFitting
void EnableSurfaceFitting(const GridFunction &s0, const Array< bool > &smarker, Coefficient &coeff, AdaptivityEvaluator &ae)
Fitting of certain DOFs to the zero level set of a function.
Definition: tmop.cpp:2750
mfem::Mesh::DeleteGeometricFactors
void DeleteGeometricFactors()
Destroy all GeometricFactors stored by the Mesh.
Definition: mesh.cpp:882
mfem::TMOP_Metric_360::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:942
s
RefCoord s[3]
Definition: ncmesh_tables.hpp:158
mfem::TMOP_Metric_211::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:651
mfem::TMOP_Combo_QualityMetric::EvalWMatrixForm
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.cpp:23
mfem::TMOP_Metric_311::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:877
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSpec
virtual void SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
Definition: tmop.cpp:1858
mfem::DiscreteAdaptTC::SetParDiscreteTargetSkew
virtual void SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
Definition: tmop.cpp:1699
mfem::TMOP_Metric_050::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:454
mfem::TMOPComboIntegrator::GetRefinementElementEnergy
virtual double GetRefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, const IntegrationRule &irule)
Definition: tmop.cpp:4238
mfem::TMOP_Integrator::exact_action
bool exact_action
Definition: tmop.hpp:1623
mfem::DiscreteAdaptTC::SetParDiscreteTargetOrientation
virtual void SetParDiscreteTargetOrientation(const ParGridFunction &tspec_)
Definition: tmop.cpp:1719
mfem::DenseMatrix::UseExternalData
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:80
mfem::ParGridFunction
Class for parallel grid function.
Definition: pgridfunc.hpp:32
mfem::InvariantsEvaluator3D::Assemble_ddI1b
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:827
mfem::TargetConstructor
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:1166
mfem::TMOP_AMetric_036::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1293
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
mfem::DiscreteAdaptTC::~DiscreteAdaptTC
virtual ~DiscreteAdaptTC()
Definition: tmop.cpp:2613
mfem::Operator
Abstract operator.
Definition: operator.hpp:24
mfem::IsoparametricTransformation::GetPointMat
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:404
mfem::TMOP_Metric_002::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:256
mfem::FiniteElement::CalcDShape
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...
mfem::AnalyticAdaptTC::SetAnalyticTargetSpec
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, TMOPMatrixCoefficient *mspec)
Definition: tmop.cpp:1543
mfem::TMOP_Metric_004::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:284
mfem::DenseTensor
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:953
mfem::ParMesh
Class for parallel meshes.
Definition: pmesh.hpp:32
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSize
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
Definition: tmop.cpp:1772
mfem::DenseMatrix::FNorm2
double FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Definition: densemat.hpp:246
mfem::IntRules
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
mfem::TMOP_Integrator::lim_func
TMOP_LimiterFunction * lim_func
Definition: tmop.hpp:1592
mfem::TMOP_Metric_315::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:933
mfem::TMOP_Metric_303::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:814
mfem::TMOP_Metric_302::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:607
mfem::TMOP_Integrator::adapt_lim_gf0
const GridFunction * adapt_lim_gf0
Definition: tmop.hpp:1597
mfem::TMOPComboIntegrator::AssembleElementVector
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: tmop.cpp:4206
mfem::TMOP_Integrator::ComputeMinDetT
double ComputeMinDetT(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:4021
mfem::TMOP_Integrator::dim
int dim
Definition: tmop.hpp:1667
mfem::TMOP_Metric_322::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:1093
mfem::TMOP_Integrator::Jpt
DenseMatrix Jpt
Definition: tmop.hpp:1638
mfem::Vector::Read
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:449
mfem::TMOP_Combo_QualityMetric::wt_arr
Array< double > wt_arr
Definition: tmop.hpp:86
mfem::TMOP_Integrator::ParEnableNormalization
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:3839
mfem::TMOP_Metric_352::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:1177
mfem::TMOP_Integrator::P
DenseMatrix P
Definition: tmop.hpp:1638
mfem::TargetConstructor::GetComm
MPI_Comm GetComm() const
Definition: tmop.hpp:1236
mfem::ElementTransformation::mesh
class Mesh * mesh
The Mesh object containing the element.
Definition: eltrans.hpp:84
mfem::TMOP_Metric_352::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:1189
mfem::TMOP_Metric_360::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:1238
mfem::TMOP_WorstCaseUntangleOptimizer_Metric::wctype
WorstCaseType wctype
Definition: tmop.hpp:152
mfem::DiscreteAdaptTC::UpdateAfterMeshTopologyChange
void UpdateAfterMeshTopologyChange()
Update all discrete fields based on tspec and update for AMR.
Definition: tmop.cpp:1847
mfem::TMOP_Metric_313::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:920
mfem::TargetConstructor::Parallel
bool Parallel() const
Definition: tmop.hpp:1235
mfem::TMOP_Integrator::fd_call_flag
bool fd_call_flag
Definition: tmop.hpp:1620
mfem::TMOP_Metric_352::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Definition: tmop.cpp:1169
mfem::AdaptivityEvaluator::fes
FiniteElementSpace * fes
Definition: tmop.hpp:1117
mfem::TMOP_Metric_077::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:590
mfem::ParGridFunction::ParFESpace
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
mfem::TMOP_Metric_304::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:649
mfem::DiscreteAdaptTC::coarse_tspec_fesv
FiniteElementSpace * coarse_tspec_fesv
Definition: tmop.hpp:1376
mfem::InvariantsEvaluator3D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:737
Generated on Wed Dec 7 2022 10:41:46 for MFEM by   doxygen 1.8.5