MFEM  v4.3.0
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-2021, 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::EvalW(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]->EvalW(Jpt);
29  }
30  return metric;
31 }
32 
33 void TMOP_Combo_QualityMetric::EvalP(const DenseMatrix &Jpt,
34  DenseMatrix &P) const
35 {
36  DenseMatrix Pt(P.Size());
37  P = 0.0;
38  for (int i = 0; i < tmop_q_arr.Size(); i++)
39  {
40  tmop_q_arr[i]->EvalP(Jpt, Pt);
41  Pt *= wt_arr[i];
42  P += Pt;
43  }
44 }
45 
46 void TMOP_Combo_QualityMetric::AssembleH(const DenseMatrix &Jpt,
47  const DenseMatrix &DS,
48  const double weight,
49  DenseMatrix &A) const
50 {
51  DenseMatrix At(A.Size());
52  for (int i = 0; i < tmop_q_arr.Size(); i++)
53  {
54  At = 0.0;
55  tmop_q_arr[i]->AssembleH(Jpt, DS, weight, At);
56  At *= wt_arr[i];
57  A += At;
58  }
59 }
60 
61 double TMOP_Metric_001::EvalW(const DenseMatrix &Jpt) const
62 {
63  ie.SetJacobian(Jpt.GetData());
64  return ie.Get_I1();
65 }
66 
67 void TMOP_Metric_001::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
68 {
69  ie.SetJacobian(Jpt.GetData());
70  P = ie.Get_dI1();
71 }
72 
73 void TMOP_Metric_001::AssembleH(const DenseMatrix &Jpt,
74  const DenseMatrix &DS,
75  const double weight,
76  DenseMatrix &A) const
77 {
78  ie.SetJacobian(Jpt.GetData());
79  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
80  ie.Assemble_ddI1(weight, A.GetData());
81 }
82 
83 double TMOP_Metric_skew2D::EvalW(const DenseMatrix &Jpt) const
84 {
85  MFEM_VERIFY(Jtr != NULL,
86  "Requires a target Jacobian, use SetTargetJacobian().");
87 
88  DenseMatrix Jpr(2, 2);
89  Mult(Jpt, *Jtr, Jpr);
90 
91  Vector col1, col2;
92  Jpr.GetColumn(0, col1);
93  Jpr.GetColumn(1, col2);
94  double norm_prod = col1.Norml2() * col2.Norml2();
95  const double cos_Jpr = (col1 * col2) / norm_prod,
96  sin_Jpr = fabs(Jpr.Det()) / norm_prod;
97 
98  Jtr->GetColumn(0, col1);
99  Jtr->GetColumn(1, col2);
100  norm_prod = col1.Norml2() * col2.Norml2();
101  const double cos_Jtr = (col1 * col2) / norm_prod,
102  sin_Jtr = fabs(Jtr->Det()) / norm_prod;
103 
104  return 0.5 * (1.0 - cos_Jpr * cos_Jtr - sin_Jpr * sin_Jtr);
105 }
106 
107 double TMOP_Metric_skew3D::EvalW(const DenseMatrix &Jpt) const
108 {
109  MFEM_VERIFY(Jtr != NULL,
110  "Requires a target Jacobian, use SetTargetJacobian().");
111 
112  DenseMatrix Jpr(3, 3);
113  Mult(Jpt, *Jtr, Jpr);
114 
115  Vector col1, col2, col3;
116  Jpr.GetColumn(0, col1);
117  Jpr.GetColumn(1, col2);
118  Jpr.GetColumn(2, col3);
119  double norm_c1 = col1.Norml2(),
120  norm_c2 = col2.Norml2(),
121  norm_c3 = col3.Norml2();
122  double cos_Jpr_12 = (col1 * col2) / (norm_c1 * norm_c2),
123  cos_Jpr_13 = (col1 * col3) / (norm_c1 * norm_c3),
124  cos_Jpr_23 = (col2 * col3) / (norm_c2 * norm_c3);
125  double sin_Jpr_12 = std::sqrt(1.0 - cos_Jpr_12 * cos_Jpr_12),
126  sin_Jpr_13 = std::sqrt(1.0 - cos_Jpr_13 * cos_Jpr_13),
127  sin_Jpr_23 = std::sqrt(1.0 - cos_Jpr_23 * cos_Jpr_23);
128 
129  Jtr->GetColumn(0, col1);
130  Jtr->GetColumn(1, col2);
131  Jtr->GetColumn(2, col3);
132  norm_c1 = col1.Norml2();
133  norm_c2 = col2.Norml2(),
134  norm_c3 = col3.Norml2();
135  double cos_Jtr_12 = (col1 * col2) / (norm_c1 * norm_c2),
136  cos_Jtr_13 = (col1 * col3) / (norm_c1 * norm_c3),
137  cos_Jtr_23 = (col2 * col3) / (norm_c2 * norm_c3);
138  double sin_Jtr_12 = std::sqrt(1.0 - cos_Jtr_12 * cos_Jtr_12),
139  sin_Jtr_13 = std::sqrt(1.0 - cos_Jtr_13 * cos_Jtr_13),
140  sin_Jtr_23 = std::sqrt(1.0 - cos_Jtr_23 * cos_Jtr_23);
141 
142  return (3.0 - cos_Jpr_12 * cos_Jtr_12 - sin_Jpr_12 * sin_Jtr_12
143  - cos_Jpr_13 * cos_Jtr_13 - sin_Jpr_13 * sin_Jtr_13
144  - cos_Jpr_23 * cos_Jtr_23 - sin_Jpr_23 * sin_Jtr_23) / 6.0;
145 }
146 
147 double TMOP_Metric_aspratio2D::EvalW(const DenseMatrix &Jpt) const
148 {
149  MFEM_VERIFY(Jtr != NULL,
150  "Requires a target Jacobian, use SetTargetJacobian().");
151 
152  DenseMatrix Jpr(2, 2);
153  Mult(Jpt, *Jtr, Jpr);
154 
155  Vector col1, col2;
156  Jpr.GetColumn(0, col1);
157  Jpr.GetColumn(1, col2);
158  const double ratio_Jpr = col2.Norml2() / col1.Norml2();
159 
160  Jtr->GetColumn(0, col1);
161  Jtr->GetColumn(1, col2);
162  const double ratio_Jtr = col2.Norml2() / col1.Norml2();
163 
164  return 0.5 * (ratio_Jpr / ratio_Jtr + ratio_Jtr / ratio_Jpr) - 1.0;
165 }
166 
167 double TMOP_Metric_aspratio3D::EvalW(const DenseMatrix &Jpt) const
168 {
169  MFEM_VERIFY(Jtr != NULL,
170  "Requires a target Jacobian, use SetTargetJacobian().");
171 
172  DenseMatrix Jpr(3, 3);
173  Mult(Jpt, *Jtr, Jpr);
174 
175  Vector col1, col2, col3;
176  Jpr.GetColumn(0, col1);
177  Jpr.GetColumn(1, col2);
178  Jpr.GetColumn(2, col3);
179  double norm_c1 = col1.Norml2(),
180  norm_c2 = col2.Norml2(),
181  norm_c3 = col3.Norml2();
182  double ratio_Jpr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
183  ratio_Jpr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
184  ratio_Jpr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
185 
186  Jtr->GetColumn(0, col1);
187  Jtr->GetColumn(1, col2);
188  Jtr->GetColumn(2, col3);
189  norm_c1 = col1.Norml2();
190  norm_c2 = col2.Norml2();
191  norm_c3 = col3.Norml2();
192  double ratio_Jtr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
193  ratio_Jtr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
194  ratio_Jtr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
195 
196  return ( 0.5 * (ratio_Jpr_1 / ratio_Jtr_1 + ratio_Jtr_1 / ratio_Jpr_1) +
197  0.5 * (ratio_Jpr_2 / ratio_Jtr_2 + ratio_Jtr_2 / ratio_Jpr_2) +
198  0.5 * (ratio_Jpr_3 / ratio_Jtr_3 + ratio_Jtr_3 / ratio_Jpr_3) - 3.0
199  ) / 3.0;
200 }
201 
202 double TMOP_Metric_002::EvalW(const DenseMatrix &Jpt) const
203 {
204  ie.SetJacobian(Jpt.GetData());
205  return 0.5 * ie.Get_I1b() - 1.0;
206 }
207 
208 void TMOP_Metric_002::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
209 {
210  ie.SetJacobian(Jpt.GetData());
211  P.Set(0.5, ie.Get_dI1b());
212 }
213 
214 void TMOP_Metric_002::AssembleH(const DenseMatrix &Jpt,
215  const DenseMatrix &DS,
216  const double weight,
217  DenseMatrix &A) const
218 {
219  ie.SetJacobian(Jpt.GetData());
220  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
221  ie.Assemble_ddI1b(0.5*weight, A.GetData());
222 }
223 
224 double TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
225 {
226  // mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
227  ie.SetJacobian(Jpt.GetData());
228  return ie.Get_I1()*(1. + 1./ie.Get_I2()) - 4.0;
229 }
230 
231 void TMOP_Metric_007::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
232 {
233  // P = d(I1*(1 + 1/I2)) = (1 + 1/I2) dI1 - I1/I2^2 dI2
234  ie.SetJacobian(Jpt.GetData());
235  const double I2 = ie.Get_I2();
236  Add(1. + 1./I2, ie.Get_dI1(), -ie.Get_I1()/(I2*I2), ie.Get_dI2(), P);
237 }
238 
239 void TMOP_Metric_007::AssembleH(const DenseMatrix &Jpt,
240  const DenseMatrix &DS,
241  const double weight,
242  DenseMatrix &A) const
243 {
244  // P = d(I1*(1 + 1/I2))
245  // = (1 + 1/I2) dI1 - I1/I2^2 dI2
246  //
247  // dP = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 -
248  // (dI2 x d(I1/I2^2)) - I1/I2^2 ddI2
249  // = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 +
250  // (-1/I2^2) (dI2 x [dI1 - 2 I1/I2 dI2]) - I1/I2^2 ddI2
251  // = (-1/I2^2) (dI1 x dI2 + dI2 x dI1) + (1 + 1/I2) ddI1 +
252  // (2 I1/I2^3) (dI2 x dI2) - I1/I2^2 ddI2
253  ie.SetJacobian(Jpt.GetData());
254  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
255  const double c1 = 1./ie.Get_I2();
256  const double c2 = weight*c1*c1;
257  const double c3 = ie.Get_I1()*c2;
258  ie.Assemble_ddI1(weight*(1. + c1), A.GetData());
259  ie.Assemble_ddI2(-c3, A.GetData());
260  ie.Assemble_TProd(-c2, ie.Get_dI1(), ie.Get_dI2(), A.GetData());
261  ie.Assemble_TProd(2*c1*c3, ie.Get_dI2(), A.GetData());
262 }
263 
264 double TMOP_Metric_009::EvalW(const DenseMatrix &Jpt) const
265 {
266  // mu_9 = det(J)*|J-J^{-t}|^2 = I1b * (I2b^2 + 1) - 4 * I2b
267  // = (I1 - 4)*I2b + I1b
268  ie.SetJacobian(Jpt.GetData());
269  return (ie.Get_I1() - 4.0)*ie.Get_I2b() + ie.Get_I1b();
270 }
271 
272 void TMOP_Metric_009::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
273 {
274  // mu_9 = (I1 - 4)*I2b + I1b
275  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
276  ie.SetJacobian(Jpt.GetData());
277  Add(ie.Get_I1() - 4.0, ie.Get_dI2b(), ie.Get_I2b(), ie.Get_dI1(), P);
278  P += ie.Get_dI1b();
279 }
280 
281 void TMOP_Metric_009::AssembleH(const DenseMatrix &Jpt,
282  const DenseMatrix &DS,
283  const double weight,
284  DenseMatrix &A) const
285 {
286  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
287  // dP = dI2b x dI1 + (I1-4)*ddI2b + dI1 x dI2b + I2b*ddI1 + ddI1b
288  // = (dI1 x dI2b + dI2b x dI1) + (I1-4)*ddI2b + I2b*ddI1 + ddI1b
289  ie.SetJacobian(Jpt.GetData());
290  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
291  ie.Assemble_TProd(weight, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
292  ie.Assemble_ddI2b(weight*(ie.Get_I1()-4.0), A.GetData());
293  ie.Assemble_ddI1(weight*ie.Get_I2b(), A.GetData());
294  ie.Assemble_ddI1b(weight, A.GetData());
295 }
296 
297 // mu_14 = |T-I|^2
298 double TMOP_Metric_014::EvalW(const DenseMatrix &Jpt) const
299 {
300  MFEM_VERIFY(Jtr != NULL,
301  "Requires a target Jacobian, use SetTargetJacobian().");
302 
303  DenseMatrix Id(2,2);
304 
305  Id(0,0) = 1; Id(0,1) = 0;
306  Id(1,0) = 0; Id(1,1) = 1;
307 
308  DenseMatrix Mat(2,2);
309  Mat = Jpt;
310  Mat.Add(-1,Id);
311  return Mat.FNorm2();
312 }
313 
314 double TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const
315 {
316  // mu_22 = (0.5*|J|^2 - det(J)) / (det(J) - tau0)
317  // = (0.5*I1 - I2b) / (I2b - tau0)
318  ie.SetJacobian(Jpt.GetData());
319  const double I2b = ie.Get_I2b();
320 
321  double d = I2b - min_detT;
322  if (d < 0.0 && min_detT == 0.0)
323  {
324  // The mesh has been untangled, but it's still possible to get negative
325  // detJ in FD calculations, as they move the nodes around with some small
326  // increments and can produce negative determinants. Thus we put a small
327  // value in the denominator. Note that here I2b < 0.
328  d = - I2b * 0.1;
329  }
330 
331  return (0.5*ie.Get_I1() - I2b) / d;
332 }
333 
334 void TMOP_Metric_022::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
335 {
336  // mu_22 = (0.5*I1 - I2b) / (I2b - tau0)
337  // P = 1/(I2b - tau0)*(0.5*dI1 - dI2b) - (0.5*I1 - I2b)/(I2b - tau0)^2*dI2b
338  // = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
339  ie.SetJacobian(Jpt.GetData());
340  const double c1 = 1.0/(ie.Get_I2b() - min_detT);
341  Add(c1/2, ie.Get_dI1(), (min_detT - ie.Get_I1()/2)*c1*c1, ie.Get_dI2b(), P);
342 }
343 
344 void TMOP_Metric_022::AssembleH(const DenseMatrix &Jpt,
345  const DenseMatrix &DS,
346  const double weight,
347  DenseMatrix &A) const
348 {
349  // P = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
350  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b) + 0.5/(I2b - tau0)*ddI1
351  // + (dI2b x dz) + z*ddI2b
352  //
353  // z = (tau0 - 0.5*I1)/(I2b - tau0)^2
354  // dz = -0.5/(I2b - tau0)^2*dI1 - 2*(tau0 - 0.5*I1)/(I2b - tau0)^3*dI2b
355  //
356  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b + dI2b x dI1)
357  // -2*z/(I2b - tau0)*(dI2b x dI2b)
358  // +0.5/(I2b - tau0)*ddI1 + z*ddI2b
359  ie.SetJacobian(Jpt.GetData());
360  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
361  const double c1 = 1.0/(ie.Get_I2b() - min_detT);
362  const double c2 = weight*c1/2;
363  const double c3 = c1*c2;
364  const double c4 = (2*min_detT - ie.Get_I1())*c3; // weight*z
365  ie.Assemble_TProd(-c3, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
366  ie.Assemble_TProd(-2*c1*c4, ie.Get_dI2b(), A.GetData());
367  ie.Assemble_ddI1(c2, A.GetData());
368  ie.Assemble_ddI2b(c4, A.GetData());
369 }
370 
371 double TMOP_Metric_050::EvalW(const DenseMatrix &Jpt) const
372 {
373  // mu_50 = 0.5*|J^t J|^2/det(J)^2 - 1
374  // = 0.5*(l1^4 + l2^4)/(l1*l2)^2 - 1
375  // = 0.5*((l1/l2)^2 + (l2/l1)^2) - 1 = 0.5*(l1/l2 - l2/l1)^2
376  // = 0.5*(l1/l2 + l2/l1)^2 - 2 = 0.5*I1b^2 - 2
377  ie.SetJacobian(Jpt.GetData());
378  const double I1b = ie.Get_I1b();
379  return 0.5*I1b*I1b - 2.0;
380 }
381 
382 void TMOP_Metric_050::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
383 {
384  // mu_50 = 0.5*I1b^2 - 2
385  // P = I1b*dI1b
386  ie.SetJacobian(Jpt.GetData());
387  P.Set(ie.Get_I1b(), ie.Get_dI1b());
388 }
389 
390 void TMOP_Metric_050::AssembleH(const DenseMatrix &Jpt,
391  const DenseMatrix &DS,
392  const double weight,
393  DenseMatrix &A) const
394 {
395  // P = I1b*dI1b
396  // dP = dI1b x dI1b + I1b*ddI1b
397  ie.SetJacobian(Jpt.GetData());
398  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
399  ie.Assemble_TProd(weight, ie.Get_dI1b(), A.GetData());
400  ie.Assemble_ddI1b(weight*ie.Get_I1b(), A.GetData());
401 }
402 
403 double TMOP_Metric_055::EvalW(const DenseMatrix &Jpt) const
404 {
405  // mu_55 = (det(J) - 1)^2 = (I2b - 1)^2
406  ie.SetJacobian(Jpt.GetData());
407  const double c1 = ie.Get_I2b() - 1.0;
408  return c1*c1;
409 }
410 
411 void TMOP_Metric_055::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
412 {
413  // mu_55 = (I2b - 1)^2
414  // P = 2*(I2b - 1)*dI2b
415  ie.SetJacobian(Jpt.GetData());
416  P.Set(2*(ie.Get_I2b() - 1.0), ie.Get_dI2b());
417 }
418 
419 void TMOP_Metric_055::AssembleH(const DenseMatrix &Jpt,
420  const DenseMatrix &DS,
421  const double weight,
422  DenseMatrix &A) const
423 {
424  // P = 2*(I2b - 1)*dI2b
425  // dP = 2*(dI2b x dI2b) + 2*(I2b - 1)*ddI2b
426  ie.SetJacobian(Jpt.GetData());
427  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
428  ie.Assemble_TProd(2*weight, ie.Get_dI2b(), A.GetData());
429  ie.Assemble_ddI2b(2*weight*(ie.Get_I2b() - 1.0), A.GetData());
430 }
431 
432 double TMOP_Metric_056::EvalW(const DenseMatrix &Jpt) const
433 {
434  // mu_56 = 0.5*(I2b + 1/I2b) - 1
435  ie.SetJacobian(Jpt.GetData());
436  const double I2b = ie.Get_I2b();
437  return 0.5*(I2b + 1.0/I2b) - 1.0;
438 }
439 
440 void TMOP_Metric_056::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
441 {
442  // mu_56 = 0.5*(I2b + 1/I2b) - 1
443  // P = 0.5*(1 - 1/I2b^2)*dI2b
444  ie.SetJacobian(Jpt.GetData());
445  P.Set(0.5 - 0.5/ie.Get_I2(), ie.Get_dI2b());
446 }
447 
448 void TMOP_Metric_056::AssembleH(const DenseMatrix &Jpt,
449  const DenseMatrix &DS,
450  const double weight,
451  DenseMatrix &A) const
452 {
453  // P = 0.5*(1 - 1/I2b^2)*dI2b
454  // dP = (1/I2b^3)*(dI2b x dI2b) + (0.5 - 0.5/I2)*ddI2b
455  ie.SetJacobian(Jpt.GetData());
456  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
457  ie.Assemble_TProd(weight/(ie.Get_I2()*ie.Get_I2b()),
458  ie.Get_dI2b(), A.GetData());
459  ie.Assemble_ddI2b(weight*(0.5 - 0.5/ie.Get_I2()), A.GetData());
460 }
461 
462 double TMOP_Metric_058::EvalW(const DenseMatrix &Jpt) const
463 {
464  // mu_58 = I1b*(I1b - 2)
465  ie.SetJacobian(Jpt.GetData());
466  const double I1b = ie.Get_I1b();
467  return I1b*(I1b - 1.0);
468 }
469 
470 void TMOP_Metric_058::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
471 {
472  // mu_58 = I1b*(I1b - 2)
473  // P = (2*I1b - 2)*dI1b
474  ie.SetJacobian(Jpt.GetData());
475  P.Set(2*ie.Get_I1b() - 2.0, ie.Get_dI1b());
476 }
477 
478 void TMOP_Metric_058::AssembleH(const DenseMatrix &Jpt,
479  const DenseMatrix &DS,
480  const double weight,
481  DenseMatrix &A) const
482 {
483  // P = (2*I1b - 2)*dI1b
484  // dP = 2*(dI1b x dI1b) + (2*I1b - 2)*ddI1b
485  ie.SetJacobian(Jpt.GetData());
486  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
487  ie.Assemble_TProd(2*weight, ie.Get_dI1b(), A.GetData());
488  ie.Assemble_ddI1b(weight*(2*ie.Get_I1b() - 2.0), A.GetData());
489 }
490 
491 double TMOP_Metric_077::EvalW(const DenseMatrix &Jpt) const
492 {
493  ie.SetJacobian(Jpt.GetData());
494  const double I2b = ie.Get_I2b();
495  return 0.5*(I2b*I2b + 1./(I2b*I2b) - 2.);
496 }
497 
498 void TMOP_Metric_077::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
499 {
500  // Using I2b^2 = I2.
501  // dmu77_dJ = 1/2 (1 - 1/I2^2) dI2_dJ.
502  ie.SetJacobian(Jpt.GetData());
503  const double I2 = ie.Get_I2();
504  P.Set(0.5 * (1.0 - 1.0 / (I2 * I2)), ie.Get_dI2());
505 }
506 
507 void TMOP_Metric_077::AssembleH(const DenseMatrix &Jpt,
508  const DenseMatrix &DS,
509  const double weight,
510  DenseMatrix &A) const
511 {
512  ie.SetJacobian(Jpt.GetData());
513  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
514  const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
515  ie.Assemble_ddI2(weight*0.5*(1.0 - I2inv_sq), A.GetData());
516  ie.Assemble_TProd(weight * I2inv_sq / I2, ie.Get_dI2(), A.GetData());
517 }
518 
519 // mu_85 = |T-T'|^2, where T'= |T|*I/sqrt(2)
520 double TMOP_Metric_085::EvalW(const DenseMatrix &Jpt) const
521 {
522  MFEM_VERIFY(Jtr != NULL,
523  "Requires a target Jacobian, use SetTargetJacobian().");
524 
525  DenseMatrix Id(2,2);
526  DenseMatrix Mat(2,2);
527  Mat = Jpt;
528 
529  Id(0,0) = 1; Id(0,1) = 0;
530  Id(1,0) = 0; Id(1,1) = 1;
531  Id *= Mat.FNorm()/pow(2,0.5);
532 
533  Mat.Add(-1.,Id);
534  return Mat.FNorm2();
535 }
536 
537 // mu_98 = 1/(tau)|T-I|^2
538 double TMOP_Metric_098::EvalW(const DenseMatrix &Jpt) const
539 {
540  MFEM_VERIFY(Jtr != NULL,
541  "Requires a target Jacobian, use SetTargetJacobian().");
542 
543  DenseMatrix Id(2,2);
544 
545  Id(0,0) = 1; Id(0,1) = 0;
546  Id(1,0) = 0; Id(1,1) = 1;
547 
548  DenseMatrix Mat(2,2);
549  Mat = Jpt;
550  Mat.Add(-1,Id);
551  return Mat.FNorm2()/Jtr->Det();
552 }
553 
554 double TMOP_Metric_211::EvalW(const DenseMatrix &Jpt) const
555 {
556  // mu_211 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
557  // = (I2b - 1)^2 - I2b + sqrt(I2b^2 + eps)
558  ie.SetJacobian(Jpt.GetData());
559  const double I2b = ie.Get_I2b();
560  return (I2b - 1.0)*(I2b - 1.0) - I2b + std::sqrt(I2b*I2b + eps);
561 }
562 
563 void TMOP_Metric_211::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
564 {
565  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
566 }
567 
568 void TMOP_Metric_211::AssembleH(const DenseMatrix &Jpt,
569  const DenseMatrix &DS,
570  const double weight,
571  DenseMatrix &A) const
572 {
573  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
574 }
575 
576 double TMOP_Metric_252::EvalW(const DenseMatrix &Jpt) const
577 {
578  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0).
579  ie.SetJacobian(Jpt.GetData());
580  const double I2b = ie.Get_I2b();
581  return 0.5*(I2b - 1.0)*(I2b - 1.0)/(I2b - tau0);
582 }
583 
584 void TMOP_Metric_252::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
585 {
586  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
587  // P = (c - 0.5*c*c ) * dI2b
588  //
589  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
590  ie.SetJacobian(Jpt.GetData());
591  const double I2b = ie.Get_I2b();
592  const double c = (I2b - 1.0)/(I2b - tau0);
593  P.Set(c - 0.5*c*c, ie.Get_dI2b());
594 }
595 
596 void TMOP_Metric_252::AssembleH(const DenseMatrix &Jpt,
597  const DenseMatrix &DS,
598  const double weight,
599  DenseMatrix &A) const
600 {
601  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
602  //
603  // P = (c - 0.5*c*c ) * dI2b
604  // dP = (1 - c)^2/(I2b - tau0)*(dI2b x dI2b) + (c - 0.5*c*c)*ddI2b
605  ie.SetJacobian(Jpt.GetData());
606  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
607  const double I2b = ie.Get_I2b();
608  const double c0 = 1.0/(I2b - tau0);
609  const double c = c0*(I2b - 1.0);
610  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI2b(), A.GetData());
611  ie.Assemble_ddI2b(weight*(c - 0.5*c*c), A.GetData());
612 }
613 
614 double TMOP_Metric_301::EvalW(const DenseMatrix &Jpt) const
615 {
616  ie.SetJacobian(Jpt.GetData());
617  return std::sqrt(ie.Get_I1b()*ie.Get_I2b())/3. - 1.;
618 }
619 
620 void TMOP_Metric_301::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
621 {
622  // W = (1/3)*sqrt(I1b*I2b) - 1
623  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
624  ie.SetJacobian(Jpt.GetData());
625  const double a = 1./(6.*std::sqrt(ie.Get_I1b()*ie.Get_I2b()));
626  Add(a*ie.Get_I2b(), ie.Get_dI1b(), a*ie.Get_I1b(), ie.Get_dI2b(), P);
627 }
628 
629 void TMOP_Metric_301::AssembleH(const DenseMatrix &Jpt,
630  const DenseMatrix &DS,
631  const double weight,
632  DenseMatrix &A) const
633 {
634  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
635  // dW = (1/6)*[z2*dI1b + z1*dI2b], z1 = sqrt(I1b/I2b), z2 = sqrt(I2b/I1b)
636  // ddW = (1/6)*[dI1b x dz2 + z2*ddI1b + dI2b x dz1 + z1*ddI2b]
637  //
638  // dz1 = (1/2)*sqrt(I2b/I1b) [ (1/I2b)*dI1b + (I1b/(I2b*I2b))*dI2b ]
639  // = (1/2)/sqrt(I1b*I2b) [ dI1b + (I1b/I2b)*dI2b ]
640  // dz2 = (1/2)/sqrt(I1b*I2b) [ (I2b/I1b)*dI1b + dI2b ]
641  //
642  // dI1b x dz2 + dI2b x dz1 =
643  // (1/2)/sqrt(I1b*I2b) dI1b x [ (I2b/I1b)*dI1b + dI2b ] +
644  // (1/2)/sqrt(I1b*I2b) dI2b x [ dI1b + (I1b/I2b)*dI2b ] =
645  // (1/2)/sqrt(I1b*I2b) [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] x
646  // [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] =
647  // (1/2)/sqrt(I1b*I2b) [ 6*dW x 6*dW ] =
648  // (1/2)*(I1b*I2b)^{-3/2} (I2b*dI1b + I1b*dI2b) x (I2b*dI1b + I1b*dI2b)
649  //
650  // z1 = I1b/sqrt(I1b*I2b), z2 = I2b/sqrt(I1b*I2b)
651 
652  ie.SetJacobian(Jpt.GetData());
653  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
654  double d_I1b_I2b_data[9];
655  DenseMatrix d_I1b_I2b(d_I1b_I2b_data, 3, 3);
656  Add(ie.Get_I2b(), ie.Get_dI1b(), ie.Get_I1b(), ie.Get_dI2b(), d_I1b_I2b);
657  const double I1b_I2b = ie.Get_I1b()*ie.Get_I2b();
658  const double a = weight/(6*std::sqrt(I1b_I2b));
659  ie.Assemble_ddI1b(a*ie.Get_I2b(), A.GetData());
660  ie.Assemble_ddI2b(a*ie.Get_I1b(), A.GetData());
661  ie.Assemble_TProd(a/(2*I1b_I2b), d_I1b_I2b_data, A.GetData());
662 }
663 
664 double TMOP_Metric_302::EvalW(const DenseMatrix &Jpt) const
665 {
666  // mu_2 = |J|^2 |J^{-1}|^2 / 9 - 1
667  // = (l1^2 + l2^2 + l3^3)*(l1^{-2} + l2^{-2} + l3^{-2}) / 9 - 1
668  // = I1*(l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/l1^2/l2^2/l3^2/9 - 1
669  // = I1*I2/det(J)^2/9 - 1 = I1b*I2b/9-1
670  ie.SetJacobian(Jpt.GetData());
671  return ie.Get_I1b()*ie.Get_I2b()/9. - 1.;
672 }
673 
674 void TMOP_Metric_302::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
675 {
676  // mu_2 = I1b*I2b/9-1
677  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
678  ie.SetJacobian(Jpt.GetData());
679  Add(ie.Get_I1b()/9, ie.Get_dI2b(), ie.Get_I2b()/9, ie.Get_dI1b(), P);
680 }
681 
682 void TMOP_Metric_302::AssembleH(const DenseMatrix &Jpt,
683  const DenseMatrix &DS,
684  const double weight,
685  DenseMatrix &A) const
686 {
687  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
688  // dP = (dI2b x dI1b)/9 + (I1b/9)*ddI2b + (dI1b x dI2b)/9 + (I2b/9)*ddI1b
689  // = (dI2b x dI1b + dI1b x dI2b)/9 + (I1b/9)*ddI2b + (I2b/9)*ddI1b
690  ie.SetJacobian(Jpt.GetData());
691  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
692  const double c1 = weight/9;
693  ie.Assemble_TProd(c1, ie.Get_dI1b(), ie.Get_dI2b(), A.GetData());
694  ie.Assemble_ddI2b(c1*ie.Get_I1b(), A.GetData());
695  ie.Assemble_ddI1b(c1*ie.Get_I2b(), A.GetData());
696 }
697 
698 double TMOP_Metric_303::EvalW(const DenseMatrix &Jpt) const
699 {
700  ie.SetJacobian(Jpt.GetData());
701  return ie.Get_I1b()/3.0 - 1.0;
702 }
703 
704 void TMOP_Metric_303::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
705 {
706  ie.SetJacobian(Jpt.GetData());
707  P.Set(1./3., ie.Get_dI1b());
708 }
709 
710 void TMOP_Metric_303::AssembleH(const DenseMatrix &Jpt,
711  const DenseMatrix &DS,
712  const double weight,
713  DenseMatrix &A) const
714 {
715  ie.SetJacobian(Jpt.GetData());
716  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
717  ie.Assemble_ddI1b(weight/3., A.GetData());
718 }
719 
720 double TMOP_Metric_311::EvalW(const DenseMatrix &Jpt) const
721 {
722  // mu_311 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
723  // = (I3b - 1)^2 - I3b + sqrt(I3b^2 + eps)
724  ie.SetJacobian(Jpt.GetData());
725  const double I3b = ie.Get_I3b();
726  return (I3b - 1.0)*(I3b - 1.0) - I3b + std::sqrt(I3b*I3b + eps);
727 }
728 
729 void TMOP_Metric_311::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
730 {
731  ie.SetJacobian(Jpt.GetData());
732  const double I3b = ie.Get_I3b();
733  const double c = 2*I3b-3+(I3b)/(std::pow((I3b*I3b+eps),0.5));
734  P.Set(c, ie.Get_dI3b());
735 }
736 
737 void TMOP_Metric_311::AssembleH(const DenseMatrix &Jpt,
738  const DenseMatrix &DS,
739  const double weight,
740  DenseMatrix &A) const
741 {
742  ie.SetJacobian(Jpt.GetData());
743  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
744  const double I3b = ie.Get_I3b();
745  const double c0 = I3b*I3b+eps;
746  const double c1 = 2 + 1/(pow(c0,0.5)) - I3b*I3b/(pow(c0,1.5));
747  const double c2 = 2*I3b - 3 + I3b/(pow(c0,0.5));
748  ie.Assemble_TProd(weight*c1, ie.Get_dI3b(), A.GetData());
749  ie.Assemble_ddI3b(c2*weight, A.GetData());
750 }
751 
752 double TMOP_Metric_313::EvalW(const DenseMatrix &Jpt) const
753 {
754  ie.SetJacobian(Jpt.GetData());
755 
756  const double I3b = ie.Get_I3b();
757  double d = I3b - min_detT;
758  if (d < 0.0 && min_detT == 0.0)
759  {
760  // The mesh has been untangled, but it's still possible to get negative
761  // detJ in FD calculations, as they move the nodes around with some small
762  // increments and can produce negative determinants. Thus we put a small
763  // value in the denominator. Note that here I3b < 0.
764  d = - I3b * 0.1;
765  }
766 
767  const double c = std::pow(d, -2.0/3.0);
768 
769  return ie.Get_I1() * c / 3.0;
770 }
771 
772 void TMOP_Metric_313::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
773 {
774  MFEM_ABORT("Metric not implemented yet.");
775 }
776 
777 void TMOP_Metric_313::AssembleH(const DenseMatrix &Jpt,
778  const DenseMatrix &DS,
779  const double weight,
780  DenseMatrix &A) const
781 {
782  MFEM_ABORT("Metric not implemented yet.");
783 }
784 
785 double TMOP_Metric_315::EvalW(const DenseMatrix &Jpt) const
786 {
787  // mu_315 = mu_15_3D = (det(J) - 1)^2
788  ie.SetJacobian(Jpt.GetData());
789  const double c1 = ie.Get_I3b() - 1.0;
790  return c1*c1;
791 }
792 
793 void TMOP_Metric_315::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
794 {
795  // mu_315 = (I3b - 1)^2
796  // P = 2*(I3b - 1)*dI3b
797  ie.SetJacobian(Jpt.GetData());
798  P.Set(2*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
799 }
800 
801 void TMOP_Metric_315::AssembleH(const DenseMatrix &Jpt,
802  const DenseMatrix &DS,
803  const double weight,
804  DenseMatrix &A) const
805 {
806  // P = 2*(I3b - 1)*dI3b
807  // dP = 2*(dI3b x dI3b) + 2*(I3b - 1)*ddI3b
808  ie.SetJacobian(Jpt.GetData());
809  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
810  ie.Assemble_TProd(2*weight, ie.Get_dI3b(), A.GetData());
811  ie.Assemble_ddI3b(2*weight*(ie.Get_I3b() - 1.0), A.GetData());
812 }
813 
814 double TMOP_Metric_316::EvalW(const DenseMatrix &Jpt) const
815 {
816  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
817  ie.SetJacobian(Jpt.GetData());
818  const double I3b = ie.Get_I3b();
819  return 0.5*(I3b + 1.0/I3b) - 1.0;
820 }
821 
822 void TMOP_Metric_316::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
823 {
824  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
825  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
826  ie.SetJacobian(Jpt.GetData());
827  P.Set(0.5 - 0.5/ie.Get_I3(), ie.Get_dI3b());
828 }
829 
830 void TMOP_Metric_316::AssembleH(const DenseMatrix &Jpt,
831  const DenseMatrix &DS,
832  const double weight,
833  DenseMatrix &A) const
834 {
835  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
836  // dP = (1/I3b^3)*(dI3b x dI3b) + (0.5 - 0.5/I3)*ddI3b
837  ie.SetJacobian(Jpt.GetData());
838  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
839  ie.Assemble_TProd(weight/(ie.Get_I3()*ie.Get_I3b()),
840  ie.Get_dI3b(), A.GetData());
841  ie.Assemble_ddI3b(weight*(0.5 - 0.5/ie.Get_I3()), A.GetData());
842 }
843 
844 double TMOP_Metric_321::EvalW(const DenseMatrix &Jpt) const
845 {
846  // mu_321 = mu_21_3D = |J - J^{-t}|^2
847  // = |J|^2 + |J^{-1}|^2 - 6
848  // = |J|^2 + (l1^{-2} + l2^{-2} + l3^{-2}) - 6
849  // = |J|^2 + (l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/det(J)^2 - 6
850  // = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
851  ie.SetJacobian(Jpt.GetData());
852  return ie.Get_I1() + ie.Get_I2()/ie.Get_I3() - 6.0;
853 }
854 
855 void TMOP_Metric_321::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
856 {
857  // mu_321 = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
858  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
859  ie.SetJacobian(Jpt.GetData());
860  const double I3 = ie.Get_I3();
861  Add(1.0/I3, ie.Get_dI2(),
862  -2*ie.Get_I2()/(I3*ie.Get_I3b()), ie.Get_dI3b(), P);
863  P += ie.Get_dI1();
864 }
865 
866 void TMOP_Metric_321::AssembleH(const DenseMatrix &Jpt,
867  const DenseMatrix &DS,
868  const double weight,
869  DenseMatrix &A) const
870 {
871  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
872  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b) + (1/I3)*ddI2 + (dI3b x dz) + z*ddI3b
873  //
874  // z = -2*I2/I3b^3
875  // dz = (-2/I3b^3)*dI2 + (2*I2)*(3/I3b^4)*dI3b
876  //
877  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b + dI3b x dI2) + (1/I3)*ddI2
878  // + (6*I2/I3b^4)*(dI3b x dI3b) + (-2*I2/I3b^3)*ddI3b
879  ie.SetJacobian(Jpt.GetData());
880  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
881  const double c0 = 1.0/ie.Get_I3b();
882  const double c1 = weight*c0*c0;
883  const double c2 = -2*c0*c1;
884  const double c3 = c2*ie.Get_I2();
885  ie.Assemble_ddI1(weight, A.GetData());
886  ie.Assemble_ddI2(c1, A.GetData());
887  ie.Assemble_ddI3b(c3, A.GetData());
888  ie.Assemble_TProd(c2, ie.Get_dI2(), ie.Get_dI3b(), A.GetData());
889  ie.Assemble_TProd(-3*c0*c3, ie.Get_dI3b(), A.GetData());
890 }
891 
892 double TMOP_Metric_352::EvalW(const DenseMatrix &Jpt) const
893 {
894  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
895  ie.SetJacobian(Jpt.GetData());
896  const double I3b = ie.Get_I3b();
897  return 0.5*(I3b - 1.0)*(I3b - 1.0)/(I3b - tau0);
898 }
899 
900 void TMOP_Metric_352::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
901 {
902  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
903  // P = (I3b - 1)/(I3b - tau0)*dI3b + 0.5*(I3b - 1)^2*(-1/(I3b - tau0)^2)*dI3b
904  // = [ (I3b - 1)/(I3b - tau0) - 0.5*(I3b - 1)^2/(I3b - tau0)^2 ] * dI3b
905  // = (c - 0.5*c*c) * dI3b
906  ie.SetJacobian(Jpt.GetData());
907  const double I3b = ie.Get_I3b();
908  const double c = (I3b - 1.0)/(I3b - tau0);
909  P.Set(c - 0.5*c*c, ie.Get_dI3b());
910 }
911 
912 void TMOP_Metric_352::AssembleH(const DenseMatrix &Jpt,
913  const DenseMatrix &DS,
914  const double weight,
915  DenseMatrix &A) const
916 {
917  // c = (I3b - 1)/(I3b - tau0)
918  //
919  // P = (c - 0.5*c*c) * dI3b
920  // dP = (1 - c)*(dI3b x dc) + (c - 0.5*c*c)*ddI3b
921  //
922  // dc = 1/(I3b - tau0)*dI3b - (I3b - 1)/(I3b - tau)^2*dI3b =
923  // = (1 - c)/(I3b - tau0)*dI3b
924  //
925  // dP = (1 - c)^2/(I3b - tau0)*(dI3b x dI3b) + (c - 0.5*c*c)*ddI3b
926  ie.SetJacobian(Jpt.GetData());
927  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
928  const double I3b = ie.Get_I3b();
929  const double c0 = 1.0/(I3b - tau0);
930  const double c = c0*(I3b - 1.0);
931  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI3b(), A.GetData());
932  ie.Assemble_ddI3b(weight*(c - 0.5*c*c), A.GetData());
933 }
934 
935 double TMOP_AMetric_011::EvalW(const DenseMatrix &Jpt) const
936 {
937  MFEM_VERIFY(Jtr != NULL,
938  "Requires a target Jacobian, use SetTargetJacobian().");
939 
940  int dim = Jpt.Size();
941 
942  DenseMatrix Jpr(dim, dim);
943  Mult(Jpt, *Jtr, Jpr);
944 
945  double alpha = Jpr.Det(),
946  omega = Jtr->Det();
947 
948  DenseMatrix AdjAt(dim), WtW(dim), WRK(dim), Jtrt(dim);
949  CalcAdjugateTranspose(Jpr, AdjAt);
950  Jtrt.Transpose(*Jtr);
951  MultAAt(Jtrt, WtW);
952  WtW *= 1./omega;
953  Mult(AdjAt, WtW, WRK);
954 
955  WRK -= Jpr;
956  WRK *= -1.;
957 
958  return (0.25/alpha)*WRK.FNorm2();
959 }
960 
961 double TMOP_AMetric_014a::EvalW(const DenseMatrix &Jpt) const
962 {
963  MFEM_VERIFY(Jtr != NULL,
964  "Requires a target Jacobian, use SetTargetJacobian().");
965 
966  int dim = Jpt.Size();
967 
968  DenseMatrix Jpr(dim, dim);
969  Mult(Jpt, *Jtr, Jpr);
970 
971  double sqalpha = pow(Jpr.Det(), 0.5),
972  sqomega = pow(Jtr->Det(), 0.5);
973 
974  return 0.5*pow(sqalpha/sqomega - sqomega/sqalpha, 2.);
975 }
976 
977 double TMOP_AMetric_036::EvalW(const DenseMatrix &Jpt) const
978 {
979  MFEM_VERIFY(Jtr != NULL,
980  "Requires a target Jacobian, use SetTargetJacobian().");
981 
982  int dim = Jpt.Size();
983 
984  DenseMatrix Jpr(dim, dim);
985  Mult(Jpt, *Jtr, Jpr); // T*W = A
986 
987  double alpha = Jpr.Det(); // det(A)
988  Jpr -= *Jtr; // A-W
989 
990  return (1./alpha)*(Jpr.FNorm2()); //(1/alpha)*(|A-W|^2)
991 }
992 
993 double TMOP_AMetric_107a::EvalW(const DenseMatrix &Jpt) const
994 {
995  MFEM_VERIFY(Jtr != NULL,
996  "Requires a target Jacobian, use SetTargetJacobian().");
997 
998  int dim = Jpt.Size();
999 
1000  DenseMatrix Jpr(dim, dim);
1001  Mult(Jpt, *Jtr, Jpr);
1002 
1003  double alpha = Jpr.Det(),
1004  aw = Jpr.FNorm()/Jtr->FNorm();
1005 
1006  DenseMatrix W = *Jtr;
1007  W *= aw;
1008  Jpr -= W;
1009 
1010  return (0.5/alpha)*Jpr.FNorm2();
1011 }
1012 
1013 
1014 void TargetConstructor::ComputeAvgVolume() const
1015 {
1016  MFEM_VERIFY(nodes, "Nodes are not given!");
1017  MFEM_ASSERT(avg_volume == 0.0, "The average volume is already computed!");
1018 
1019  Mesh *mesh = nodes->FESpace()->GetMesh();
1020  const int NE = mesh->GetNE();
1021  IsoparametricTransformation Tr;
1022  double volume = 0.0;
1023 
1024  for (int i = 0; i < NE; i++)
1025  {
1026  mesh->GetElementTransformation(i, *nodes, &Tr);
1027  const IntegrationRule &ir =
1028  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr.OrderJ());
1029  for (int j = 0; j < ir.GetNPoints(); j++)
1030  {
1031  const IntegrationPoint &ip = ir.IntPoint(j);
1032  Tr.SetIntPoint(&ip);
1033  volume += ip.weight * Tr.Weight();
1034  }
1035  }
1036 
1037  NCMesh *ncmesh = mesh->ncmesh;
1038  if (Parallel() == false)
1039  {
1040  avg_volume = (ncmesh == NULL) ?
1041  volume / NE : volume / ncmesh->GetNumRootElements();
1042 
1043  }
1044 #ifdef MFEM_USE_MPI
1045  else
1046  {
1047  double area_NE[4];
1048  area_NE[0] = volume; area_NE[1] = NE;
1049  MPI_Allreduce(area_NE, area_NE + 2, 2, MPI_DOUBLE, MPI_SUM, comm);
1050  avg_volume = (ncmesh == NULL) ?
1051  area_NE[2] / area_NE[3] : area_NE[2] / ncmesh->GetNumRootElements();
1052  }
1053 #endif
1054 }
1055 
1056 void TargetConstructor::ComputeAllElementTargets_Fallback(
1057  const FiniteElementSpace &fes,
1058  const IntegrationRule &ir,
1059  const Vector &xe,
1060  DenseTensor &Jtr) const
1061 {
1062  // Fallback to the 1-element method, ComputeElementTargets()
1063 
1064  // When UsesPhysicalCoordinates() == true, we assume 'xe' uses
1065  // ElementDofOrdering::LEXICOGRAPHIC iff 'fe' is a TensorFiniteElement.
1066 
1067  const Mesh *mesh = fes.GetMesh();
1068  const int NE = mesh->GetNE();
1069  // Quick return for empty processors:
1070  if (NE == 0) { return; }
1071  const int dim = mesh->Dimension();
1072  MFEM_VERIFY(mesh->GetNumGeometries(dim) <= 1,
1073  "mixed meshes are not supported");
1074  MFEM_VERIFY(!fes.IsVariableOrder(), "variable orders are not supported");
1075  const FiniteElement &fe = *fes.GetFE(0);
1076  const int sdim = fes.GetVDim();
1077  const int nvdofs = sdim*fe.GetDof();
1078  MFEM_VERIFY(!UsesPhysicalCoordinates() ||
1079  xe.Size() == NE*nvdofs, "invalid input Vector 'xe'!");
1080  const int NQ = ir.GetNPoints();
1081  const Array<int> *dof_map = nullptr;
1082  if (UsesPhysicalCoordinates())
1083  {
1084  const TensorBasisElement *tfe =
1085  dynamic_cast<const TensorBasisElement *>(&fe);
1086  if (tfe)
1087  {
1088  dof_map = &tfe->GetDofMap();
1089  if (dof_map->Size() == 0) { dof_map = nullptr; }
1090  }
1091  }
1092 
1093  Vector elfun_lex, elfun_nat;
1094  DenseTensor J;
1095  xe.HostRead();
1096  Jtr.HostWrite();
1097  if (UsesPhysicalCoordinates() && dof_map != nullptr)
1098  {
1099  elfun_nat.SetSize(nvdofs);
1100  }
1101  for (int e = 0; e < NE; e++)
1102  {
1103  if (UsesPhysicalCoordinates())
1104  {
1105  if (!dof_map)
1106  {
1107  elfun_nat.SetDataAndSize(xe.GetData()+e*nvdofs, nvdofs);
1108  }
1109  else
1110  {
1111  elfun_lex.SetDataAndSize(xe.GetData()+e*nvdofs, nvdofs);
1112  const int ndofs = fe.GetDof();
1113  for (int d = 0; d < sdim; d++)
1114  {
1115  for (int i_lex = 0; i_lex < ndofs; i_lex++)
1116  {
1117  elfun_nat[(*dof_map)[i_lex]+d*ndofs] =
1118  elfun_lex[i_lex+d*ndofs];
1119  }
1120  }
1121  }
1122  }
1123  J.UseExternalData(Jtr(e*NQ).Data(), sdim, dim, NQ);
1124  ComputeElementTargets(e, fe, ir, elfun_nat, J);
1125  }
1126 }
1127 
1128 bool TargetConstructor::ContainsVolumeInfo() const
1129 {
1130  switch (target_type)
1131  {
1132  case IDEAL_SHAPE_UNIT_SIZE: return false;
1133  case IDEAL_SHAPE_EQUAL_SIZE:
1134  case IDEAL_SHAPE_GIVEN_SIZE:
1135  case GIVEN_SHAPE_AND_SIZE:
1136  case GIVEN_FULL: return true;
1137  default: MFEM_ABORT("TargetType not added to ContainsVolumeInfo.");
1138  }
1139  return false;
1140 }
1141 
1142 void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
1143  const IntegrationRule &ir,
1144  const Vector &elfun,
1145  DenseTensor &Jtr) const
1146 {
1147  MFEM_CONTRACT_VAR(elfun);
1148  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1149 
1150  const FiniteElement *nfe = (target_type != IDEAL_SHAPE_UNIT_SIZE) ?
1151  nodes->FESpace()->GetFE(e_id) : NULL;
1152  const DenseMatrix &Wideal =
1153  Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
1154  MFEM_ASSERT(Wideal.Height() == Jtr.SizeI(), "");
1155  MFEM_ASSERT(Wideal.Width() == Jtr.SizeJ(), "");
1156 
1157  switch (target_type)
1158  {
1159  case IDEAL_SHAPE_UNIT_SIZE:
1160  {
1161  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = Wideal; }
1162  break;
1163  }
1164  case IDEAL_SHAPE_EQUAL_SIZE:
1165  {
1166  if (avg_volume == 0.0) { ComputeAvgVolume(); }
1167  DenseMatrix W(Wideal.Height());
1168 
1169  NCMesh *ncmesh = nodes->FESpace()->GetMesh()->ncmesh;
1170  double el_volume = avg_volume;
1171  if (ncmesh)
1172  {
1173  el_volume = avg_volume / ncmesh->GetElementSizeReduction(e_id);
1174  }
1175 
1176  W.Set(std::pow(volume_scale * el_volume / Wideal.Det(),
1177  1./W.Height()), Wideal);
1178  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = W; }
1179  break;
1180  }
1181  case IDEAL_SHAPE_GIVEN_SIZE:
1182  case GIVEN_SHAPE_AND_SIZE:
1183  {
1184  const int dim = nfe->GetDim(), dof = nfe->GetDof();
1185  MFEM_ASSERT(dim == nodes->FESpace()->GetVDim(), "");
1186  DenseMatrix dshape(dof, dim), pos(dof, dim);
1187  Array<int> xdofs(dof * dim);
1188  Vector posV(pos.Data(), dof * dim);
1189  double detW;
1190 
1191  // always initialize detW to suppress a warning:
1192  detW = (target_type == IDEAL_SHAPE_GIVEN_SIZE) ? Wideal.Det() : 0.0;
1193  nodes->FESpace()->GetElementVDofs(e_id, xdofs);
1194  nodes->GetSubVector(xdofs, posV);
1195  for (int i = 0; i < ir.GetNPoints(); i++)
1196  {
1197  nfe->CalcDShape(ir.IntPoint(i), dshape);
1198  MultAtB(pos, dshape, Jtr(i));
1199  if (target_type == IDEAL_SHAPE_GIVEN_SIZE)
1200  {
1201  const double det = Jtr(i).Det();
1202  MFEM_VERIFY(det > 0.0, "The given mesh is inverted!");
1203  Jtr(i).Set(std::pow(det / detW, 1./dim), Wideal);
1204  }
1205  }
1206  break;
1207  }
1208  default:
1209  MFEM_ABORT("invalid target type!");
1210  }
1211 }
1212 
1213 void TargetConstructor::ComputeElementTargetsGradient(
1214  const IntegrationRule &ir,
1215  const Vector &elfun,
1216  IsoparametricTransformation &Tpr,
1217  DenseTensor &dJtr) const
1218 {
1219  MFEM_CONTRACT_VAR(elfun);
1220  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1221 
1222  // TODO: Compute derivative for targets with GIVEN_SHAPE or/and GIVEN_SIZE
1223  for (int i = 0; i < Tpr.GetFE()->GetDim()*ir.GetNPoints(); i++)
1224  { dJtr(i) = 0.; }
1225 }
1226 
1227 void AnalyticAdaptTC::SetAnalyticTargetSpec(Coefficient *sspec,
1228  VectorCoefficient *vspec,
1229  TMOPMatrixCoefficient *mspec)
1230 {
1231  scalar_tspec = sspec;
1232  vector_tspec = vspec;
1233  matrix_tspec = mspec;
1234 }
1235 
1236 void AnalyticAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
1237  const IntegrationRule &ir,
1238  const Vector &elfun,
1239  DenseTensor &Jtr) const
1240 {
1241  DenseMatrix point_mat;
1242  point_mat.UseExternalData(elfun.GetData(), fe.GetDof(), fe.GetDim());
1243 
1244  switch (target_type)
1245  {
1246  case GIVEN_FULL:
1247  {
1248  MFEM_VERIFY(matrix_tspec != NULL,
1249  "Target type GIVEN_FULL requires a MatrixCoefficient.");
1250 
1251  IsoparametricTransformation Tpr;
1252  Tpr.SetFE(&fe);
1253  Tpr.ElementNo = e_id;
1254  Tpr.ElementType = ElementTransformation::ELEMENT;
1255  Tpr.GetPointMat().Transpose(point_mat);
1256 
1257  for (int i = 0; i < ir.GetNPoints(); i++)
1258  {
1259  const IntegrationPoint &ip = ir.IntPoint(i);
1260  Tpr.SetIntPoint(&ip);
1261  matrix_tspec->Eval(Jtr(i), Tpr, ip);
1262  }
1263  break;
1264  }
1265  default:
1266  MFEM_ABORT("Incompatible target type for analytic adaptation!");
1267  }
1268 }
1269 
1270 void AnalyticAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
1271  const Vector &elfun,
1272  IsoparametricTransformation &Tpr,
1273  DenseTensor &dJtr) const
1274 {
1275  const FiniteElement *fe = Tpr.GetFE();
1276  DenseMatrix point_mat;
1277  point_mat.UseExternalData(elfun.GetData(), fe->GetDof(), fe->GetDim());
1278 
1279  switch (target_type)
1280  {
1281  case GIVEN_FULL:
1282  {
1283  MFEM_VERIFY(matrix_tspec != NULL,
1284  "Target type GIVEN_FULL requires a TMOPMatrixCoefficient.");
1285 
1286  for (int d = 0; d < fe->GetDim(); d++)
1287  {
1288  for (int i = 0; i < ir.GetNPoints(); i++)
1289  {
1290  const IntegrationPoint &ip = ir.IntPoint(i);
1291  Tpr.SetIntPoint(&ip);
1292  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1293  matrix_tspec->EvalGrad(dJtr_i, Tpr, ip, d);
1294  }
1295  }
1296  break;
1297  }
1298  default:
1299  MFEM_ABORT("Incompatible target type for analytic adaptation!");
1300  }
1301 }
1302 
1303 namespace internal
1304 {
1305 
1306 // MFEM_FORALL-based copy kernel -- used by protected methods below.
1307 // Needed as a workaround for the nvcc restriction that methods with MFEM_FORALL
1308 // in them must to be public.
1309 static inline void device_copy(double *d_dest, const double *d_src, int size)
1310 {
1311  MFEM_FORALL(i, size, d_dest[i] = d_src[i];);
1312 }
1313 
1314 } // namespace internal
1315 
1316 #ifdef MFEM_USE_MPI
1317 void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction
1318  &tspec_)
1319 {
1320  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1321  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1322 
1323  ParFiniteElementSpace *ptspec_fes = tspec_.ParFESpace();
1324 
1325  adapt_eval->SetParMetaInfo(*ptspec_fes->GetParMesh(),
1326  *ptspec_fes->FEColl(), ncomp);
1327  adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
1328 
1329  tspec_sav = tspec;
1330 
1331  delete tspec_fesv;
1332  tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
1333  tspec_fes->FEColl(), ncomp);
1334 }
1335 
1336 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const ParGridFunction &tspec_)
1337 {
1338  const int vdim = tspec_.FESpace()->GetVDim(),
1339  dof_cnt = tspec_.Size()/vdim;
1340  const auto tspec__d = tspec_.Read();
1341  auto tspec_d = tspec.ReadWrite();
1342  const int offset = idx*dof_cnt;
1343  internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
1344  FinalizeParDiscreteTargetSpec(tspec_);
1345 }
1346 
1347 void DiscreteAdaptTC::SetParDiscreteTargetSize(const ParGridFunction &tspec_)
1348 {
1349  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1350  sizeidx = ncomp;
1351  SetDiscreteTargetBase(tspec_);
1352  FinalizeParDiscreteTargetSpec(tspec_);
1353 }
1354 
1355 void DiscreteAdaptTC::SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
1356 {
1357  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1358  skewidx = ncomp;
1359  SetDiscreteTargetBase(tspec_);
1360  FinalizeParDiscreteTargetSpec(tspec_);
1361 }
1362 
1363 void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction
1364  &tspec_)
1365 {
1366  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, tspec_); return; }
1367  aspectratioidx = ncomp;
1368  SetDiscreteTargetBase(tspec_);
1369  FinalizeParDiscreteTargetSpec(tspec_);
1370 }
1371 
1372 void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction
1373  &tspec_)
1374 {
1375  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, tspec_); return; }
1376  orientationidx = ncomp;
1377  SetDiscreteTargetBase(tspec_);
1378  FinalizeParDiscreteTargetSpec(tspec_);
1379 }
1380 
1381 void DiscreteAdaptTC::SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
1382 {
1383  SetParDiscreteTargetSize(tspec_);
1384  FinalizeParDiscreteTargetSpec(tspec_);
1385 }
1386 #endif // MFEM_USE_MPI
1387 
1388 void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
1389 {
1390  const int vdim = tspec_.FESpace()->GetVDim(),
1391  dof_cnt = tspec_.Size()/vdim;
1392 
1393  ncomp += vdim;
1394 
1395  delete tspec_fes;
1396  tspec_fes = new FiniteElementSpace(tspec_.FESpace()->GetMesh(),
1397  tspec_.FESpace()->FEColl(), 1);
1398 
1399  // need to append data to tspec
1400  // make a copy of tspec->tspec_temp, increase its size, and
1401  // copy data from tspec_temp -> tspec, then add new entries
1402  Vector tspec_temp = tspec;
1403  tspec.UseDevice(true);
1404  tspec_sav.UseDevice(true);
1405  tspec.SetSize(ncomp*dof_cnt);
1406 
1407  const auto tspec_temp_d = tspec_temp.Read();
1408  auto tspec_d = tspec.ReadWrite();
1409  internal::device_copy(tspec_d, tspec_temp_d, tspec_temp.Size());
1410 
1411  const auto tspec__d = tspec_.Read();
1412  const int offset = (ncomp-vdim)*dof_cnt;
1413  internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
1414 }
1415 
1416 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
1417 {
1418  const int vdim = tspec_.FESpace()->GetVDim(),
1419  dof_cnt = tspec_.Size()/vdim;
1420 
1421  const auto tspec__d = tspec_.Read();
1422  auto tspec_d = tspec.ReadWrite();
1423  const int offset = idx*dof_cnt;
1424  internal::device_copy(tspec_d + offset, tspec__d, dof_cnt*vdim);
1425  FinalizeSerialDiscreteTargetSpec();
1426 }
1427 
1428 void DiscreteAdaptTC::SetSerialDiscreteTargetSize(const GridFunction &tspec_)
1429 {
1430 
1431  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1432  sizeidx = ncomp;
1433  SetDiscreteTargetBase(tspec_);
1434  FinalizeSerialDiscreteTargetSpec();
1435 }
1436 
1437 void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
1438 {
1439  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1440  skewidx = ncomp;
1441  SetDiscreteTargetBase(tspec_);
1442  FinalizeSerialDiscreteTargetSpec();
1443 }
1444 
1445 void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(
1446  const GridFunction &tspec_)
1447 {
1448  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, tspec_); return; }
1449  aspectratioidx = ncomp;
1450  SetDiscreteTargetBase(tspec_);
1451  FinalizeSerialDiscreteTargetSpec();
1452 }
1453 
1454 void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(
1455  const GridFunction &tspec_)
1456 {
1457  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, tspec_); return; }
1458  orientationidx = ncomp;
1459  SetDiscreteTargetBase(tspec_);
1460  FinalizeSerialDiscreteTargetSpec();
1461 }
1462 
1463 void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec()
1464 {
1465  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1466  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1467 
1468  adapt_eval->SetSerialMetaInfo(*tspec_fes->GetMesh(),
1469  *tspec_fes->FEColl(), ncomp);
1470  adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
1471 
1472  tspec_sav = tspec;
1473 
1474  delete tspec_fesv;
1475  tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
1476  tspec_fes->FEColl(), ncomp);
1477 }
1478 
1479 void DiscreteAdaptTC::SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
1480 {
1481  SetSerialDiscreteTargetSize(tspec_);
1482  FinalizeSerialDiscreteTargetSpec();
1483 }
1484 
1485 
1486 void DiscreteAdaptTC::UpdateTargetSpecification(const Vector &new_x,
1487  bool use_flag)
1488 {
1489  if (use_flag && good_tspec) { return; }
1490 
1491  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1492  adapt_eval->ComputeAtNewPosition(new_x, tspec);
1493  tspec_sav = tspec;
1494 
1495  good_tspec = use_flag;
1496 }
1497 
1498 void DiscreteAdaptTC::UpdateTargetSpecification(Vector &new_x,
1499  Vector &IntData)
1500 {
1501  adapt_eval->ComputeAtNewPosition(new_x, IntData);
1502 }
1503 
1504 void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
1505  ElementTransformation &T,
1506  int dofidx, int dir,
1507  const Vector &IntData)
1508 {
1509  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1510 
1511  Array<int> dofs;
1512  tspec_fes->GetElementDofs(T.ElementNo, dofs);
1513  const int cnt = tspec.Size()/ncomp; // dofs per scalar-field
1514 
1515  for (int i = 0; i < ncomp; i++)
1516  {
1517  tspec(dofs[dofidx]+i*cnt) = IntData(dofs[dofidx] + i*cnt + dir*cnt*ncomp);
1518  }
1519 }
1520 
1521 void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
1522  int dofidx)
1523 {
1524  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1525 
1526  Array<int> dofs;
1527  tspec_fes->GetElementDofs(T.ElementNo, dofs);
1528  const int cnt = tspec.Size()/ncomp;
1529  for (int i = 0; i < ncomp; i++)
1530  {
1531  tspec(dofs[dofidx] + i*cnt) = tspec_sav(dofs[dofidx] + i*cnt);
1532  }
1533 }
1534 
1535 void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
1536  const IntegrationRule &ir,
1537  const Vector &elfun,
1538  DenseTensor &Jtr) const
1539 {
1540  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
1541  const int dim = fe.GetDim(),
1542  nqp = ir.GetNPoints();
1543  Jtrcomp.SetSize(dim, dim, 4*nqp);
1544 
1545  switch (target_type)
1546  {
1547  case IDEAL_SHAPE_GIVEN_SIZE:
1548  case GIVEN_SHAPE_AND_SIZE:
1549  {
1550  const DenseMatrix &Wideal =
1551  Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
1552  const int dim = Wideal.Height(),
1553  ndofs = tspec_fes->GetFE(e_id)->GetDof(),
1554  ntspec_dofs = ndofs*ncomp;
1555 
1556  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
1557  par_vals_c1, par_vals_c2, par_vals_c3;
1558 
1559  Array<int> dofs;
1560  DenseMatrix D_rho(dim), Q_phi(dim), R_theta(dim);
1561  tspec_fesv->GetElementVDofs(e_id, dofs);
1562  tspec.UseDevice(true);
1563  tspec.GetSubVector(dofs, tspec_vals);
1564 
1565  for (int q = 0; q < nqp; q++)
1566  {
1567  const IntegrationPoint &ip = ir.IntPoint(q);
1568  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1569  Jtr(q) = Wideal; // Initialize to identity
1570  for (int d = 0; d < 4; d++)
1571  {
1572  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(d + 4*q), dim, dim);
1573  Jtrcomp_q = Wideal; // Initialize to identity
1574  }
1575 
1576  if (sizeidx != -1) // Set size
1577  {
1578  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
1579  const double min_size = par_vals.Min();
1580  MFEM_VERIFY(min_size > 0.0,
1581  "Non-positive size propagated in the target definition.");
1582  const double size = std::max(shape * par_vals, min_size);
1583  Jtr(q).Set(std::pow(size, 1.0/dim), Jtr(q));
1584  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(0 + 4*q), dim, dim);
1585  Jtrcomp_q = Jtr(q);
1586  } // Done size
1587 
1588  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
1589 
1590  if (aspectratioidx != -1) // Set aspect ratio
1591  {
1592  if (dim == 2)
1593  {
1594  par_vals.SetDataAndSize(tspec_vals.GetData()+
1595  aspectratioidx*ndofs, ndofs);
1596 
1597  const double aspectratio = shape * par_vals;
1598  D_rho = 0.;
1599  D_rho(0,0) = 1./pow(aspectratio,0.5);
1600  D_rho(1,1) = pow(aspectratio,0.5);
1601  }
1602  else
1603  {
1604  par_vals.SetDataAndSize(tspec_vals.GetData()+
1605  aspectratioidx*ndofs, ndofs*3);
1606  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
1607  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
1608  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
1609 
1610  const double rho1 = shape * par_vals_c1;
1611  const double rho2 = shape * par_vals_c2;
1612  const double rho3 = shape * par_vals_c3;
1613  D_rho = 0.;
1614  D_rho(0,0) = pow(rho1,2./3.);
1615  D_rho(1,1) = pow(rho2,2./3.);
1616  D_rho(2,2) = pow(rho3,2./3.);
1617  }
1618  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(1 + 4*q), dim, dim);
1619  Jtrcomp_q = D_rho;
1620  DenseMatrix Temp = Jtr(q);
1621  Mult(D_rho, Temp, Jtr(q));
1622  } // Done aspect ratio
1623 
1624  if (skewidx != -1) // Set skew
1625  {
1626  if (dim == 2)
1627  {
1628  par_vals.SetDataAndSize(tspec_vals.GetData()+
1629  skewidx*ndofs, ndofs);
1630 
1631  const double skew = shape * par_vals;
1632 
1633  Q_phi = 0.;
1634  Q_phi(0,0) = 1.;
1635  Q_phi(0,1) = cos(skew);
1636  Q_phi(1,1) = sin(skew);
1637  }
1638  else
1639  {
1640  par_vals.SetDataAndSize(tspec_vals.GetData()+
1641  skewidx*ndofs, ndofs*3);
1642  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
1643  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
1644  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
1645 
1646  const double phi12 = shape * par_vals_c1;
1647  const double phi13 = shape * par_vals_c2;
1648  const double chi = shape * par_vals_c3;
1649 
1650  Q_phi = 0.;
1651  Q_phi(0,0) = 1.;
1652  Q_phi(0,1) = cos(phi12);
1653  Q_phi(0,2) = cos(phi13);
1654 
1655  Q_phi(1,1) = sin(phi12);
1656  Q_phi(1,2) = sin(phi13)*cos(chi);
1657 
1658  Q_phi(2,2) = sin(phi13)*sin(chi);
1659  }
1660  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*q), dim, dim);
1661  Jtrcomp_q = Q_phi;
1662  DenseMatrix Temp = Jtr(q);
1663  Mult(Q_phi, Temp, Jtr(q));
1664  } // Done skew
1665 
1666  if (orientationidx != -1) // Set orientation
1667  {
1668  if (dim == 2)
1669  {
1670  par_vals.SetDataAndSize(tspec_vals.GetData()+
1671  orientationidx*ndofs, ndofs);
1672 
1673  const double theta = shape * par_vals;
1674  R_theta(0,0) = cos(theta);
1675  R_theta(0,1) = -sin(theta);
1676  R_theta(1,0) = sin(theta);
1677  R_theta(1,1) = cos(theta);
1678  }
1679  else
1680  {
1681  par_vals.SetDataAndSize(tspec_vals.GetData()+
1682  orientationidx*ndofs, ndofs*3);
1683  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
1684  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
1685  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
1686 
1687  const double theta = shape * par_vals_c1;
1688  const double psi = shape * par_vals_c2;
1689  const double beta = shape * par_vals_c3;
1690 
1691  double ct = cos(theta), st = sin(theta),
1692  cp = cos(psi), sp = sin(psi),
1693  cb = cos(beta), sb = sin(beta);
1694 
1695  R_theta = 0.;
1696  R_theta(0,0) = ct*sp;
1697  R_theta(1,0) = st*sp;
1698  R_theta(2,0) = cp;
1699 
1700  R_theta(0,1) = -st*cb + ct*cp*sb;
1701  R_theta(1,1) = ct*cb + st*cp*sb;
1702  R_theta(2,1) = -sp*sb;
1703 
1704  R_theta(0,0) = -st*sb - ct*cp*cb;
1705  R_theta(1,0) = ct*sb - st*cp*cb;
1706  R_theta(2,0) = sp*cb;
1707  }
1708  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(3 + 4*q), dim, dim);
1709  Jtrcomp_q = R_theta;
1710  DenseMatrix Temp = Jtr(q);
1711  Mult(R_theta, Temp, Jtr(q));
1712  } // Done orientation
1713  }
1714  break;
1715  }
1716  default:
1717  MFEM_ABORT("Incompatible target type for discrete adaptation!");
1718  }
1719 }
1720 
1721 void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
1722  const Vector &elfun,
1723  IsoparametricTransformation &Tpr,
1724  DenseTensor &dJtr) const
1725 {
1726  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1727 
1728  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
1729 
1730  dJtr = 0.;
1731  const int e_id = Tpr.ElementNo;
1732  const FiniteElement *fe = Tpr.GetFE();
1733 
1734  switch (target_type)
1735  {
1736  case IDEAL_SHAPE_GIVEN_SIZE:
1737  case GIVEN_SHAPE_AND_SIZE:
1738  {
1739  const DenseMatrix &Wideal =
1740  Geometries.GetGeomToPerfGeomJac(fe->GetGeomType());
1741  const int dim = Wideal.Height(),
1742  ndofs = fe->GetDof(),
1743  ntspec_dofs = ndofs*ncomp;
1744 
1745  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
1746  par_vals_c1(ndofs), par_vals_c2(ndofs), par_vals_c3(ndofs);
1747 
1748  Array<int> dofs;
1749  DenseMatrix dD_rho(dim), dQ_phi(dim), dR_theta(dim);
1750  DenseMatrix dQ_phi13(dim), dQ_phichi(dim); // dQ_phi is used for dQ/dphi12 in 3D
1751  DenseMatrix dR_psi(dim), dR_beta(dim);
1752  tspec_fesv->GetElementVDofs(e_id, dofs);
1753  tspec.GetSubVector(dofs, tspec_vals);
1754 
1755  DenseMatrix grad_e_c1(ndofs, dim),
1756  grad_e_c2(ndofs, dim),
1757  grad_e_c3(ndofs, dim);
1758  Vector grad_ptr_c1(grad_e_c1.GetData(), ndofs*dim),
1759  grad_ptr_c2(grad_e_c2.GetData(), ndofs*dim),
1760  grad_ptr_c3(grad_e_c3.GetData(), ndofs*dim);
1761 
1762  DenseMatrix grad_phys; // This will be (dof x dim, dof).
1763  fe->ProjectGrad(*fe, Tpr, grad_phys);
1764 
1765  for (int i = 0; i < ir.GetNPoints(); i++)
1766  {
1767  const IntegrationPoint &ip = ir.IntPoint(i);
1768  DenseMatrix Jtrcomp_s(Jtrcomp.GetData(0 + 4*i), dim, dim); // size
1769  DenseMatrix Jtrcomp_d(Jtrcomp.GetData(1 + 4*i), dim, dim); // aspect-ratio
1770  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*i), dim, dim); // skew
1771  DenseMatrix Jtrcomp_r(Jtrcomp.GetData(3 + 4*i), dim, dim); // orientation
1772  DenseMatrix work1(dim), work2(dim), work3(dim);
1773 
1774  if (sizeidx != -1) // Set size
1775  {
1776  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
1777 
1778  grad_phys.Mult(par_vals, grad_ptr_c1);
1779  Vector grad_q(dim);
1780  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1781  grad_e_c1.MultTranspose(shape, grad_q);
1782 
1783  const double min_size = par_vals.Min();
1784  MFEM_VERIFY(min_size > 0.0,
1785  "Non-positive size propagated in the target definition.");
1786  const double size = std::max(shape * par_vals, min_size);
1787  double dz_dsize = (1./dim)*pow(size, 1./dim - 1.);
1788 
1789  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
1790  Mult(Jtrcomp_r, work1, work2); // R*Q*D
1791 
1792  for (int d = 0; d < dim; d++)
1793  {
1794  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1795  work1 = Wideal;
1796  work1.Set(dz_dsize, work1); // dz/dsize
1797  work1 *= grad_q(d); // dz/dsize*dsize/dx
1798  AddMult(work1, work2, dJtr_i); // dz/dx*R*Q*D
1799  }
1800  } // Done size
1801 
1802  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
1803 
1804  if (aspectratioidx != -1) // Set aspect ratio
1805  {
1806  if (dim == 2)
1807  {
1808  par_vals.SetDataAndSize(tspec_vals.GetData()+
1809  aspectratioidx*ndofs, ndofs);
1810 
1811  grad_phys.Mult(par_vals, grad_ptr_c1);
1812  Vector grad_q(dim);
1813  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1814  grad_e_c1.MultTranspose(shape, grad_q);
1815 
1816  const double aspectratio = shape * par_vals;
1817  dD_rho = 0.;
1818  dD_rho(0,0) = -0.5*pow(aspectratio,-1.5);
1819  dD_rho(1,1) = 0.5*pow(aspectratio,-0.5);
1820 
1821  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
1822  Mult(work1, Jtrcomp_q, work2); // z*R*Q
1823 
1824  for (int d = 0; d < dim; d++)
1825  {
1826  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1827  work1 = dD_rho;
1828  work1 *= grad_q(d); // work1 = dD/drho*drho/dx
1829  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
1830  }
1831  }
1832  else // 3D
1833  {
1834  par_vals.SetDataAndSize(tspec_vals.GetData()+
1835  aspectratioidx*ndofs, ndofs*3);
1836  par_vals_c1.SetData(par_vals.GetData());
1837  par_vals_c2.SetData(par_vals.GetData()+ndofs);
1838  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
1839 
1840  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
1841  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
1842  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
1843  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
1844  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1845  grad_e_c1.MultTranspose(shape, grad_q1);
1846  grad_e_c2.MultTranspose(shape, grad_q2);
1847  grad_e_c3.MultTranspose(shape, grad_q3);
1848 
1849  const double rho1 = shape * par_vals_c1;
1850  const double rho2 = shape * par_vals_c2;
1851  const double rho3 = shape * par_vals_c3;
1852  dD_rho = 0.;
1853  dD_rho(0,0) = (2./3.)*pow(rho1,-1./3.);
1854  dD_rho(1,1) = (2./3.)*pow(rho2,-1./3.);
1855  dD_rho(2,2) = (2./3.)*pow(rho3,-1./3.);
1856 
1857  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
1858  Mult(work1, Jtrcomp_q, work2); // z*R*Q
1859 
1860 
1861  for (int d = 0; d < dim; d++)
1862  {
1863  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1864  work1 = dD_rho;
1865  work1(0,0) *= grad_q1(d);
1866  work1(1,2) *= grad_q2(d);
1867  work1(2,2) *= grad_q3(d);
1868  // work1 = dD/dx = dD/drho1*drho1/dx + dD/drho2*drho2/dx
1869  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
1870  }
1871  }
1872  } // Done aspect ratio
1873 
1874  if (skewidx != -1) // Set skew
1875  {
1876  if (dim == 2)
1877  {
1878  par_vals.SetDataAndSize(tspec_vals.GetData()+
1879  skewidx*ndofs, ndofs);
1880 
1881  grad_phys.Mult(par_vals, grad_ptr_c1);
1882  Vector grad_q(dim);
1883  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1884  grad_e_c1.MultTranspose(shape, grad_q);
1885 
1886  const double skew = shape * par_vals;
1887 
1888  dQ_phi = 0.;
1889  dQ_phi(0,0) = 1.;
1890  dQ_phi(0,1) = -sin(skew);
1891  dQ_phi(1,1) = cos(skew);
1892 
1893  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
1894 
1895  for (int d = 0; d < dim; d++)
1896  {
1897  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1898  work1 = dQ_phi;
1899  work1 *= grad_q(d); // work1 = dQ/dphi*dphi/dx
1900  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
1901  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
1902  }
1903  }
1904  else
1905  {
1906  par_vals.SetDataAndSize(tspec_vals.GetData()+
1907  skewidx*ndofs, ndofs*3);
1908  par_vals_c1.SetData(par_vals.GetData());
1909  par_vals_c2.SetData(par_vals.GetData()+ndofs);
1910  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
1911 
1912  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
1913  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
1914  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
1915  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
1916  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1917  grad_e_c1.MultTranspose(shape, grad_q1);
1918  grad_e_c2.MultTranspose(shape, grad_q2);
1919  grad_e_c3.MultTranspose(shape, grad_q3);
1920 
1921  const double phi12 = shape * par_vals_c1;
1922  const double phi13 = shape * par_vals_c2;
1923  const double chi = shape * par_vals_c3;
1924 
1925  dQ_phi = 0.;
1926  dQ_phi(0,0) = 1.;
1927  dQ_phi(0,1) = -sin(phi12);
1928  dQ_phi(1,1) = cos(phi12);
1929 
1930  dQ_phi13 = 0.;
1931  dQ_phi13(0,2) = -sin(phi13);
1932  dQ_phi13(1,2) = cos(phi13)*cos(chi);
1933  dQ_phi13(2,2) = cos(phi13)*sin(chi);
1934 
1935  dQ_phichi = 0.;
1936  dQ_phichi(1,2) = -sin(phi13)*sin(chi);
1937  dQ_phichi(2,2) = sin(phi13)*cos(chi);
1938 
1939  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
1940 
1941  for (int d = 0; d < dim; d++)
1942  {
1943  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1944  work1 = dQ_phi;
1945  work1 *= grad_q1(d); // work1 = dQ/dphi12*dphi12/dx
1946  work1.Add(grad_q2(d), dQ_phi13); // + dQ/dphi13*dphi13/dx
1947  work1.Add(grad_q3(d), dQ_phichi); // + dQ/dchi*dchi/dx
1948  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
1949  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
1950  }
1951  }
1952  } // Done skew
1953 
1954  if (orientationidx != -1) // Set orientation
1955  {
1956  if (dim == 2)
1957  {
1958  par_vals.SetDataAndSize(tspec_vals.GetData()+
1959  orientationidx*ndofs, ndofs);
1960 
1961  grad_phys.Mult(par_vals, grad_ptr_c1);
1962  Vector grad_q(dim);
1963  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1964  grad_e_c1.MultTranspose(shape, grad_q);
1965 
1966  const double theta = shape * par_vals;
1967  dR_theta(0,0) = -sin(theta);
1968  dR_theta(0,1) = -cos(theta);
1969  dR_theta(1,0) = cos(theta);
1970  dR_theta(1,1) = -sin(theta);
1971 
1972  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
1973  Mult(Jtrcomp_s, work1, work2); // z*Q*D
1974  for (int d = 0; d < dim; d++)
1975  {
1976  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1977  work1 = dR_theta;
1978  work1 *= grad_q(d); // work1 = dR/dtheta*dtheta/dx
1979  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
1980  }
1981  }
1982  else
1983  {
1984  par_vals.SetDataAndSize(tspec_vals.GetData()+
1985  orientationidx*ndofs, ndofs*3);
1986  par_vals_c1.SetData(par_vals.GetData());
1987  par_vals_c2.SetData(par_vals.GetData()+ndofs);
1988  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
1989 
1990  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
1991  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
1992  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
1993  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
1994  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
1995  grad_e_c1.MultTranspose(shape, grad_q1);
1996  grad_e_c2.MultTranspose(shape, grad_q2);
1997  grad_e_c3.MultTranspose(shape, grad_q3);
1998 
1999  const double theta = shape * par_vals_c1;
2000  const double psi = shape * par_vals_c2;
2001  const double beta = shape * par_vals_c3;
2002 
2003  const double ct = cos(theta), st = sin(theta),
2004  cp = cos(psi), sp = sin(psi),
2005  cb = cos(beta), sb = sin(beta);
2006 
2007  dR_theta = 0.;
2008  dR_theta(0,0) = -st*sp;
2009  dR_theta(1,0) = ct*sp;
2010  dR_theta(2,0) = 0;
2011 
2012  dR_theta(0,1) = -ct*cb - st*cp*sb;
2013  dR_theta(1,1) = -st*cb + ct*cp*sb;
2014  dR_theta(2,1) = 0.;
2015 
2016  dR_theta(0,0) = -ct*sb + st*cp*cb;
2017  dR_theta(1,0) = -st*sb - ct*cp*cb;
2018  dR_theta(2,0) = 0.;
2019 
2020  dR_beta = 0.;
2021  dR_beta(0,0) = 0.;
2022  dR_beta(1,0) = 0.;
2023  dR_beta(2,0) = 0.;
2024 
2025  dR_beta(0,1) = st*sb + ct*cp*cb;
2026  dR_beta(1,1) = -ct*sb + st*cp*cb;
2027  dR_beta(2,1) = -sp*cb;
2028 
2029  dR_beta(0,0) = -st*cb + ct*cp*sb;
2030  dR_beta(1,0) = ct*cb + st*cp*sb;
2031  dR_beta(2,0) = 0.;
2032 
2033  dR_psi = 0.;
2034  dR_psi(0,0) = ct*cp;
2035  dR_psi(1,0) = st*cp;
2036  dR_psi(2,0) = -sp;
2037 
2038  dR_psi(0,1) = 0. - ct*sp*sb;
2039  dR_psi(1,1) = 0. + st*sp*sb;
2040  dR_psi(2,1) = -cp*sb;
2041 
2042  dR_psi(0,0) = 0. + ct*sp*cb;
2043  dR_psi(1,0) = 0. + st*sp*cb;
2044  dR_psi(2,0) = cp*cb;
2045 
2046  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2047  Mult(Jtrcomp_s, work1, work2); // z*Q*D
2048  for (int d = 0; d < dim; d++)
2049  {
2050  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2051  work1 = dR_theta;
2052  work1 *= grad_q1(d); // work1 = dR/dtheta*dtheta/dx
2053  work1.Add(grad_q2(d), dR_psi); // +dR/dpsi*dpsi/dx
2054  work1.Add(grad_q3(d), dR_beta); // +dR/dbeta*dbeta/dx
2055  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
2056  }
2057  }
2058  } // Done orientation
2059  }
2060  break;
2061  }
2062  default:
2063  MFEM_ABORT("Incompatible target type for discrete adaptation!");
2064  }
2065  Jtrcomp.Clear();
2066 }
2067 
2068 void DiscreteAdaptTC::UpdateGradientTargetSpecification(const Vector &x,
2069  const double dx,
2070  bool use_flag)
2071 {
2072  if (use_flag && good_tspec_grad) { return; }
2073 
2074  const int dim = tspec_fes->GetFE(0)->GetDim(),
2075  cnt = x.Size()/dim;
2076 
2077  tspec_pert1h.SetSize(x.Size()*ncomp);
2078 
2079  Vector TSpecTemp;
2080  Vector xtemp = x;
2081  for (int j = 0; j < dim; j++)
2082  {
2083  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) += dx; }
2084 
2085  TSpecTemp.NewDataAndSize(tspec_pert1h.GetData() + j*cnt*ncomp, cnt*ncomp);
2086  UpdateTargetSpecification(xtemp, TSpecTemp);
2087 
2088  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) -= dx; }
2089  }
2090 
2091  good_tspec_grad = use_flag;
2092 }
2093 
2094 void DiscreteAdaptTC::UpdateHessianTargetSpecification(const Vector &x,
2095  double dx, bool use_flag)
2096 {
2097 
2098  if (use_flag && good_tspec_hess) { return; }
2099 
2100  const int dim = tspec_fes->GetFE(0)->GetDim(),
2101  cnt = x.Size()/dim,
2102  totmix = 1+2*(dim-2);
2103 
2104  tspec_pert2h.SetSize(cnt*dim*ncomp);
2105  tspec_pertmix.SetSize(cnt*totmix*ncomp);
2106 
2107  Vector TSpecTemp;
2108  Vector xtemp = x;
2109 
2110  // T(x+2h)
2111  for (int j = 0; j < dim; j++)
2112  {
2113  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) += 2*dx; }
2114 
2115  TSpecTemp.NewDataAndSize(tspec_pert2h.GetData() + j*cnt*ncomp, cnt*ncomp);
2116  UpdateTargetSpecification(xtemp, TSpecTemp);
2117 
2118  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) -= 2*dx; }
2119  }
2120 
2121  // T(x+h,y+h)
2122  int j = 0;
2123  for (int k1 = 0; k1 < dim; k1++)
2124  {
2125  for (int k2 = 0; (k1 != k2) && (k2 < dim); k2++)
2126  {
2127  for (int i = 0; i < cnt; i++)
2128  {
2129  xtemp(k1*cnt+i) += dx;
2130  xtemp(k2*cnt+i) += dx;
2131  }
2132 
2133  TSpecTemp.NewDataAndSize(tspec_pertmix.GetData() + j*cnt*ncomp, cnt*ncomp);
2134  UpdateTargetSpecification(xtemp, TSpecTemp);
2135 
2136  for (int i = 0; i < cnt; i++)
2137  {
2138  xtemp(k1*cnt+i) -= dx;
2139  xtemp(k2*cnt+i) -= dx;
2140  }
2141  j++;
2142  }
2143  }
2144 
2145  good_tspec_hess = use_flag;
2146 }
2147 
2148 void AdaptivityEvaluator::SetSerialMetaInfo(const Mesh &m,
2149  const FiniteElementCollection &fec,
2150  int num_comp)
2151 {
2152  delete fes;
2153  delete mesh;
2154  mesh = new Mesh(m, true);
2155  fes = new FiniteElementSpace(mesh, &fec, num_comp);
2156  dim = fes->GetFE(0)->GetDim();
2157  ncomp = num_comp;
2158 }
2159 
2160 #ifdef MFEM_USE_MPI
2161 void AdaptivityEvaluator::SetParMetaInfo(const ParMesh &m,
2162  const FiniteElementCollection &fec,
2163  int num_comp)
2164 {
2165  delete pfes;
2166  delete pmesh;
2167  pmesh = new ParMesh(m, true);
2168  pfes = new ParFiniteElementSpace(pmesh, &fec, num_comp);
2169  dim = pfes->GetFE(0)->GetDim();
2170  ncomp = num_comp;
2171 }
2172 #endif
2173 
2174 void AdaptivityEvaluator::ClearGeometricFactors()
2175 {
2176 #ifdef MFEM_USE_MPI
2177  if (pmesh) { pmesh->DeleteGeometricFactors(); }
2178 #else
2179  if (mesh) { mesh->DeleteGeometricFactors(); }
2180 #endif
2181 }
2182 
2183 AdaptivityEvaluator::~AdaptivityEvaluator()
2184 {
2185  delete fes;
2186  delete mesh;
2187 #ifdef MFEM_USE_MPI
2188  delete pfes;
2189  delete pmesh;
2190 #endif
2191 }
2192 
2193 void TMOP_Integrator::ReleasePADeviceMemory()
2194 {
2195  if (PA.enabled)
2196  {
2197  PA.H.GetMemory().DeleteDevice();
2198  PA.H0.GetMemory().DeleteDevice();
2199  PA.Jtr.GetMemory().DeleteDevice();
2200  }
2201 }
2202 
2203 TMOP_Integrator::~TMOP_Integrator()
2204 {
2205  delete lim_func;
2206  delete zeta;
2207  for (int i = 0; i < ElemDer.Size(); i++)
2208  {
2209  delete ElemDer[i];
2210  delete ElemPertEnergy[i];
2211  }
2212 }
2213 
2214 void TMOP_Integrator::EnableLimiting(const GridFunction &n0,
2215  const GridFunction &dist, Coefficient &w0,
2216  TMOP_LimiterFunction *lfunc)
2217 {
2218  nodes0 = &n0;
2219  coeff0 = &w0;
2220  lim_dist = &dist;
2221  MFEM_VERIFY(lim_dist->FESpace()->GetVDim() == 1,
2222  "'dist' must be a scalar GridFunction!");
2223 
2224  delete lim_func;
2225  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2226 }
2227 
2228 void TMOP_Integrator::EnableLimiting(const GridFunction &n0, Coefficient &w0,
2229  TMOP_LimiterFunction *lfunc)
2230 {
2231  nodes0 = &n0;
2232  coeff0 = &w0;
2233  lim_dist = NULL;
2234 
2235  delete lim_func;
2236  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2237 }
2238 
2239 void TMOP_Integrator::EnableAdaptiveLimiting(const GridFunction &z0,
2240  Coefficient &coeff,
2241  AdaptivityEvaluator &ae)
2242 {
2243  zeta_0 = &z0;
2244  delete zeta;
2245  zeta = new GridFunction(z0);
2246  coeff_zeta = &coeff;
2247  adapt_eval = &ae;
2248 
2249  adapt_eval->SetSerialMetaInfo(*zeta->FESpace()->GetMesh(),
2250  *zeta->FESpace()->FEColl(), 1);
2251  adapt_eval->SetInitialField
2252  (*zeta->FESpace()->GetMesh()->GetNodes(), *zeta);
2253 }
2254 
2255 #ifdef MFEM_USE_MPI
2256 void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
2257  Coefficient &coeff,
2258  AdaptivityEvaluator &ae)
2259 {
2260  zeta_0 = &z0;
2261  delete zeta;
2262  zeta = new GridFunction(z0);
2263  coeff_zeta = &coeff;
2264  adapt_eval = &ae;
2265 
2266  adapt_eval->SetParMetaInfo(*z0.ParFESpace()->GetParMesh(),
2267  *z0.ParFESpace()->FEColl(), 1);
2268  adapt_eval->SetInitialField
2269  (*zeta->FESpace()->GetMesh()->GetNodes(), *zeta);
2270 }
2271 #endif
2272 
2273 double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
2274  ElementTransformation &T,
2275  const Vector &elfun)
2276 {
2277  const int dof = el.GetDof(), dim = el.GetDim();
2278  double energy;
2279 
2280  // No adaptive limiting terms if this is a FD computation.
2281  const bool adaptive_limiting = (zeta && fd_call_flag == false);
2282 
2283  DSh.SetSize(dof, dim);
2284  Jrt.SetSize(dim);
2285  Jpr.SetSize(dim);
2286  Jpt.SetSize(dim);
2287  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2288 
2289  const IntegrationRule &ir = EnergyIntegrationRule(el);
2290 
2291  energy = 0.0;
2292  DenseTensor Jtr(dim, dim, ir.GetNPoints());
2293  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2294 
2295  // Limited case.
2296  Vector shape, p, p0, d_vals;
2297  DenseMatrix pos0;
2298  if (coeff0)
2299  {
2300  shape.SetSize(dof);
2301  p.SetSize(dim);
2302  p0.SetSize(dim);
2303  pos0.SetSize(dof, dim);
2304  Vector pos0V(pos0.Data(), dof * dim);
2305  Array<int> pos_dofs;
2306  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
2307  nodes0->GetSubVector(pos_dofs, pos0V);
2308  if (lim_dist)
2309  {
2310  lim_dist->GetValues(T.ElementNo, ir, d_vals);
2311  }
2312  else
2313  {
2314  d_vals.SetSize(ir.GetNPoints()); d_vals = 1.0;
2315  }
2316  }
2317 
2318  // Define ref->physical transformation, when a Coefficient is specified.
2319  IsoparametricTransformation *Tpr = NULL;
2320  if (coeff1 || coeff0 || adaptive_limiting)
2321  {
2322  Tpr = new IsoparametricTransformation;
2323  Tpr->SetFE(&el);
2324  Tpr->ElementNo = T.ElementNo;
2325  Tpr->ElementType = ElementTransformation::ELEMENT;
2326  Tpr->Attribute = T.Attribute;
2327  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2328  }
2329  // TODO: computing the coefficients 'coeff1' and 'coeff0' in physical
2330  // coordinates means that, generally, the gradient and Hessian of the
2331  // TMOP_Integrator will depend on the derivatives of the coefficients.
2332  //
2333  // In some cases the coefficients are independent of any movement of
2334  // the physical coordinates (i.e. changes in 'elfun'), e.g. when the
2335  // coefficient is a ConstantCoefficient or a GridFunctionCoefficient.
2336 
2337  Vector zeta_q, zeta0_q;
2338  if (adaptive_limiting)
2339  {
2340  zeta->GetValues(T.ElementNo, ir, zeta_q);
2341  zeta_0->GetValues(T.ElementNo, ir, zeta0_q);
2342  }
2343 
2344  for (int i = 0; i < ir.GetNPoints(); i++)
2345  {
2346  const IntegrationPoint &ip = ir.IntPoint(i);
2347  const DenseMatrix &Jtr_i = Jtr(i);
2348  metric->SetTargetJacobian(Jtr_i);
2349  CalcInverse(Jtr_i, Jrt);
2350  const double weight = ip.weight * Jtr_i.Det();
2351 
2352  el.CalcDShape(ip, DSh);
2353  MultAtB(PMatI, DSh, Jpr);
2354  Mult(Jpr, Jrt, Jpt);
2355 
2356  double val = metric_normal * metric->EvalW(Jpt);
2357  if (coeff1) { val *= coeff1->Eval(*Tpr, ip); }
2358 
2359  if (coeff0)
2360  {
2361  el.CalcShape(ip, shape);
2362  PMatI.MultTranspose(shape, p);
2363  pos0.MultTranspose(shape, p0);
2364  val += lim_normal *
2365  lim_func->Eval(p, p0, d_vals(i)) *
2366  coeff0->Eval(*Tpr, ip);
2367  }
2368 
2369  if (adaptive_limiting)
2370  {
2371  const double diff = zeta_q(i) - zeta0_q(i);
2372  val += coeff_zeta->Eval(*Tpr, ip) * lim_normal * diff * diff;
2373  }
2374 
2375  energy += weight * val;
2376  }
2377  delete Tpr;
2378 
2379  return energy;
2380 }
2381 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
2382  ElementTransformation &T,
2383  const Vector &elfun, Vector &elvect)
2384 {
2385  if (!fdflag)
2386  {
2387  AssembleElementVectorExact(el, T, elfun, elvect);
2388  }
2389  else
2390  {
2391  AssembleElementVectorFD(el, T, elfun, elvect);
2392  }
2393 }
2394 
2395 void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,
2396  ElementTransformation &T,
2397  const Vector &elfun,
2398  DenseMatrix &elmat)
2399 {
2400  if (!fdflag)
2401  {
2402  AssembleElementGradExact(el, T, elfun, elmat);
2403  }
2404  else
2405  {
2406  AssembleElementGradFD(el, T, elfun, elmat);
2407  }
2408 }
2409 
2410 void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
2411  ElementTransformation &T,
2412  const Vector &elfun,
2413  Vector &elvect)
2414 {
2415  const int dof = el.GetDof(), dim = el.GetDim();
2416 
2417  DenseMatrix Amat(dim), work1(dim), work2(dim);
2418  DSh.SetSize(dof, dim);
2419  DS.SetSize(dof, dim);
2420  Jrt.SetSize(dim);
2421  Jpt.SetSize(dim);
2422  P.SetSize(dim);
2423  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2424  elvect.SetSize(dof*dim);
2425  PMatO.UseExternalData(elvect.GetData(), dof, dim);
2426 
2427  const IntegrationRule &ir = ActionIntegrationRule(el);
2428  const int nqp = ir.GetNPoints();
2429 
2430  elvect = 0.0;
2431  Vector weights(nqp);
2432  DenseTensor Jtr(dim, dim, nqp);
2433  DenseTensor dJtr(dim, dim, dim*nqp);
2434  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2435 
2436  // Limited case.
2437  DenseMatrix pos0;
2438  Vector shape, p, p0, d_vals, grad;
2439  shape.SetSize(dof);
2440  if (coeff0)
2441  {
2442  p.SetSize(dim);
2443  p0.SetSize(dim);
2444  pos0.SetSize(dof, dim);
2445  Vector pos0V(pos0.Data(), dof * dim);
2446  Array<int> pos_dofs;
2447  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
2448  nodes0->GetSubVector(pos_dofs, pos0V);
2449  if (lim_dist)
2450  {
2451  lim_dist->GetValues(T.ElementNo, ir, d_vals);
2452  }
2453  else
2454  {
2455  d_vals.SetSize(nqp); d_vals = 1.0;
2456  }
2457  }
2458 
2459  // Define ref->physical transformation, when a Coefficient is specified.
2460  IsoparametricTransformation *Tpr = NULL;
2461  if (coeff1 || coeff0 || zeta || exact_action)
2462  {
2463  Tpr = new IsoparametricTransformation;
2464  Tpr->SetFE(&el);
2465  Tpr->ElementNo = T.ElementNo;
2466  Tpr->ElementType = ElementTransformation::ELEMENT;
2467  Tpr->Attribute = T.Attribute;
2468  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2469  if (exact_action)
2470  {
2471  targetC->ComputeElementTargetsGradient(ir, elfun, *Tpr, dJtr);
2472  }
2473  }
2474 
2475 
2476  Vector d_detW_dx(dim);
2477  Vector d_Winv_dx(dim);
2478 
2479  for (int q = 0; q < nqp; q++)
2480  {
2481  const IntegrationPoint &ip = ir.IntPoint(q);
2482  const DenseMatrix &Jtr_q = Jtr(q);
2483  metric->SetTargetJacobian(Jtr_q);
2484  CalcInverse(Jtr_q, Jrt);
2485  weights(q) = ip.weight * Jtr_q.Det();
2486  double weight_m = weights(q) * metric_normal;
2487 
2488  el.CalcDShape(ip, DSh);
2489  Mult(DSh, Jrt, DS);
2490  MultAtB(PMatI, DS, Jpt);
2491 
2492  metric->EvalP(Jpt, P);
2493 
2494  if (coeff1) { weight_m *= coeff1->Eval(*Tpr, ip); }
2495 
2496  P *= weight_m;
2497  AddMultABt(DS, P, PMatO); // w_q det(W) dmu/dx : dA/dx Winv
2498 
2499  if (exact_action)
2500  {
2501  el.CalcShape(ip, shape);
2502  // Derivatives of adaptivity-based targets.
2503  // First term: w_q d*(Det W)/dx * mu(T)
2504  // d(Det W)/dx = det(W)*Tr[Winv*dW/dx]
2505  DenseMatrix dwdx(dim);
2506  for (int d = 0; d < dim; d++)
2507  {
2508  const DenseMatrix &dJtr_q = dJtr(q + d * nqp);
2509  Mult(Jrt, dJtr_q, dwdx );
2510  d_detW_dx(d) = dwdx.Trace();
2511  }
2512  d_detW_dx *= weight_m*metric->EvalW(Jpt); // *[w_q*det(W)]*mu(T)
2513 
2514  // Second term: w_q det(W) dmu/dx : AdWinv/dx
2515  // dWinv/dx = -Winv*dW/dx*Winv
2516  MultAtB(PMatI, DSh, Amat);
2517  for (int d = 0; d < dim; d++)
2518  {
2519  const DenseMatrix &dJtr_q = dJtr(q + d*nqp);
2520  Mult(Jrt, dJtr_q, work1); // Winv*dw/dx
2521  Mult(work1, Jrt, work2); // Winv*dw/dx*Winv
2522  Mult(Amat, work2, work1); // A*Winv*dw/dx*Winv
2523  MultAtB(P, work1, work2); // dmu/dT^T*A*Winv*dw/dx*Winv
2524  d_Winv_dx(d) = work2.Trace(); // Tr[dmu/dT : AWinv*dw/dx*Winv]
2525  }
2526  d_Winv_dx *= -weight_m; // Include (-) factor as well
2527 
2528  d_detW_dx += d_Winv_dx;
2529  AddMultVWt(shape, d_detW_dx, PMatO);
2530  }
2531 
2532  if (coeff0)
2533  {
2534  if (!exact_action) { el.CalcShape(ip, shape); }
2535  PMatI.MultTranspose(shape, p);
2536  pos0.MultTranspose(shape, p0);
2537  lim_func->Eval_d1(p, p0, d_vals(q), grad);
2538  grad *= weights(q) * lim_normal * coeff0->Eval(*Tpr, ip);
2539  AddMultVWt(shape, grad, PMatO);
2540  }
2541  }
2542 
2543  if (zeta) { AssembleElemVecAdaptLim(el, weights, *Tpr, ir, PMatO); }
2544 
2545  delete Tpr;
2546 }
2547 
2548 void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
2549  ElementTransformation &T,
2550  const Vector &elfun,
2551  DenseMatrix &elmat)
2552 {
2553  const int dof = el.GetDof(), dim = el.GetDim();
2554 
2555  DSh.SetSize(dof, dim);
2556  DS.SetSize(dof, dim);
2557  Jrt.SetSize(dim);
2558  Jpt.SetSize(dim);
2559  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2560  elmat.SetSize(dof*dim);
2561 
2562  const IntegrationRule &ir = GradientIntegrationRule(el);
2563  const int nqp = ir.GetNPoints();
2564 
2565  elmat = 0.0;
2566  Vector weights(nqp);
2567  DenseTensor Jtr(dim, dim, nqp);
2568  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2569 
2570  // Limited case.
2571  DenseMatrix pos0, grad_grad;
2572  Vector shape, p, p0, d_vals;
2573  if (coeff0)
2574  {
2575  shape.SetSize(dof);
2576  p.SetSize(dim);
2577  p0.SetSize(dim);
2578  pos0.SetSize(dof, dim);
2579  Vector pos0V(pos0.Data(), dof * dim);
2580  Array<int> pos_dofs;
2581  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
2582  nodes0->GetSubVector(pos_dofs, pos0V);
2583  if (lim_dist)
2584  {
2585  lim_dist->GetValues(T.ElementNo, ir, d_vals);
2586  }
2587  else
2588  {
2589  d_vals.SetSize(nqp); d_vals = 1.0;
2590  }
2591  }
2592 
2593  // Define ref->physical transformation, when a Coefficient is specified.
2594  IsoparametricTransformation *Tpr = NULL;
2595  if (coeff1 || coeff0 || zeta)
2596  {
2597  Tpr = new IsoparametricTransformation;
2598  Tpr->SetFE(&el);
2599  Tpr->ElementNo = T.ElementNo;
2600  Tpr->ElementType = ElementTransformation::ELEMENT;
2601  Tpr->Attribute = T.Attribute;
2602  Tpr->GetPointMat().Transpose(PMatI);
2603  }
2604 
2605  for (int q = 0; q < nqp; q++)
2606  {
2607  const IntegrationPoint &ip = ir.IntPoint(q);
2608  const DenseMatrix &Jtr_q = Jtr(q);
2609  metric->SetTargetJacobian(Jtr_q);
2610  CalcInverse(Jtr_q, Jrt);
2611  weights(q) = ip.weight * Jtr_q.Det();
2612  double weight_m = weights(q) * metric_normal;
2613 
2614  el.CalcDShape(ip, DSh);
2615  Mult(DSh, Jrt, DS);
2616  MultAtB(PMatI, DS, Jpt);
2617 
2618  if (coeff1) { weight_m *= coeff1->Eval(*Tpr, ip); }
2619 
2620  metric->AssembleH(Jpt, DS, weight_m, elmat);
2621 
2622  // TODO: derivatives of adaptivity-based targets.
2623 
2624  if (coeff0)
2625  {
2626  el.CalcShape(ip, shape);
2627  PMatI.MultTranspose(shape, p);
2628  pos0.MultTranspose(shape, p0);
2629  weight_m = weights(q) * lim_normal * coeff0->Eval(*Tpr, ip);
2630  lim_func->Eval_d2(p, p0, d_vals(q), grad_grad);
2631  for (int i = 0; i < dof; i++)
2632  {
2633  const double w_shape_i = weight_m * shape(i);
2634  for (int j = 0; j < dof; j++)
2635  {
2636  const double w = w_shape_i * shape(j);
2637  for (int d1 = 0; d1 < dim; d1++)
2638  {
2639  for (int d2 = 0; d2 < dim; d2++)
2640  {
2641  elmat(d1*dof + i, d2*dof + j) += w * grad_grad(d1, d2);
2642  }
2643  }
2644  }
2645  }
2646  }
2647  }
2648 
2649  if (zeta) { AssembleElemGradAdaptLim(el, weights, *Tpr, ir, elmat); }
2650 
2651  delete Tpr;
2652 }
2653 
2654 void TMOP_Integrator::AssembleElemVecAdaptLim(const FiniteElement &el,
2655  const Vector &weights,
2656  IsoparametricTransformation &Tpr,
2657  const IntegrationRule &ir,
2658  DenseMatrix &mat)
2659 {
2660  if (zeta == NULL) { return; }
2661 
2662  const int dof = el.GetDof(), dim = el.GetDim();
2663  Vector shape(dof), zeta_e, zeta_q, zeta0_q;
2664 
2665  Array<int> dofs;
2666  zeta->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
2667  zeta->GetSubVector(dofs, zeta_e);
2668  zeta->GetValues(Tpr.ElementNo, ir, zeta_q);
2669  zeta_0->GetValues(Tpr.ElementNo, ir, zeta0_q);
2670 
2671  // Project the gradient of zeta in the same space.
2672  // The FE coefficients of the gradient go in zeta_grad_e.
2673  DenseMatrix zeta_grad_e(dof, dim);
2674  DenseMatrix grad_phys; // This will be (dof x dim, dof).
2675  el.ProjectGrad(el, Tpr, grad_phys);
2676  Vector grad_ptr(zeta_grad_e.GetData(), dof*dim);
2677  grad_phys.Mult(zeta_e, grad_ptr);
2678 
2679  Vector zeta_grad_q(dim);
2680 
2681  const int nqp = weights.Size();
2682  for (int q = 0; q < nqp; q++)
2683  {
2684  const IntegrationPoint &ip = ir.IntPoint(q);
2685  el.CalcShape(ip, shape);
2686  zeta_grad_e.MultTranspose(shape, zeta_grad_q);
2687  zeta_grad_q *= 2.0 * (zeta_q(q) - zeta0_q(q));
2688  zeta_grad_q *= weights(q) * lim_normal * coeff_zeta->Eval(Tpr, ip);
2689  AddMultVWt(shape, zeta_grad_q, mat);
2690  }
2691 }
2692 
2693 void TMOP_Integrator::AssembleElemGradAdaptLim(const FiniteElement &el,
2694  const Vector &weights,
2695  IsoparametricTransformation &Tpr,
2696  const IntegrationRule &ir,
2697  DenseMatrix &mat)
2698 {
2699  if (zeta == NULL) { return; }
2700 
2701  const int dof = el.GetDof(), dim = el.GetDim();
2702  Vector shape(dof), zeta_e, zeta_q, zeta0_q;
2703 
2704  Array<int> dofs;
2705  zeta->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
2706  zeta->GetSubVector(dofs, zeta_e);
2707  zeta->GetValues(Tpr.ElementNo, ir, zeta_q);
2708  zeta_0->GetValues(Tpr.ElementNo, ir, zeta0_q);
2709 
2710  // Project the gradient of zeta in the same space.
2711  // The FE coefficients of the gradient go in zeta_grad_e.
2712  DenseMatrix zeta_grad_e(dof, dim);
2713  DenseMatrix grad_phys; // This will be (dof x dim, dof).
2714  el.ProjectGrad(el, Tpr, grad_phys);
2715  Vector grad_ptr(zeta_grad_e.GetData(), dof*dim);
2716  grad_phys.Mult(zeta_e, grad_ptr);
2717 
2718  // Project the gradient of each gradient of zeta in the same space.
2719  // The FE coefficients of the second derivatives go in zeta_grad_grad_e.
2720  DenseMatrix zeta_grad_grad_e(dof*dim, dim);
2721  Mult(grad_phys, zeta_grad_e, zeta_grad_grad_e);
2722  // Reshape to be more convenient later (no change in the data).
2723  zeta_grad_grad_e.SetSize(dof, dim*dim);
2724 
2725  Vector zeta_grad_q(dim);
2726  DenseMatrix zeta_grad_grad_q(dim, dim);
2727 
2728  const int nqp = weights.Size();
2729  for (int q = 0; q < nqp; q++)
2730  {
2731  const IntegrationPoint &ip = ir.IntPoint(q);
2732  el.CalcShape(ip, shape);
2733 
2734  zeta_grad_e.MultTranspose(shape, zeta_grad_q);
2735  Vector gg_ptr(zeta_grad_grad_q.GetData(), dim*dim);
2736  zeta_grad_grad_e.MultTranspose(shape, gg_ptr);
2737 
2738  const double w = weights(q) * lim_normal * coeff_zeta->Eval(Tpr, ip);
2739  for (int i = 0; i < dof * dim; i++)
2740  {
2741  const int idof = i % dof, idim = i / dof;
2742  for (int j = 0; j <= i; j++)
2743  {
2744  const int jdof = j % dof, jdim = j / dof;
2745  const double entry =
2746  w * ( 2.0 * zeta_grad_q(idim) * shape(idof) *
2747  /* */ zeta_grad_q(jdim) * shape(jdof) +
2748  2.0 * (zeta_q(q) - zeta0_q(q)) *
2749  zeta_grad_grad_q(idim, jdim) * shape(idof) * shape(jdof));
2750  mat(i, j) += entry;
2751  if (i != j) { mat(j, i) += entry; }
2752  }
2753  }
2754  }
2755 }
2756 
2757 double TMOP_Integrator::GetFDDerivative(const FiniteElement &el,
2758  ElementTransformation &T,
2759  Vector &elfun, const int dofidx,
2760  const int dir, const double e_fx,
2761  bool update_stored)
2762 {
2763  int dof = el.GetDof();
2764  int idx = dir*dof+dofidx;
2765  elfun[idx] += dx;
2766  double e_fxph = GetElementEnergy(el, T, elfun);
2767  elfun[idx] -= dx;
2768  double dfdx = (e_fxph-e_fx)/dx;
2769 
2770  if (update_stored)
2771  {
2772  (*(ElemPertEnergy[T.ElementNo]))(idx) = e_fxph;
2773  (*(ElemDer[T.ElementNo]))(idx) = dfdx;
2774  }
2775 
2776  return dfdx;
2777 }
2778 
2779 void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
2780  ElementTransformation &T,
2781  const Vector &elfun,
2782  Vector &elvect)
2783 {
2784  const int dof = el.GetDof(), dim = el.GetDim(), elnum = T.ElementNo;
2785  if (elnum>=ElemDer.Size())
2786  {
2787  ElemDer.Append(new Vector);
2788  ElemPertEnergy.Append(new Vector);
2789  ElemDer[elnum]->SetSize(dof*dim);
2790  ElemPertEnergy[elnum]->SetSize(dof*dim);
2791  }
2792 
2793  elvect.SetSize(dof*dim);
2794  Vector elfunmod(elfun);
2795 
2796  // In GetElementEnergy(), skip terms that have exact derivative calculations.
2797  fd_call_flag = true;
2798 
2799  // Energy for unperturbed configuration.
2800  const double e_fx = GetElementEnergy(el, T, elfun);
2801 
2802  for (int j = 0; j < dim; j++)
2803  {
2804  for (int i = 0; i < dof; i++)
2805  {
2806  if (discr_tc)
2807  {
2808  discr_tc->UpdateTargetSpecificationAtNode(
2809  el, T, i, j, discr_tc->GetTspecPert1H());
2810  }
2811  elvect(j*dof+i) = GetFDDerivative(el, T, elfunmod, i, j, e_fx, true);
2812  if (discr_tc) { discr_tc->RestoreTargetSpecificationAtNode(T, i); }
2813  }
2814  }
2815  fd_call_flag = false;
2816 
2817  // Contributions from adaptive limiting (exact derivatives).
2818  if (zeta)
2819  {
2820  const IntegrationRule &ir = ActionIntegrationRule(el);
2821  const int nqp = ir.GetNPoints();
2822  DenseTensor Jtr(dim, dim, nqp);
2823  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2824 
2825  IsoparametricTransformation Tpr;
2826  Tpr.SetFE(&el);
2827  Tpr.ElementNo = T.ElementNo;
2828  Tpr.Attribute = T.Attribute;
2829  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2830  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2831 
2832  Vector weights(nqp);
2833  for (int q = 0; q < nqp; q++)
2834  {
2835  weights(q) = ir.IntPoint(q).weight * Jtr(q).Det();
2836  }
2837 
2838  PMatO.UseExternalData(elvect.GetData(), dof, dim);
2839  AssembleElemVecAdaptLim(el, weights, Tpr, ir, PMatO);
2840  }
2841 }
2842 
2843 void TMOP_Integrator::AssembleElementGradFD(const FiniteElement &el,
2844  ElementTransformation &T,
2845  const Vector &elfun,
2846  DenseMatrix &elmat)
2847 {
2848  const int dof = el.GetDof(), dim = el.GetDim();
2849 
2850  elmat.SetSize(dof*dim);
2851  Vector elfunmod(elfun);
2852 
2853  const Vector &ElemDerLoc = *(ElemDer[T.ElementNo]);
2854  const Vector &ElemPertLoc = *(ElemPertEnergy[T.ElementNo]);
2855 
2856  // In GetElementEnergy(), skip terms that have exact derivative calculations.
2857  fd_call_flag = true;
2858  for (int i = 0; i < dof; i++)
2859  {
2860  for (int j = 0; j < i+1; j++)
2861  {
2862  for (int k1 = 0; k1 < dim; k1++)
2863  {
2864  for (int k2 = 0; k2 < dim; k2++)
2865  {
2866  elfunmod(k2*dof+j) += dx;
2867 
2868  if (discr_tc)
2869  {
2870  discr_tc->UpdateTargetSpecificationAtNode(
2871  el, T, j, k2, discr_tc->GetTspecPert1H());
2872  if (j != i)
2873  {
2874  discr_tc->UpdateTargetSpecificationAtNode(
2875  el, T, i, k1, discr_tc->GetTspecPert1H());
2876  }
2877  else // j==i
2878  {
2879  if (k1 != k2)
2880  {
2881  int idx = k1+k2-1;
2882  discr_tc->UpdateTargetSpecificationAtNode(
2883  el, T, i, idx, discr_tc->GetTspecPertMixH());
2884  }
2885  else // j==i && k1==k2
2886  {
2887  discr_tc->UpdateTargetSpecificationAtNode(
2888  el, T, i, k1, discr_tc->GetTspecPert2H());
2889  }
2890  }
2891  }
2892 
2893  double e_fx = ElemPertLoc(k2*dof+j);
2894  double e_fpxph = GetFDDerivative(el, T, elfunmod, i, k1, e_fx,
2895  false);
2896  elfunmod(k2*dof+j) -= dx;
2897  double e_fpx = ElemDerLoc(k1*dof+i);
2898 
2899  elmat(k1*dof+i, k2*dof+j) = (e_fpxph - e_fpx) / dx;
2900  elmat(k2*dof+j, k1*dof+i) = (e_fpxph - e_fpx) / dx;
2901 
2902  if (discr_tc)
2903  {
2904  discr_tc->RestoreTargetSpecificationAtNode(T, i);
2905  discr_tc->RestoreTargetSpecificationAtNode(T, j);
2906  }
2907  }
2908  }
2909  }
2910  }
2911  fd_call_flag = false;
2912 
2913  // Contributions from adaptive limiting.
2914  if (zeta)
2915  {
2916  const IntegrationRule &ir = GradientIntegrationRule(el);
2917  const int nqp = ir.GetNPoints();
2918  DenseTensor Jtr(dim, dim, nqp);
2919  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2920 
2921  IsoparametricTransformation Tpr;
2922  Tpr.SetFE(&el);
2923  Tpr.ElementNo = T.ElementNo;
2924  Tpr.Attribute = T.Attribute;
2925  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2926  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2927 
2928  Vector weights(nqp);
2929  for (int q = 0; q < nqp; q++)
2930  {
2931  weights(q) = ir.IntPoint(q).weight * Jtr(q).Det();
2932  }
2933 
2934  AssembleElemGradAdaptLim(el, weights, Tpr, ir, elmat);
2935  }
2936 }
2937 
2938 void TMOP_Integrator::EnableNormalization(const GridFunction &x)
2939 {
2940  ComputeNormalizationEnergies(x, metric_normal, lim_normal);
2941  metric_normal = 1.0 / metric_normal;
2942  lim_normal = 1.0 / lim_normal;
2943 }
2944 
2945 #ifdef MFEM_USE_MPI
2946 void TMOP_Integrator::ParEnableNormalization(const ParGridFunction &x)
2947 {
2948  double loc[2];
2949  ComputeNormalizationEnergies(x, loc[0], loc[1]);
2950  double rdc[2];
2951  MPI_Allreduce(loc, rdc, 2, MPI_DOUBLE, MPI_SUM, x.ParFESpace()->GetComm());
2952  metric_normal = 1.0 / rdc[0];
2953  lim_normal = 1.0 / rdc[1];
2954 }
2955 #endif
2956 
2957 void TMOP_Integrator::ComputeNormalizationEnergies(const GridFunction &x,
2958  double &metric_energy,
2959  double &lim_energy)
2960 {
2961  Array<int> vdofs;
2962  Vector x_vals;
2963  const FiniteElementSpace* const fes = x.FESpace();
2964 
2965  const int dim = fes->GetMesh()->Dimension();
2966  Jrt.SetSize(dim);
2967  Jpr.SetSize(dim);
2968  Jpt.SetSize(dim);
2969 
2970  metric_energy = 0.0;
2971  lim_energy = 0.0;
2972  for (int i = 0; i < fes->GetNE(); i++)
2973  {
2974  const FiniteElement *fe = fes->GetFE(i);
2975  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
2976  const int nqp = ir.GetNPoints();
2977  DenseTensor Jtr(dim, dim, nqp);
2978  const int dof = fe->GetDof();
2979  DSh.SetSize(dof, dim);
2980 
2981  fes->GetElementVDofs(i, vdofs);
2982  x.GetSubVector(vdofs, x_vals);
2983  PMatI.UseExternalData(x_vals.GetData(), dof, dim);
2984 
2985  targetC->ComputeElementTargets(i, *fe, ir, x_vals, Jtr);
2986 
2987  for (int q = 0; q < nqp; q++)
2988  {
2989  const IntegrationPoint &ip = ir.IntPoint(q);
2990  metric->SetTargetJacobian(Jtr(q));
2991  CalcInverse(Jtr(q), Jrt);
2992  const double weight = ip.weight * Jtr(q).Det();
2993 
2994  fe->CalcDShape(ip, DSh);
2995  MultAtB(PMatI, DSh, Jpr);
2996  Mult(Jpr, Jrt, Jpt);
2997 
2998  metric_energy += weight * metric->EvalW(Jpt);
2999  lim_energy += weight;
3000  }
3001  }
3002  if (targetC->ContainsVolumeInfo() == false)
3003  {
3004  // Special case when the targets don't contain volumetric information.
3005  lim_energy = fes->GetNE();
3006  }
3007 }
3008 
3009 void TMOP_Integrator::ComputeMinJac(const Vector &x,
3010  const FiniteElementSpace &fes)
3011 {
3012  const FiniteElement *fe = fes.GetFE(0);
3013  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
3014  const int NE = fes.GetMesh()->GetNE(), dim = fe->GetDim(),
3015  dof = fe->GetDof(), nsp = ir.GetNPoints();
3016 
3017  Array<int> xdofs(dof * dim);
3018  DenseMatrix Jpr(dim), dshape(dof, dim), pos(dof, dim);
3019  Vector posV(pos.Data(), dof * dim);
3020 
3021  dx = std::numeric_limits<float>::max();
3022 
3023  double detv_sum;
3024  double detv_avg_min = std::numeric_limits<float>::max();
3025  for (int i = 0; i < NE; i++)
3026  {
3027  fes.GetElementVDofs(i, xdofs);
3028  x.GetSubVector(xdofs, posV);
3029  detv_sum = 0.;
3030  for (int j = 0; j < nsp; j++)
3031  {
3032  fes.GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
3033  MultAtB(pos, dshape, Jpr);
3034  detv_sum += std::fabs(Jpr.Det());
3035  }
3036  double detv_avg = pow(detv_sum/nsp, 1./dim);
3037  detv_avg_min = std::min(detv_avg, detv_avg_min);
3038  }
3039  dx = detv_avg_min / dxscale;
3040 }
3041 
3042 void TMOP_Integrator::UpdateAfterMeshChange(const Vector &new_x)
3043 {
3044  if (discr_tc)
3045  {
3046  PA.Jtr_needs_update = true;
3047  }
3048  // Update zeta if adaptive limiting is enabled.
3049  if (zeta) { adapt_eval->ComputeAtNewPosition(new_x, *zeta); }
3050 }
3051 
3052 void TMOP_Integrator::ComputeFDh(const Vector &x, const FiniteElementSpace &fes)
3053 {
3054  if (!fdflag) { return; }
3055  ComputeMinJac(x, fes);
3056 }
3057 
3058 #ifdef MFEM_USE_MPI
3059 void TMOP_Integrator::ComputeFDh(const Vector &x,
3060  const ParFiniteElementSpace &pfes)
3061 {
3062  if (!fdflag) { return; }
3063  ComputeMinJac(x, pfes);
3064  double min_jac_all;
3065  MPI_Allreduce(&dx, &min_jac_all, 1, MPI_DOUBLE, MPI_MIN, pfes.GetComm());
3066  dx = min_jac_all;
3067 }
3068 #endif
3069 
3070 void TMOP_Integrator::EnableFiniteDifferences(const GridFunction &x)
3071 {
3072  fdflag = true;
3073  const FiniteElementSpace *fes = x.FESpace();
3074  ComputeFDh(x,*fes);
3075  if (discr_tc)
3076  {
3077  discr_tc->UpdateTargetSpecification(x);
3078  discr_tc->UpdateGradientTargetSpecification(x, dx);
3079  discr_tc->UpdateHessianTargetSpecification(x, dx);
3080  }
3081 }
3082 
3083 #ifdef MFEM_USE_MPI
3084 void TMOP_Integrator::EnableFiniteDifferences(const ParGridFunction &x)
3085 {
3086  fdflag = true;
3087  const ParFiniteElementSpace *pfes = x.ParFESpace();
3088  ComputeFDh(x,*pfes);
3089  if (discr_tc)
3090  {
3091  discr_tc->UpdateTargetSpecification(x);
3092  discr_tc->UpdateGradientTargetSpecification(x, dx);
3093  discr_tc->UpdateHessianTargetSpecification(x, dx);
3094  }
3095 }
3096 #endif
3097 
3098 void TMOPComboIntegrator::EnableLimiting(const GridFunction &n0,
3099  const GridFunction &dist,
3100  Coefficient &w0,
3101  TMOP_LimiterFunction *lfunc)
3102 {
3103  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3104 
3105  tmopi[0]->EnableLimiting(n0, dist, w0, lfunc);
3106  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
3107 }
3108 
3109 void TMOPComboIntegrator::EnableLimiting(const GridFunction &n0,
3110  Coefficient &w0,
3111  TMOP_LimiterFunction *lfunc)
3112 {
3113  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3114 
3115  tmopi[0]->EnableLimiting(n0, w0, lfunc);
3116  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
3117 }
3118 
3119 void TMOPComboIntegrator::SetLimitingNodes(const GridFunction &n0)
3120 {
3121  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3122 
3123  tmopi[0]->SetLimitingNodes(n0);
3124  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
3125 }
3126 
3127 double TMOPComboIntegrator::GetElementEnergy(const FiniteElement &el,
3128  ElementTransformation &T,
3129  const Vector &elfun)
3130 {
3131  double energy= 0.0;
3132  for (int i = 0; i < tmopi.Size(); i++)
3133  {
3134  energy += tmopi[i]->GetElementEnergy(el, T, elfun);
3135  }
3136  return energy;
3137 }
3138 
3139 void TMOPComboIntegrator::AssembleElementVector(const FiniteElement &el,
3140  ElementTransformation &T,
3141  const Vector &elfun,
3142  Vector &elvect)
3143 {
3144  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3145 
3146  tmopi[0]->AssembleElementVector(el, T, elfun, elvect);
3147  for (int i = 1; i < tmopi.Size(); i++)
3148  {
3149  Vector elvect_i;
3150  tmopi[i]->AssembleElementVector(el, T, elfun, elvect_i);
3151  elvect += elvect_i;
3152  }
3153 }
3154 
3155 void TMOPComboIntegrator::AssembleElementGrad(const FiniteElement &el,
3156  ElementTransformation &T,
3157  const Vector &elfun,
3158  DenseMatrix &elmat)
3159 {
3160  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3161 
3162  tmopi[0]->AssembleElementGrad(el, T, elfun, elmat);
3163  for (int i = 1; i < tmopi.Size(); i++)
3164  {
3165  DenseMatrix elmat_i;
3166  tmopi[i]->AssembleElementGrad(el, T, elfun, elmat_i);
3167  elmat += elmat_i;
3168  }
3169 }
3170 
3171 void TMOPComboIntegrator::EnableNormalization(const GridFunction &x)
3172 {
3173  const int cnt = tmopi.Size();
3174  double total_integral = 0.0;
3175  for (int i = 0; i < cnt; i++)
3176  {
3177  tmopi[i]->EnableNormalization(x);
3178  total_integral += 1.0 / tmopi[i]->metric_normal;
3179  }
3180  for (int i = 0; i < cnt; i++)
3181  {
3182  tmopi[i]->metric_normal = 1.0 / total_integral;
3183  }
3184 }
3185 
3186 #ifdef MFEM_USE_MPI
3187 void TMOPComboIntegrator::ParEnableNormalization(const ParGridFunction &x)
3188 {
3189  const int cnt = tmopi.Size();
3190  double total_integral = 0.0;
3191  for (int i = 0; i < cnt; i++)
3192  {
3193  tmopi[i]->ParEnableNormalization(x);
3194  total_integral += 1.0 / tmopi[i]->metric_normal;
3195  }
3196  for (int i = 0; i < cnt; i++)
3197  {
3198  tmopi[i]->metric_normal = 1.0 / total_integral;
3199  }
3200 }
3201 #endif
3202 
3203 void TMOPComboIntegrator::AssemblePA(const FiniteElementSpace &fes)
3204 {
3205  for (int i = 0; i < tmopi.Size(); i++)
3206  {
3207  tmopi[i]->AssemblePA(fes);
3208  }
3209 }
3210 
3211 void TMOPComboIntegrator::AssembleGradPA(const Vector &xe,
3212  const FiniteElementSpace &fes)
3213 {
3214  for (int i = 0; i < tmopi.Size(); i++)
3215  {
3216  tmopi[i]->AssembleGradPA(xe,fes);
3217  }
3218 }
3219 
3220 void TMOPComboIntegrator::AssembleGradDiagonalPA(Vector &de) const
3221 {
3222  for (int i = 0; i < tmopi.Size(); i++)
3223  {
3224  tmopi[i]->AssembleGradDiagonalPA(de);
3225  }
3226 }
3227 
3228 void TMOPComboIntegrator::AddMultPA(const Vector &xe, Vector &ye) const
3229 {
3230  for (int i = 0; i < tmopi.Size(); i++)
3231  {
3232  tmopi[i]->AddMultPA(xe, ye);
3233  }
3234 }
3235 
3236 void TMOPComboIntegrator::AddMultGradPA(const Vector &re, Vector &ce) const
3237 {
3238  for (int i = 0; i < tmopi.Size(); i++)
3239  {
3240  tmopi[i]->AddMultGradPA(re, ce);
3241  }
3242 }
3243 
3244 double TMOPComboIntegrator::GetLocalStateEnergyPA(const Vector &xe) const
3245 {
3246  double energy = 0.0;
3247  for (int i = 0; i < tmopi.Size(); i++)
3248  {
3249  energy += tmopi[i]->GetLocalStateEnergyPA(xe);
3250  }
3251  return energy;
3252 }
3253 
3254 void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric,
3255  const TargetConstructor &tc,
3256  const Mesh &mesh, GridFunction &metric_gf)
3257 {
3258  const int NE = mesh.GetNE();
3259  const GridFunction &nodes = *mesh.GetNodes();
3260  const int dim = mesh.Dimension();
3261  DenseMatrix Winv(dim), T(dim), A(dim), dshape, pos;
3262  Array<int> pos_dofs, gf_dofs;
3263  DenseTensor W;
3264  Vector posV;
3265 
3266  for (int i = 0; i < NE; i++)
3267  {
3268  const FiniteElement &fe_pos = *nodes.FESpace()->GetFE(i);
3269  const IntegrationRule &ir = metric_gf.FESpace()->GetFE(i)->GetNodes();
3270  const int nsp = ir.GetNPoints(), dof = fe_pos.GetDof();
3271 
3272  dshape.SetSize(dof, dim);
3273  pos.SetSize(dof, dim);
3274  posV.SetDataAndSize(pos.Data(), dof * dim);
3275 
3276  metric_gf.FESpace()->GetElementDofs(i, gf_dofs);
3277  nodes.FESpace()->GetElementVDofs(i, pos_dofs);
3278  nodes.GetSubVector(pos_dofs, posV);
3279 
3280  W.SetSize(dim, dim, nsp);
3281  tc.ComputeElementTargets(i, fe_pos, ir, posV, W);
3282 
3283  for (int j = 0; j < nsp; j++)
3284  {
3285  const DenseMatrix &Wj = W(j);
3286  metric.SetTargetJacobian(Wj);
3287  CalcInverse(Wj, Winv);
3288 
3289  const IntegrationPoint &ip = ir.IntPoint(j);
3290  fe_pos.CalcDShape(ip, dshape);
3291  MultAtB(pos, dshape, A);
3292  Mult(A, Winv, T);
3293 
3294  metric_gf(gf_dofs[j]) = metric.EvalW(T);
3295  }
3296  }
3297 }
3298 
3299 } // 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.hpp:243
mfem::TMOP_Metric_077::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:343
mfem::DenseMatrix::Size
int Size() const
For backward compatibility define Size to be synonym of Width()
Definition: densemat.hpp:102
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:344
mfem::TMOP_QualityMetric::Id
virtual int Id() const
Return the metric ID.
Definition: tmop.hpp:74
mfem::DiscreteAdaptTC::SetParDiscreteTargetSize
virtual void SetParDiscreteTargetSize(const ParGridFunction &tspec_)
Definition: tmop.cpp:1347
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:73
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
mfem::TMOP_AMetric_014a::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:961
mfem::TMOP_Metric_022::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:314
mfem::TMOP_Metric_311::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:506
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.hpp:317
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:2776
mfem::AdaptivityEvaluator::SetSerialMetaInfo
void SetSerialMetaInfo(const Mesh &m, const FiniteElementCollection &fec, int num_comp)
Definition: tmop.cpp:2148
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:1009
mfem::TMOP_Integrator::EnableFiniteDifferences
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:3070
mfem::Vector::NewDataAndSize
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:153
mfem::TargetConstructor::target_type
const TargetType target_type
Definition: tmop.hpp:895
mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
mfem::TMOP_Integrator::UpdateAfterMeshChange
void UpdateAfterMeshChange(const Vector &new_x)
Definition: tmop.cpp:3042
mfem::TMOP_Integrator::PMatO
DenseMatrix PMatO
Definition: tmop.hpp:1260
mfem::TMOPComboIntegrator::AddMultPA
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
Definition: tmop.cpp:3228
mfem::TMOP_Integrator::Jtr
DenseTensor Jtr
Definition: tmop.hpp:1290
mfem::DenseTensor::UseExternalData
void UseExternalData(double *ext_data, int i, int j, int k)
Definition: densemat.hpp:804
mfem::TMOP_Metric_022::min_detT
double & min_detT
Definition: tmop.hpp:256
mfem::TMOP_Metric_313::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:525
mfem::TMOP_QualityMetric::Jtr
const DenseMatrix * Jtr
Definition: tmop.hpp:26
mfem::TMOP_Metric_321::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:581
mfem::FiniteElementSpace::IsVariableOrder
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:432
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:920
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:674
mfem::ParFiniteElementSpace::GetParMesh
ParMesh * GetParMesh() const
Definition: pfespace.hpp:267
mfem::VectorCoefficient
Base class for vector Coefficients that optionally depend on time and space.
Definition: coefficient.hpp:342
mfem::InvariantsEvaluator3D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:1087
mfem::TMOP_Integrator::Jrt
DenseMatrix Jrt
Definition: tmop.hpp:1260
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:1222
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetOrientation
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
Definition: tmop.cpp:1454
mfem::TMOPComboIntegrator::GetElementEnergy
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Compute the local energy.
Definition: tmop.cpp:3127
mfem::TMOP_Integrator::coeff0
Coefficient * coeff0
Definition: tmop.hpp:1220
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:5428
mfem::TMOP_Metric_055::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:291
mfem::TMOP_Metric_311::eps
const double eps
Definition: tmop.hpp:505
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:2938
mfem::InvariantsEvaluator3D::Get_dI3b
const scalar_t * Get_dI3b()
Definition: invariants.hpp:780
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:3254
mfem::FiniteElementSpace::GetElementVDofs
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:262
mfem::TMOP_Integrator::metric
TMOP_QualityMetric * metric
Definition: tmop.hpp:1204
mfem::TMOP_Metric_056::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:308
mfem::DenseMatrix::Det
double Det() const
Definition: densemat.cpp:451
mfem::TMOP_LimiterFunction
Base class for limiting functions to be used in class TMOP_Integrator.
Definition: tmop.hpp:758
mfem::TMOP_QualityMetric::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
mfem::TMOP_AMetric_011::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:935
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:1721
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
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:866
mfem::Operator::Width
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
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:793
mfem::TMOP_Integrator::lim_normal
double lim_normal
Definition: tmop.hpp:1226
mfem::Vector::Norml2
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:783
mfem::TMOP_Metric_002::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:202
mfem::TMOP_Integrator::Jpr
DenseMatrix Jpr
Definition: tmop.hpp:1260
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:584
mfem::ElementTransformation::SetIntPoint
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
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:525
mfem::TMOP_AMetric_107a::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:993
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:111
mfem::AnalyticAdaptTC::ComputeElementTargetsGradient
virtual void ComputeElementTargetsGradient(const IntegrationRule &ir, const Vector &elfun, IsoparametricTransformation &Tpr, DenseTensor &dJtr) const
Definition: tmop.cpp:1270
mfem::AdaptivityEvaluator::pfes
ParFiniteElementSpace * pfes
Definition: tmop.hpp:823
mfem::TMOP_Metric_302::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:664
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
mfem::TMOP_Integrator::AssembleElemGradAdaptLim
void AssembleElemGradAdaptLim(const FiniteElement &el, const Vector &weights, IsoparametricTransformation &Tpr, const IntegrationRule &ir, DenseMatrix &m)
Definition: tmop.cpp:2693
mfem::DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec
void FinalizeSerialDiscreteTargetSpec()
Definition: tmop.cpp:1463
mfem::TMOP_Integrator::metric_normal
double metric_normal
Definition: tmop.hpp:1214
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:2239
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:737
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
mfem::TMOPComboIntegrator::EnableNormalization
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
Definition: tmop.cpp:3171
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:596
mfem::TMOP_Integrator::coeff1
Coefficient * coeff1
Definition: tmop.hpp:1212
mfem::TMOP_Metric_313::min_detT
double & min_detT
Definition: tmop.hpp:524
mfem::TMOP_Metric_085::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:520
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:3203
mfem::IsoparametricTransformation::GetFE
const FiniteElement * GetFE() const
Get the current element used to compute the transformations.
Definition: eltrans.hpp:376
mfem::TMOP_Metric_009::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:225
mfem::TMOP_Integrator::ir
const IntegrationRule * ir
Definition: tmop.hpp:1298
mfem::TMOP_Metric_058::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:326
mfem::FiniteElementSpace::GetElementDofs
virtual void GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2298
mfem::TMOP_Integrator::AssembleElementVectorFD
void AssembleElementVectorFD(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Definition: tmop.cpp:2779
mfem::TargetConstructor::nodes
const GridFunction * nodes
Definition: tmop.hpp:892
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSkew
virtual void SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
Definition: tmop.cpp:1437
mfem::TMOP_Integrator::PA
struct mfem::TMOP_Integrator::@13 PA
mfem::DiscreteAdaptTC::UpdateGradientTargetSpecification
void UpdateGradientTargetSpecification(const Vector &x, double dx, bool use_flag=false)
Definition: tmop.cpp:2068
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:430
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.cpp:161
mfem::Vector::GetData
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:199
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:855
mfem::TMOP_QuadraticLimiter
Default limiter function in TMOP_Integrator.
Definition: tmop.hpp:779
mfem::Mesh::GetElementBaseGeometry
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:971
mfem::Vector::UseDevice
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:108
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:535
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:1535
mfem::TMOP_Integrator::AssembleElementGradExact
void AssembleElementGradExact(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Definition: tmop.cpp:2548
mfem::TMOP_Integrator::ElemDer
Array< Vector * > ElemDer
Definition: tmop.hpp:1247
mfem::DiscreteAdaptTC::GetTspecPertMixH
const Vector & GetTspecPertMixH()
Definition: tmop.hpp:1166
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:3098
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).
Definition: tmop.cpp:432
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:1236
mfem::TMOP_Integrator::zeta
GridFunction * zeta
Definition: tmop.hpp:1230
mfem::TMOP_Integrator::adapt_eval
AdaptivityEvaluator * adapt_eval
Definition: tmop.hpp:1232
mfem::TMOP_Integrator::coeff_zeta
Coefficient * coeff_zeta
Definition: tmop.hpp:1231
mfem::AnalyticAdaptTC::scalar_tspec
Coefficient * scalar_tspec
Definition: tmop.hpp:1008
mfem::TMOP_Metric_252::tau0
double & tau0
Definition: tmop.hpp:433
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:620
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio
virtual void SetSerialDiscreteTargetAspectRatio(const GridFunction &tspec_)
Definition: tmop.cpp:1445
mfem::TMOP_Integrator::ReleasePADeviceMemory
void ReleasePADeviceMemory()
Definition: tmop.cpp:2193
mfem::Geometry::GetGeomToPerfGeomJac
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:97
mfem::DiscreteAdaptTC::UpdateTargetSpecification
void UpdateTargetSpecification(const Vector &new_x, bool use_flag=false)
Definition: tmop.cpp:1486
mfem::DiscreteAdaptTC::tspec_fesv
const FiniteElementSpace * tspec_fesv
Definition: tmop.hpp:1067
mfem::TMOP_Metric_311::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:720
mfem::TMOP_Metric_skew3D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:107
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:710
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:448
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:33
mfem::TMOPComboIntegrator::AddMultGradPA
virtual void AddMultGradPA(const Vector &, Vector &) const
Method for partially assembled gradient action.
Definition: tmop.cpp:3236
mfem::TargetConstructor::ComputeElementTargetsGradient
virtual void ComputeElementTargetsGradient(const IntegrationRule &ir, const Vector &elfun, IsoparametricTransformation &Tpr, DenseTensor &dJtr) const
Definition: tmop.cpp:1213
mfem::FiniteElement::GetGeomType
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:320
mfem::Add
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1931
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:567
mfem::TMOP_QualityMetric
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
mfem::Geometries
Geometry Geometries
Definition: fe.cpp:13507
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
mfem::DenseMatrix::MultTranspose
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:203
mfem::AdaptivityEvaluator::~AdaptivityEvaluator
virtual ~AdaptivityEvaluator()
Definition: tmop.cpp:2183
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:2161
mfem::InvariantsEvaluator3D::Assemble_ddI3b
void Assemble_ddI3b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:1250
mfem::TMOP_Metric_315::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:545
mfem::TMOPComboIntegrator::SetLimitingNodes
void SetLimitingNodes(const GridFunction &n0)
Update the original/reference nodes used for limiting.
Definition: tmop.cpp:3119
mfem::DiscreteAdaptTC::SetDiscreteTargetBase
void SetDiscreteTargetBase(const GridFunction &tspec_)
Definition: tmop.cpp:1388
mfem::ParFiniteElementSpace::GetComm
MPI_Comm GetComm() const
Definition: pfespace.hpp:263
mfem::ParFiniteElementSpace::GetFE
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:511
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:822
mfem::AddMult
void AddMult(const DenseMatrix &b, const DenseMatrix &c, DenseMatrix &a)
Matrix matrix multiplication. A += B * C.
Definition: densemat.cpp:2057
mfem::TMOP_Integrator::discr_tc
DiscreteAdaptTC * discr_tc
Definition: tmop.hpp:1234
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:2214
mfem::TMOP_Metric_211::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:554
mfem::TMOP_Metric_352::tau0
double & tau0
Definition: tmop.hpp:645
mfem::DiscreteAdaptTC::UpdateTargetSpecificationAtNode
void UpdateTargetSpecificationAtNode(const FiniteElement &el, ElementTransformation &T, int nodenum, int idir, const Vector &IntData)
Definition: tmop.cpp:1504
mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
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:390
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:46
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:410
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:419
mfem::TMOP_Metric_009::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:264
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:544
mfem::TMOP_Metric_211::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:415
mfem::TMOP_Metric_303::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:698
mfem::TMOP_Integrator::AssembleElementVectorExact
void AssembleElementVectorExact(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Definition: tmop.cpp:2410
mfem::TMOPComboIntegrator::AssembleGradDiagonalPA
virtual void AssembleGradDiagonalPA(Vector &) const
Method for computing the diagonal of the gradient with partial assembly.
Definition: tmop.cpp:3220
mfem::TMOP_Integrator::targetC
const TargetConstructor * targetC
Definition: tmop.hpp:1205
mfem::TMOP_Integrator::ElemPertEnergy
Array< Vector * > ElemPertEnergy
Definition: tmop.hpp:1248
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:563
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::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:478
mfem::TMOP_Metric_022::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:257
mfem::FiniteElement::GetNodes
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe.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:737
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:3211
mfem::TMOP_Metric_aspratio2D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:147
mfem::DiscreteAdaptTC::SetParDiscreteTargetAspectRatio
virtual void SetParDiscreteTargetAspectRatio(const ParGridFunction &tspec_)
Definition: tmop.cpp:1363
mfem::DenseTensor::GetData
double * GetData(int k)
Definition: densemat.hpp:843
mfem::ElementTransformation::Weight
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
mfem::TMOP_Metric_077::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:491
mfem::DenseTensor::SetSize
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:795
mfem::TMOPComboIntegrator::ParEnableNormalization
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:3187
mfem::TMOP_Metric_252::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:434
mfem::NCMesh::GetNumRootElements
int GetNumRootElements()
Return the number of root elements.
Definition: ncmesh.hpp:346
mfem::GridFunction::GetValues
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:466
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:2381
mfem::TMOPComboIntegrator::GetLocalStateEnergyPA
virtual double GetLocalStateEnergyPA(const Vector &) const
Compute the local (to the MPI rank) energy with partial assembly.
Definition: tmop.cpp:3244
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:777
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:534
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:239
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:109
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:67
mfem::Vector::SetData
void SetData(double *d)
Definition: vector.hpp:140
mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629
mfem::TargetConstructor::ContainsVolumeInfo
virtual bool ContainsVolumeInfo() const
Checks if the target matrices contain non-trivial size specification.
Definition: tmop.cpp:1128
mfem::TMOP_Combo_QualityMetric::tmop_q_arr
Array< TMOP_QualityMetric * > tmop_q_arr
Definition: tmop.hpp:81
mfem::Mesh::Dimension
int Dimension() const
Definition: mesh.hpp:911
mfem::DenseMatrix::GetColumn
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1285
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:430
mfem::DiscreteAdaptTC::RestoreTargetSpecificationAtNode
void RestoreTargetSpecificationAtNode(ElementTransformation &T, int nodenum)
Definition: tmop.cpp:1521
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:231
mfem::DiscreteAdaptTC::SetTspecAtIndex
void SetTspecAtIndex(int idx, const GridFunction &tspec_)
Definition: tmop.cpp:1416
mfem::TargetConstructor::ComputeAvgVolume
void ComputeAvgVolume() const
Definition: tmop.cpp:1014
mfem::DenseTensor::SizeI
int SizeI() const
Definition: densemat.hpp:789
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2197
mfem::TMOP_Metric_321::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:844
mfem::TMOP_Integrator::~TMOP_Integrator
~TMOP_Integrator()
Definition: tmop.cpp:2203
omega
double omega
Definition: ex25.cpp:141
mfem::TMOP_Integrator::dxscale
double dxscale
Definition: tmop.hpp:1239
tmop_tools.hpp
mfem::Vector::Min
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1173
mfem::TMOP_Combo_QualityMetric::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:23
mfem::DenseMatrix::Data
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
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::DiscreteAdaptTC::UpdateHessianTargetSpecification
void UpdateHessianTargetSpecification(const Vector &x, double dx, bool use_flag=false)
Definition: tmop.cpp:2094
mfem::DenseMatrix::Transpose
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1374
mfem::DenseMatrix::Trace
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:427
mfem::AddMultABt
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2496
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:2395
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
mfem::FiniteElement::GetDof
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:323
mfem::InvariantsEvaluator3D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:948
mfem::TMOP_Metric_211::eps
const double eps
Definition: tmop.hpp:414
mfem::InvariantsEvaluator3D::Get_dI1
const scalar_t * Get_dI1()
Definition: invariants.hpp:760
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:498
mfem::Coefficient
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
mfem::TMOP_Metric_352::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:646
mfem::FiniteElementCollection
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
mfem::AnalyticAdaptTC::matrix_tspec
TMOPMatrixCoefficient * matrix_tspec
Definition: tmop.hpp:1010
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:1142
mfem::TMOP_Integrator::ActionIntegrationRule
const IntegrationRule & ActionIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1357
mfem::TMOP_Metric_050::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:275
mfem::MultAAt
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2329
mfem::TMOP_Metric_skew2D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:83
mfem::TMOP_Integrator::AssembleElemVecAdaptLim
void AssembleElemVecAdaptLim(const FiniteElement &el, const Vector &weights, IsoparametricTransformation &Tpr, const IntegrationRule &ir, DenseMatrix &m)
Definition: tmop.cpp:2654
mfem::TargetConstructor::UsesPhysicalCoordinates
bool UsesPhysicalCoordinates() const
Return true if the methods ComputeElementTargets(), ComputeAllElementTargets(), and ComputeElementTar...
Definition: tmop.hpp:957
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:214
mfem::DenseTensor::HostWrite
double * HostWrite()
Shortcut for mfem::Write(GetMemory(), TotalSize(), false).
Definition: densemat.hpp:878
mfem::TMOPComboIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: tmop.cpp:3155
mfem::TMOP_Metric_050::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:371
mfem::TMOP_Metric_055::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:403
mfem::InvariantsEvaluator2D::Get_dI1b
const scalar_t * Get_dI1b()
Definition: invariants.hpp:208
mfem::DiscreteAdaptTC::Jtrcomp
DenseTensor Jtrcomp
Definition: tmop.hpp:1062
mfem::Vector::SetDataAndSize
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:147
mfem::TMOP_Metric_252::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:576
mfem::TMOP_Metric_058::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:462
mfem::InvariantsEvaluator3D::Get_dI1b
const scalar_t * Get_dI1b()
Definition: invariants.hpp:764
mfem::TMOP_Metric_313::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:752
mfem::TMOP_Metric_014::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:298
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.hpp:2154
mfem::TMOP_Metric_007::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:224
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:281
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:801
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:470
mfem::TMOP_Integrator::PMatI
DenseMatrix PMatI
Definition: tmop.hpp:1260
mfem::TMOP_Metric_316::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:563
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:272
mfem::DenseTensor::SizeJ
int SizeJ() const
Definition: densemat.hpp:790
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:3052
mfem::InvariantsEvaluator3D::Get_dI2b
const scalar_t * Get_dI2b()
Definition: invariants.hpp:772
mfem::Vector::ReadWrite
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:442
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::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::DiscreteAdaptTC::adapt_eval
AdaptivityEvaluator * adapt_eval
Definition: tmop.hpp:1075
mfem::TMOP_Integrator::dx
double dx
Definition: tmop.hpp:1238
mfem::TMOP_Metric_301::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:453
mfem::InvariantsEvaluator3D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition: invariants.hpp:1294
mfem::TMOP_Integrator::fes
const FiniteElementSpace * fes
Definition: tmop.hpp:1297
mfem::TMOPComboIntegrator::tmopi
Array< TMOP_Integrator * > tmopi
Definition: tmop.hpp:1530
mfem::IsoparametricTransformation
A standard isoparametric element transformation.
Definition: eltrans.hpp:348
mfem::InvariantsEvaluator3D::Get_dI2
const scalar_t * Get_dI2()
Definition: invariants.hpp:768
mfem::TMOP_Integrator::ComputeMinJac
void ComputeMinJac(const Vector &x, const FiniteElementSpace &fes)
Definition: tmop.cpp:3009
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:682
mfem::Mesh::GetElementTransformation
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:347
mfem::TMOP_Integrator::ComputeNormalizationEnergies
void ComputeNormalizationEnergies(const GridFunction &x, double &metric_energy, double &lim_energy)
Definition: tmop.cpp:2957
mfem::TMOP_Integrator::EnergyIntegrationRule
const IntegrationRule & EnergyIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1348
mfem::TMOP_Integrator::GradientIntegrationRule
const IntegrationRule & GradientIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1362
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:2757
mfem::TMOP_Integrator::fdflag
bool fdflag
Definition: tmop.hpp:1237
mfem::AdaptivityEvaluator::SetParMetaInfo
void SetParMetaInfo(const ParMesh &m, const FiniteElementCollection &fec, int num_comp)
Parallel version of SetSerialMetaInfo.
Definition: tmop.cpp:2161
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:440
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:2273
mfem::TMOP_Metric_001::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:61
mfem::IsoparametricTransformation::SetFE
void SetFE(const FiniteElement *FE)
Set the element that will be used to compute the transformations.
Definition: eltrans.hpp:368
mfem::InvariantsEvaluator3D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:790
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:2388
mfem::TMOP_Integrator::zeta_0
const GridFunction * zeta_0
Definition: tmop.hpp:1229
mfem::TargetConstructor::ComputeAllElementTargets_Fallback
void ComputeAllElementTargets_Fallback(const FiniteElementSpace &fes, const IntegrationRule &ir, const Vector &xe, DenseTensor &Jtr) const
Definition: tmop.cpp:1056
mfem::DiscreteAdaptTC::GetTspecPert2H
const Vector & GetTspecPert2H()
Definition: tmop.hpp:1165
mfem::DiscreteAdaptTC::GetTspecPert1H
const Vector & GetTspecPert1H()
Definition: tmop.hpp:1164
mfem::TMOP_Integrator::AssembleElementGradFD
void AssembleElementGradFD(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Definition: tmop.cpp:2843
mfem::TMOP_Metric_316::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:814
mfem::TMOP_Metric_002::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:189
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:629
mfem::DiscreteAdaptTC::SetParDiscreteTargetSpec
virtual void SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
Definition: tmop.cpp:1381
mfem::TMOP_Metric_007::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:207
mfem::Mesh::ncmesh
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:207
mfem::TMOP_Metric_aspratio3D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:167
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:334
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::InvariantsEvaluator2D::Get_dI2
const scalar_t * Get_dI2()
Definition: invariants.hpp:212
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:830
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).
Definition: tmop.cpp:614
mfem::DiscreteAdaptTC::FinalizeParDiscreteTargetSpec
void FinalizeParDiscreteTargetSpec(const ParGridFunction &tspec_)
Definition: tmop.cpp:1317
mfem::MultAtB
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2648
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::FiniteElementSpace::FEColl
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:554
mfem::DenseMatrix::Mult
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:175
mfem::TMOP_Integrator::DSh
DenseMatrix DSh
Definition: tmop.hpp:1260
mfem::Mesh::GetNodes
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
mfem::TMOP_Integrator::DS
DenseMatrix DS
Definition: tmop.hpp:1260
mfem::TMOP_Metric_303::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:487
mfem::TMOP_Integrator::nodes0
const GridFunction * nodes0
Definition: tmop.hpp:1219
mfem::DiscreteAdaptTC::tspec_fes
const FiniteElementSpace * tspec_fes
Definition: tmop.hpp:1066
mfem::TMOP_Metric_098::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:538
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:411
mfem::Mesh::DeleteGeometricFactors
void DeleteGeometricFactors()
Destroy all GeometricFactors stored by the Mesh.
Definition: mesh.cpp:825
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:568
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:729
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSpec
virtual void SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
Definition: tmop.cpp:1479
mfem::DiscreteAdaptTC::SetParDiscreteTargetSkew
virtual void SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
Definition: tmop.cpp:1355
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:382
mfem::TMOP_Integrator::exact_action
bool exact_action
Definition: tmop.hpp:1245
mfem::DiscreteAdaptTC::SetParDiscreteTargetOrientation
virtual void SetParDiscreteTargetOrientation(const ParGridFunction &tspec_)
Definition: tmop.cpp:1372
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:824
mfem::TargetConstructor
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:868
mfem::TMOP_AMetric_036::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:977
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
mfem::IsoparametricTransformation::GetPointMat
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:391
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:208
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:1227
mfem::DenseTensor
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:745
mfem::ParMesh
Class for parallel meshes.
Definition: pmesh.hpp:32
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSize
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
Definition: tmop.cpp:1428
mfem::DenseMatrix::FNorm2
double FNorm2() const
Compute the square of the Frobenius norm of the matrix.
Definition: densemat.hpp:242
mfem::IntRules
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:377
mfem::TMOP_Integrator::lim_func
TMOP_LimiterFunction * lim_func
Definition: tmop.hpp:1224
mfem::TMOP_Metric_315::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:785
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:704
mfem::TMOP_Metric_302::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:469
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:3139
mfem::TMOP_Integrator::dim
int dim
Definition: tmop.hpp:1289
mfem::TMOP_Integrator::Jpt
DenseMatrix Jpt
Definition: tmop.hpp:1260
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:426
mfem::TMOP_Combo_QualityMetric::wt_arr
Array< double > wt_arr
Definition: tmop.hpp:82
mfem::TMOP_Integrator::ParEnableNormalization
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:2946
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:900
mfem::TMOP_Integrator::P
DenseMatrix P
Definition: tmop.hpp:1260
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:912
mfem::AdaptivityEvaluator::ComputeAtNewPosition
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field)=0
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:772
mfem::TargetConstructor::Parallel
bool Parallel() const
Definition: tmop.hpp:900
mfem::TMOP_Integrator::fd_call_flag
bool fd_call_flag
Definition: tmop.hpp:1242
mfem::TMOP_Metric_352::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:892
mfem::AdaptivityEvaluator::fes
FiniteElementSpace * fes
Definition: tmop.hpp:818
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:507
mfem::ParGridFunction::ParFESpace
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:108
mfem::InvariantsEvaluator3D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:734
Generated on Thu Jul 29 2021 15:58:38 for MFEM by   doxygen 1.8.5