MFEM  v4.4.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-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "tmop.hpp"
13 #include "linearform.hpp"
14 #include "pgridfunc.hpp"
15 #include "tmop_tools.hpp"
16 #include "../general/forall.hpp"
17 
18 namespace mfem
19 {
20 
21 // Target-matrix optimization paradigm (TMOP) mesh quality metrics.
22 
23 double TMOP_Combo_QualityMetric::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 &t)
1318 {
1319  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1320  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1321 
1322  ParFiniteElementSpace *ptspec_fes = t.ParFESpace();
1323 
1324  adapt_eval->SetParMetaInfo(*ptspec_fes->GetParMesh(),
1325  *ptspec_fes->FEColl(), ncomp);
1326  adapt_eval->SetInitialField(*ptspec_fes->GetMesh()->GetNodes(), tspec);
1327 
1328  tspec_sav = tspec;
1329 
1330  delete tspec_fesv;
1331  tspec_fesv = new FiniteElementSpace(ptspec_fes->GetMesh(),
1332  ptspec_fes->FEColl(), ncomp);
1333 
1334  delete ptspec_fesv;
1335  ptspec_fesv = new ParFiniteElementSpace(ptspec_fes->GetParMesh(),
1336  ptspec_fes->FEColl(), ncomp);
1337 
1338  delete tspec_pgf;
1339  tspec_pgf = new ParGridFunction(ptspec_fesv, tspec);
1340  tspec_gf = tspec_pgf;
1341 }
1342 
1343 void DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange()
1344 {
1345  ptspec_fesv->Update();
1346  if (tspec_fesv)
1347  {
1348  delete tspec_fesv;
1349  tspec_fesv = new FiniteElementSpace(ptspec_fesv->GetMesh(),
1350  ptspec_fesv->FEColl(), ncomp);
1351  }
1352  tspec_pgf->Update();
1353  tspec_gf = tspec_pgf;
1354  tspec.SetDataAndSize(tspec_pgf->GetData(), tspec_pgf->Size());
1355  tspec_sav = tspec;
1356 
1357  adapt_eval->SetParMetaInfo(*ptspec_fesv->GetParMesh(),
1358  *ptspec_fesv->FEColl(), ncomp);
1359  adapt_eval->SetInitialField(*ptspec_fesv->GetMesh()->GetNodes(), tspec);
1360 }
1361 
1362 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const ParGridFunction &tspec_)
1363 {
1364  const int vdim = tspec_.FESpace()->GetVDim(),
1365  ndof = tspec_.FESpace()->GetNDofs();
1366  MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
1367 
1368  const auto tspec__d = tspec_.Read();
1369  auto tspec_d = tspec.ReadWrite();
1370  const int offset = idx*ndof;
1371  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1372  FinalizeParDiscreteTargetSpec(tspec_);
1373 }
1374 
1375 void DiscreteAdaptTC::SetParDiscreteTargetSize(const ParGridFunction &tspec_)
1376 {
1377  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1378  sizeidx = ncomp;
1379  SetDiscreteTargetBase(tspec_);
1380  FinalizeParDiscreteTargetSpec(tspec_);
1381 }
1382 
1383 void DiscreteAdaptTC::SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
1384 {
1385  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1386  skewidx = ncomp;
1387  SetDiscreteTargetBase(tspec_);
1388  FinalizeParDiscreteTargetSpec(tspec_);
1389 }
1390 
1391 void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction &ar)
1392 {
1393  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
1394  aspectratioidx = ncomp;
1395  SetDiscreteTargetBase(ar);
1396  FinalizeParDiscreteTargetSpec(ar);
1397 }
1398 
1399 void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction &o)
1400 {
1401  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
1402  orientationidx = ncomp;
1403  SetDiscreteTargetBase(o);
1404  FinalizeParDiscreteTargetSpec(o);
1405 }
1406 
1407 void DiscreteAdaptTC::SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
1408 {
1409  SetParDiscreteTargetSize(tspec_);
1410 }
1411 #endif // MFEM_USE_MPI
1412 
1413 void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
1414 {
1415  const int vdim = tspec_.FESpace()->GetVDim(),
1416  ndof = tspec_.FESpace()->GetNDofs();
1417 
1418  ncomp += vdim;
1419 
1420  // need to append data to tspec
1421  // make a copy of tspec->tspec_temp, increase its size, and
1422  // copy data from tspec_temp -> tspec, then add new entries
1423  Vector tspec_temp = tspec;
1424  tspec.UseDevice(true);
1425  tspec_sav.UseDevice(true);
1426  tspec.SetSize(ncomp*ndof);
1427 
1428  const auto tspec_temp_d = tspec_temp.Read();
1429  auto tspec_d = tspec.ReadWrite();
1430  internal::device_copy(tspec_d, tspec_temp_d, tspec_temp.Size());
1431 
1432  const auto tspec__d = tspec_.Read();
1433  const int offset = (ncomp-vdim)*ndof;
1434  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1435 }
1436 
1437 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
1438 {
1439  const int vdim = tspec_.FESpace()->GetVDim(),
1440  ndof = tspec_.FESpace()->GetNDofs();
1441  MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTargetSpec.");
1442 
1443  const auto tspec__d = tspec_.Read();
1444  auto tspec_d = tspec.ReadWrite();
1445  const int offset = idx*ndof;
1446  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1447  FinalizeSerialDiscreteTargetSpec(tspec_);
1448 }
1449 
1450 void DiscreteAdaptTC::SetSerialDiscreteTargetSize(const GridFunction &tspec_)
1451 {
1452  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1453  sizeidx = ncomp;
1454  SetDiscreteTargetBase(tspec_);
1455  FinalizeSerialDiscreteTargetSpec(tspec_);
1456 }
1457 
1458 void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
1459 {
1460  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1461  skewidx = ncomp;
1462  SetDiscreteTargetBase(tspec_);
1463  FinalizeSerialDiscreteTargetSpec(tspec_);
1464 }
1465 
1466 void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(const GridFunction &ar)
1467 {
1468  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
1469  aspectratioidx = ncomp;
1470  SetDiscreteTargetBase(ar);
1471  FinalizeSerialDiscreteTargetSpec(ar);
1472 }
1473 
1474 void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(const GridFunction &o)
1475 {
1476  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
1477  orientationidx = ncomp;
1478  SetDiscreteTargetBase(o);
1479  FinalizeSerialDiscreteTargetSpec(o);
1480 }
1481 
1482 void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec(const GridFunction &t)
1483 {
1484  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1485  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1486 
1487  const FiniteElementSpace *tspec_fes = t.FESpace();
1488  adapt_eval->SetSerialMetaInfo(*tspec_fes->GetMesh(),
1489  *tspec_fes->FEColl(), ncomp);
1490  adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
1491 
1492  tspec_sav = tspec;
1493 
1494  delete tspec_fesv;
1495  tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
1496  tspec_fes->FEColl(), ncomp);
1497 
1498  delete tspec_gf;
1499  tspec_gf = new GridFunction(tspec_fesv, tspec);
1500 }
1501 
1502 void DiscreteAdaptTC::GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
1503 {
1504  if (idx < 0) { return; }
1505  const int ndof = tspec_.FESpace()->GetNDofs(),
1506  vdim = tspec_.FESpace()->GetVDim();
1507  MFEM_VERIFY(ndof == tspec.Size()/ncomp,
1508  "Inconsistency in GetSerialDiscreteTargetSpec.");
1509 
1510  for (int i = 0; i < ndof*vdim; i++)
1511  {
1512  tspec_(i) = tspec(i + idx*ndof);
1513  }
1514 }
1515 
1516 void DiscreteAdaptTC::UpdateAfterMeshTopologyChange()
1517 {
1518  tspec_fesv->Update();
1519  tspec_gf->Update();
1520  tspec.SetDataAndSize(tspec_gf->GetData(), tspec_gf->Size());
1521  tspec_sav = tspec;
1522 
1523  adapt_eval->SetSerialMetaInfo(*tspec_fesv->GetMesh(),
1524  *tspec_fesv->FEColl(), ncomp);
1525  adapt_eval->SetInitialField(*tspec_fesv->GetMesh()->GetNodes(), tspec);
1526 }
1527 
1528 void DiscreteAdaptTC::SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
1529 {
1530  SetSerialDiscreteTargetSize(tspec_);
1531 }
1532 
1533 
1534 void DiscreteAdaptTC::UpdateTargetSpecification(const Vector &new_x,
1535  bool use_flag)
1536 {
1537  if (use_flag && good_tspec) { return; }
1538 
1539  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1540  adapt_eval->ComputeAtNewPosition(new_x, tspec);
1541  tspec_sav = tspec;
1542 
1543  good_tspec = use_flag;
1544 }
1545 
1546 void DiscreteAdaptTC::UpdateTargetSpecification(Vector &new_x,
1547  Vector &IntData)
1548 {
1549  adapt_eval->ComputeAtNewPosition(new_x, IntData);
1550 }
1551 
1552 void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
1553  ElementTransformation &T,
1554  int dofidx, int dir,
1555  const Vector &IntData)
1556 {
1557  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1558 
1559  Array<int> dofs;
1560  tspec_fesv->GetElementDofs(T.ElementNo, dofs);
1561  const int cnt = tspec.Size()/ncomp; // dofs per scalar-field
1562 
1563  for (int i = 0; i < ncomp; i++)
1564  {
1565  tspec(dofs[dofidx]+i*cnt) = IntData(dofs[dofidx] + i*cnt + dir*cnt*ncomp);
1566  }
1567 }
1568 
1569 void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
1570  int dofidx)
1571 {
1572  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
1573 
1574  Array<int> dofs;
1575  tspec_fesv->GetElementDofs(T.ElementNo, dofs);
1576  const int cnt = tspec.Size()/ncomp;
1577  for (int i = 0; i < ncomp; i++)
1578  {
1579  tspec(dofs[dofidx] + i*cnt) = tspec_sav(dofs[dofidx] + i*cnt);
1580  }
1581 }
1582 
1583 void DiscreteAdaptTC::SetTspecFromIntRule(int e_id,
1584  const IntegrationRule &intrule)
1585 {
1586  switch (target_type)
1587  {
1588  case IDEAL_SHAPE_GIVEN_SIZE:
1589  case GIVEN_SHAPE_AND_SIZE:
1590  {
1591  const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
1592  ntspec_dofs = ndofs*ncomp;
1593 
1594  Vector tspec_vals(ntspec_dofs);
1595 
1596  Array<int> dofs;
1597  tspec_fesv->GetElementVDofs(e_id, dofs);
1598  tspec.GetSubVector(dofs, tspec_vals);
1599  DenseMatrix tr;
1600  tspec_gf->GetVectorValues(e_id, intrule, tspec_refine, tr);
1601  tspec_refine.Transpose();
1602  break;
1603  }
1604  default:
1605  MFEM_ABORT("Incompatible target type for discrete adaptation!");
1606  }
1607 }
1608 
1609 void DiscreteAdaptTC::SetTspecDataForDerefinement(FiniteElementSpace *fes)
1610 {
1611  coarse_tspec_fesv = fes;
1612  const Operator *c_op = fes->GetUpdateOperator();
1613  tspec_derefine.SetSize(c_op->Height());
1614  c_op->Mult(tspec, tspec_derefine);
1615 }
1616 
1617 void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
1618  const IntegrationRule &ir,
1619  const Vector &elfun,
1620  DenseTensor &Jtr) const
1621 {
1622  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
1623  const int dim = fe.GetDim(),
1624  nqp = ir.GetNPoints();
1625  Jtrcomp.SetSize(dim, dim, 4*nqp);
1626 
1627  FiniteElementSpace *src_fes = tspec_fesv;
1628 
1629  switch (target_type)
1630  {
1631  case IDEAL_SHAPE_GIVEN_SIZE:
1632  case GIVEN_SHAPE_AND_SIZE:
1633  {
1634  const DenseMatrix &Wideal =
1635  Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
1636  const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
1637  ntspec_dofs = ndofs*ncomp;
1638 
1639  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
1640  par_vals_c1, par_vals_c2, par_vals_c3;
1641 
1642  Array<int> dofs;
1643  DenseMatrix D_rho(dim), Q_phi(dim), R_theta(dim);
1644  tspec_fesv->GetElementVDofs(e_id, dofs);
1645  tspec.UseDevice(true);
1646  tspec.GetSubVector(dofs, tspec_vals);
1647  if (tspec_refine.NumCols() > 0) // Refinement
1648  {
1649  MFEM_VERIFY(amr_el >= 0, " Target being constructed for an AMR element.");
1650  for (int i = 0; i < ncomp; i++)
1651  {
1652  for (int j = 0; j < ndofs; j++)
1653  {
1654  tspec_vals(j + i*ndofs) = tspec_refine(j + amr_el*ndofs, i);
1655  }
1656  }
1657  }
1658  else if (tspec_derefine.Size() > 0) // Derefinement
1659  {
1660  dofs.SetSize(0);
1661  coarse_tspec_fesv->GetElementVDofs(e_id, dofs);
1662  tspec_derefine.GetSubVector(dofs, tspec_vals);
1663  src_fes = coarse_tspec_fesv;
1664  }
1665 
1666  for (int q = 0; q < nqp; q++)
1667  {
1668  const IntegrationPoint &ip = ir.IntPoint(q);
1669  src_fes->GetFE(e_id)->CalcShape(ip, shape);
1670  Jtr(q) = Wideal; // Initialize to identity
1671  for (int d = 0; d < 4; d++)
1672  {
1673  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(d + 4*q), dim, dim);
1674  Jtrcomp_q = Wideal; // Initialize to identity
1675  }
1676 
1677  if (sizeidx != -1) // Set size
1678  {
1679  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
1680  double min_size = par_vals.Min();//0.001; //
1681  if (lim_min_size > 0.)
1682  {
1683  min_size = lim_min_size;
1684  }
1685  else
1686  {
1687  MFEM_VERIFY(min_size > 0.0,
1688  "Non-positive size propagated in the target definition.");
1689  }
1690  const double size = std::max(shape * par_vals, min_size);
1691  Jtr(q).Set(std::pow(size, 1.0/dim), Jtr(q));
1692  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(0 + 4*q), dim, dim);
1693  Jtrcomp_q = Jtr(q);
1694  } // Done size
1695 
1696  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
1697 
1698  if (aspectratioidx != -1) // Set aspect ratio
1699  {
1700  if (dim == 2)
1701  {
1702  par_vals.SetDataAndSize(tspec_vals.GetData()+
1703  aspectratioidx*ndofs, ndofs);
1704  const double min_size = par_vals.Min();
1705  MFEM_VERIFY(min_size > 0.0,
1706  "Non-positive aspect-ratio propagated in the target definition.");
1707 
1708  const double aspectratio = shape * par_vals;
1709  D_rho = 0.;
1710  D_rho(0,0) = 1./pow(aspectratio,0.5);
1711  D_rho(1,1) = pow(aspectratio,0.5);
1712  }
1713  else
1714  {
1715  par_vals.SetDataAndSize(tspec_vals.GetData()+
1716  aspectratioidx*ndofs, ndofs*3);
1717  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
1718  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
1719  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
1720 
1721  const double rho1 = shape * par_vals_c1;
1722  const double rho2 = shape * par_vals_c2;
1723  const double rho3 = shape * par_vals_c3;
1724  D_rho = 0.;
1725  D_rho(0,0) = pow(rho1,2./3.);
1726  D_rho(1,1) = pow(rho2,2./3.);
1727  D_rho(2,2) = pow(rho3,2./3.);
1728  }
1729  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(1 + 4*q), dim, dim);
1730  Jtrcomp_q = D_rho;
1731  DenseMatrix Temp = Jtr(q);
1732  Mult(D_rho, Temp, Jtr(q));
1733  } // Done aspect ratio
1734 
1735  if (skewidx != -1) // Set skew
1736  {
1737  if (dim == 2)
1738  {
1739  par_vals.SetDataAndSize(tspec_vals.GetData()+
1740  skewidx*ndofs, ndofs);
1741 
1742  const double skew = shape * par_vals;
1743 
1744  Q_phi = 0.;
1745  Q_phi(0,0) = 1.;
1746  Q_phi(0,1) = cos(skew);
1747  Q_phi(1,1) = sin(skew);
1748  }
1749  else
1750  {
1751  par_vals.SetDataAndSize(tspec_vals.GetData()+
1752  skewidx*ndofs, ndofs*3);
1753  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
1754  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
1755  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
1756 
1757  const double phi12 = shape * par_vals_c1;
1758  const double phi13 = shape * par_vals_c2;
1759  const double chi = shape * par_vals_c3;
1760 
1761  Q_phi = 0.;
1762  Q_phi(0,0) = 1.;
1763  Q_phi(0,1) = cos(phi12);
1764  Q_phi(0,2) = cos(phi13);
1765 
1766  Q_phi(1,1) = sin(phi12);
1767  Q_phi(1,2) = sin(phi13)*cos(chi);
1768 
1769  Q_phi(2,2) = sin(phi13)*sin(chi);
1770  }
1771  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*q), dim, dim);
1772  Jtrcomp_q = Q_phi;
1773  DenseMatrix Temp = Jtr(q);
1774  Mult(Q_phi, Temp, Jtr(q));
1775  } // Done skew
1776 
1777  if (orientationidx != -1) // Set orientation
1778  {
1779  if (dim == 2)
1780  {
1781  par_vals.SetDataAndSize(tspec_vals.GetData()+
1782  orientationidx*ndofs, ndofs);
1783 
1784  const double theta = shape * par_vals;
1785  R_theta(0,0) = cos(theta);
1786  R_theta(0,1) = -sin(theta);
1787  R_theta(1,0) = sin(theta);
1788  R_theta(1,1) = cos(theta);
1789  }
1790  else
1791  {
1792  par_vals.SetDataAndSize(tspec_vals.GetData()+
1793  orientationidx*ndofs, ndofs*3);
1794  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
1795  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
1796  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
1797 
1798  const double theta = shape * par_vals_c1;
1799  const double psi = shape * par_vals_c2;
1800  const double beta = shape * par_vals_c3;
1801 
1802  double ct = cos(theta), st = sin(theta),
1803  cp = cos(psi), sp = sin(psi),
1804  cb = cos(beta), sb = sin(beta);
1805 
1806  R_theta = 0.;
1807  R_theta(0,0) = ct*sp;
1808  R_theta(1,0) = st*sp;
1809  R_theta(2,0) = cp;
1810 
1811  R_theta(0,1) = -st*cb + ct*cp*sb;
1812  R_theta(1,1) = ct*cb + st*cp*sb;
1813  R_theta(2,1) = -sp*sb;
1814 
1815  R_theta(0,0) = -st*sb - ct*cp*cb;
1816  R_theta(1,0) = ct*sb - st*cp*cb;
1817  R_theta(2,0) = sp*cb;
1818  }
1819  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(3 + 4*q), dim, dim);
1820  Jtrcomp_q = R_theta;
1821  DenseMatrix Temp = Jtr(q);
1822  Mult(R_theta, Temp, Jtr(q));
1823  } // Done orientation
1824  }
1825  break;
1826  }
1827  default:
1828  MFEM_ABORT("Incompatible target type for discrete adaptation!");
1829  }
1830 }
1831 
1832 void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
1833  const Vector &elfun,
1834  IsoparametricTransformation &Tpr,
1835  DenseTensor &dJtr) const
1836 {
1837  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1838 
1839  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
1840 
1841  dJtr = 0.;
1842  const int e_id = Tpr.ElementNo;
1843  const FiniteElement *fe = Tpr.GetFE();
1844 
1845  switch (target_type)
1846  {
1847  case IDEAL_SHAPE_GIVEN_SIZE:
1848  case GIVEN_SHAPE_AND_SIZE:
1849  {
1850  const DenseMatrix &Wideal =
1851  Geometries.GetGeomToPerfGeomJac(fe->GetGeomType());
1852  const int dim = Wideal.Height(),
1853  ndofs = fe->GetDof(),
1854  ntspec_dofs = ndofs*ncomp;
1855 
1856  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
1857  par_vals_c1(ndofs), par_vals_c2(ndofs), par_vals_c3(ndofs);
1858 
1859  Array<int> dofs;
1860  DenseMatrix dD_rho(dim), dQ_phi(dim), dR_theta(dim);
1861  DenseMatrix dQ_phi13(dim), dQ_phichi(dim); // dQ_phi is used for dQ/dphi12 in 3D
1862  DenseMatrix dR_psi(dim), dR_beta(dim);
1863  tspec_fesv->GetElementVDofs(e_id, dofs);
1864  tspec.GetSubVector(dofs, tspec_vals);
1865 
1866  DenseMatrix grad_e_c1(ndofs, dim),
1867  grad_e_c2(ndofs, dim),
1868  grad_e_c3(ndofs, dim);
1869  Vector grad_ptr_c1(grad_e_c1.GetData(), ndofs*dim),
1870  grad_ptr_c2(grad_e_c2.GetData(), ndofs*dim),
1871  grad_ptr_c3(grad_e_c3.GetData(), ndofs*dim);
1872 
1873  DenseMatrix grad_phys; // This will be (dof x dim, dof).
1874  fe->ProjectGrad(*fe, Tpr, grad_phys);
1875 
1876  for (int i = 0; i < ir.GetNPoints(); i++)
1877  {
1878  const IntegrationPoint &ip = ir.IntPoint(i);
1879  DenseMatrix Jtrcomp_s(Jtrcomp.GetData(0 + 4*i), dim, dim); // size
1880  DenseMatrix Jtrcomp_d(Jtrcomp.GetData(1 + 4*i), dim, dim); // aspect-ratio
1881  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*i), dim, dim); // skew
1882  DenseMatrix Jtrcomp_r(Jtrcomp.GetData(3 + 4*i), dim, dim); // orientation
1883  DenseMatrix work1(dim), work2(dim), work3(dim);
1884 
1885  if (sizeidx != -1) // Set size
1886  {
1887  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
1888 
1889  grad_phys.Mult(par_vals, grad_ptr_c1);
1890  Vector grad_q(dim);
1891  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
1892  grad_e_c1.MultTranspose(shape, grad_q);
1893 
1894  const double min_size = par_vals.Min();
1895  MFEM_VERIFY(min_size > 0.0,
1896  "Non-positive size propagated in the target definition.");
1897  const double size = std::max(shape * par_vals, min_size);
1898  double dz_dsize = (1./dim)*pow(size, 1./dim - 1.);
1899 
1900  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
1901  Mult(Jtrcomp_r, work1, work2); // R*Q*D
1902 
1903  for (int d = 0; d < dim; d++)
1904  {
1905  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1906  work1 = Wideal;
1907  work1.Set(dz_dsize, work1); // dz/dsize
1908  work1 *= grad_q(d); // dz/dsize*dsize/dx
1909  AddMult(work1, work2, dJtr_i); // dz/dx*R*Q*D
1910  }
1911  } // Done size
1912 
1913  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
1914 
1915  if (aspectratioidx != -1) // Set aspect ratio
1916  {
1917  if (dim == 2)
1918  {
1919  par_vals.SetDataAndSize(tspec_vals.GetData()+
1920  aspectratioidx*ndofs, ndofs);
1921 
1922  grad_phys.Mult(par_vals, grad_ptr_c1);
1923  Vector grad_q(dim);
1924  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
1925  grad_e_c1.MultTranspose(shape, grad_q);
1926 
1927  const double aspectratio = shape * par_vals;
1928  dD_rho = 0.;
1929  dD_rho(0,0) = -0.5*pow(aspectratio,-1.5);
1930  dD_rho(1,1) = 0.5*pow(aspectratio,-0.5);
1931 
1932  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
1933  Mult(work1, Jtrcomp_q, work2); // z*R*Q
1934 
1935  for (int d = 0; d < dim; d++)
1936  {
1937  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1938  work1 = dD_rho;
1939  work1 *= grad_q(d); // work1 = dD/drho*drho/dx
1940  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
1941  }
1942  }
1943  else // 3D
1944  {
1945  par_vals.SetDataAndSize(tspec_vals.GetData()+
1946  aspectratioidx*ndofs, ndofs*3);
1947  par_vals_c1.SetData(par_vals.GetData());
1948  par_vals_c2.SetData(par_vals.GetData()+ndofs);
1949  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
1950 
1951  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
1952  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
1953  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
1954  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
1955  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
1956  grad_e_c1.MultTranspose(shape, grad_q1);
1957  grad_e_c2.MultTranspose(shape, grad_q2);
1958  grad_e_c3.MultTranspose(shape, grad_q3);
1959 
1960  const double rho1 = shape * par_vals_c1;
1961  const double rho2 = shape * par_vals_c2;
1962  const double rho3 = shape * par_vals_c3;
1963  dD_rho = 0.;
1964  dD_rho(0,0) = (2./3.)*pow(rho1,-1./3.);
1965  dD_rho(1,1) = (2./3.)*pow(rho2,-1./3.);
1966  dD_rho(2,2) = (2./3.)*pow(rho3,-1./3.);
1967 
1968  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
1969  Mult(work1, Jtrcomp_q, work2); // z*R*Q
1970 
1971 
1972  for (int d = 0; d < dim; d++)
1973  {
1974  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1975  work1 = dD_rho;
1976  work1(0,0) *= grad_q1(d);
1977  work1(1,2) *= grad_q2(d);
1978  work1(2,2) *= grad_q3(d);
1979  // work1 = dD/dx = dD/drho1*drho1/dx + dD/drho2*drho2/dx
1980  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
1981  }
1982  }
1983  } // Done aspect ratio
1984 
1985  if (skewidx != -1) // Set skew
1986  {
1987  if (dim == 2)
1988  {
1989  par_vals.SetDataAndSize(tspec_vals.GetData()+
1990  skewidx*ndofs, ndofs);
1991 
1992  grad_phys.Mult(par_vals, grad_ptr_c1);
1993  Vector grad_q(dim);
1994  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
1995  grad_e_c1.MultTranspose(shape, grad_q);
1996 
1997  const double skew = shape * par_vals;
1998 
1999  dQ_phi = 0.;
2000  dQ_phi(0,0) = 1.;
2001  dQ_phi(0,1) = -sin(skew);
2002  dQ_phi(1,1) = cos(skew);
2003 
2004  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
2005 
2006  for (int d = 0; d < dim; d++)
2007  {
2008  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2009  work1 = dQ_phi;
2010  work1 *= grad_q(d); // work1 = dQ/dphi*dphi/dx
2011  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
2012  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
2013  }
2014  }
2015  else
2016  {
2017  par_vals.SetDataAndSize(tspec_vals.GetData()+
2018  skewidx*ndofs, ndofs*3);
2019  par_vals_c1.SetData(par_vals.GetData());
2020  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2021  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2022 
2023  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2024  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2025  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2026  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2027  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2028  grad_e_c1.MultTranspose(shape, grad_q1);
2029  grad_e_c2.MultTranspose(shape, grad_q2);
2030  grad_e_c3.MultTranspose(shape, grad_q3);
2031 
2032  const double phi12 = shape * par_vals_c1;
2033  const double phi13 = shape * par_vals_c2;
2034  const double chi = shape * par_vals_c3;
2035 
2036  dQ_phi = 0.;
2037  dQ_phi(0,0) = 1.;
2038  dQ_phi(0,1) = -sin(phi12);
2039  dQ_phi(1,1) = cos(phi12);
2040 
2041  dQ_phi13 = 0.;
2042  dQ_phi13(0,2) = -sin(phi13);
2043  dQ_phi13(1,2) = cos(phi13)*cos(chi);
2044  dQ_phi13(2,2) = cos(phi13)*sin(chi);
2045 
2046  dQ_phichi = 0.;
2047  dQ_phichi(1,2) = -sin(phi13)*sin(chi);
2048  dQ_phichi(2,2) = sin(phi13)*cos(chi);
2049 
2050  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
2051 
2052  for (int d = 0; d < dim; d++)
2053  {
2054  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2055  work1 = dQ_phi;
2056  work1 *= grad_q1(d); // work1 = dQ/dphi12*dphi12/dx
2057  work1.Add(grad_q2(d), dQ_phi13); // + dQ/dphi13*dphi13/dx
2058  work1.Add(grad_q3(d), dQ_phichi); // + dQ/dchi*dchi/dx
2059  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
2060  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
2061  }
2062  }
2063  } // Done skew
2064 
2065  if (orientationidx != -1) // Set orientation
2066  {
2067  if (dim == 2)
2068  {
2069  par_vals.SetDataAndSize(tspec_vals.GetData()+
2070  orientationidx*ndofs, ndofs);
2071 
2072  grad_phys.Mult(par_vals, grad_ptr_c1);
2073  Vector grad_q(dim);
2074  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2075  grad_e_c1.MultTranspose(shape, grad_q);
2076 
2077  const double theta = shape * par_vals;
2078  dR_theta(0,0) = -sin(theta);
2079  dR_theta(0,1) = -cos(theta);
2080  dR_theta(1,0) = cos(theta);
2081  dR_theta(1,1) = -sin(theta);
2082 
2083  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2084  Mult(Jtrcomp_s, work1, work2); // z*Q*D
2085  for (int d = 0; d < dim; d++)
2086  {
2087  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2088  work1 = dR_theta;
2089  work1 *= grad_q(d); // work1 = dR/dtheta*dtheta/dx
2090  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
2091  }
2092  }
2093  else
2094  {
2095  par_vals.SetDataAndSize(tspec_vals.GetData()+
2096  orientationidx*ndofs, ndofs*3);
2097  par_vals_c1.SetData(par_vals.GetData());
2098  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2099  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2100 
2101  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2102  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2103  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2104  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2105  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2106  grad_e_c1.MultTranspose(shape, grad_q1);
2107  grad_e_c2.MultTranspose(shape, grad_q2);
2108  grad_e_c3.MultTranspose(shape, grad_q3);
2109 
2110  const double theta = shape * par_vals_c1;
2111  const double psi = shape * par_vals_c2;
2112  const double beta = shape * par_vals_c3;
2113 
2114  const double ct = cos(theta), st = sin(theta),
2115  cp = cos(psi), sp = sin(psi),
2116  cb = cos(beta), sb = sin(beta);
2117 
2118  dR_theta = 0.;
2119  dR_theta(0,0) = -st*sp;
2120  dR_theta(1,0) = ct*sp;
2121  dR_theta(2,0) = 0;
2122 
2123  dR_theta(0,1) = -ct*cb - st*cp*sb;
2124  dR_theta(1,1) = -st*cb + ct*cp*sb;
2125  dR_theta(2,1) = 0.;
2126 
2127  dR_theta(0,0) = -ct*sb + st*cp*cb;
2128  dR_theta(1,0) = -st*sb - ct*cp*cb;
2129  dR_theta(2,0) = 0.;
2130 
2131  dR_beta = 0.;
2132  dR_beta(0,0) = 0.;
2133  dR_beta(1,0) = 0.;
2134  dR_beta(2,0) = 0.;
2135 
2136  dR_beta(0,1) = st*sb + ct*cp*cb;
2137  dR_beta(1,1) = -ct*sb + st*cp*cb;
2138  dR_beta(2,1) = -sp*cb;
2139 
2140  dR_beta(0,0) = -st*cb + ct*cp*sb;
2141  dR_beta(1,0) = ct*cb + st*cp*sb;
2142  dR_beta(2,0) = 0.;
2143 
2144  dR_psi = 0.;
2145  dR_psi(0,0) = ct*cp;
2146  dR_psi(1,0) = st*cp;
2147  dR_psi(2,0) = -sp;
2148 
2149  dR_psi(0,1) = 0. - ct*sp*sb;
2150  dR_psi(1,1) = 0. + st*sp*sb;
2151  dR_psi(2,1) = -cp*sb;
2152 
2153  dR_psi(0,0) = 0. + ct*sp*cb;
2154  dR_psi(1,0) = 0. + st*sp*cb;
2155  dR_psi(2,0) = cp*cb;
2156 
2157  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2158  Mult(Jtrcomp_s, work1, work2); // z*Q*D
2159  for (int d = 0; d < dim; d++)
2160  {
2161  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2162  work1 = dR_theta;
2163  work1 *= grad_q1(d); // work1 = dR/dtheta*dtheta/dx
2164  work1.Add(grad_q2(d), dR_psi); // +dR/dpsi*dpsi/dx
2165  work1.Add(grad_q3(d), dR_beta); // +dR/dbeta*dbeta/dx
2166  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
2167  }
2168  }
2169  } // Done orientation
2170  }
2171  break;
2172  }
2173  default:
2174  MFEM_ABORT("Incompatible target type for discrete adaptation!");
2175  }
2176  Jtrcomp.Clear();
2177 }
2178 
2179 void DiscreteAdaptTC::UpdateGradientTargetSpecification(const Vector &x,
2180  const double dx,
2181  bool use_flag)
2182 {
2183  if (use_flag && good_tspec_grad) { return; }
2184 
2185  const int dim = tspec_fesv->GetFE(0)->GetDim(),
2186  cnt = x.Size()/dim;
2187 
2188  tspec_pert1h.SetSize(x.Size()*ncomp);
2189 
2190  Vector TSpecTemp;
2191  Vector xtemp = x;
2192  for (int j = 0; j < dim; j++)
2193  {
2194  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) += dx; }
2195 
2196  TSpecTemp.NewDataAndSize(tspec_pert1h.GetData() + j*cnt*ncomp, cnt*ncomp);
2197  UpdateTargetSpecification(xtemp, TSpecTemp);
2198 
2199  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) -= dx; }
2200  }
2201 
2202  good_tspec_grad = use_flag;
2203 }
2204 
2205 void DiscreteAdaptTC::UpdateHessianTargetSpecification(const Vector &x,
2206  double dx, bool use_flag)
2207 {
2208 
2209  if (use_flag && good_tspec_hess) { return; }
2210 
2211  const int dim = tspec_fesv->GetFE(0)->GetDim(),
2212  cnt = x.Size()/dim,
2213  totmix = 1+2*(dim-2);
2214 
2215  tspec_pert2h.SetSize(cnt*dim*ncomp);
2216  tspec_pertmix.SetSize(cnt*totmix*ncomp);
2217 
2218  Vector TSpecTemp;
2219  Vector xtemp = x;
2220 
2221  // T(x+2h)
2222  for (int j = 0; j < dim; j++)
2223  {
2224  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) += 2*dx; }
2225 
2226  TSpecTemp.NewDataAndSize(tspec_pert2h.GetData() + j*cnt*ncomp, cnt*ncomp);
2227  UpdateTargetSpecification(xtemp, TSpecTemp);
2228 
2229  for (int i = 0; i < cnt; i++) { xtemp(j*cnt+i) -= 2*dx; }
2230  }
2231 
2232  // T(x+h,y+h)
2233  int j = 0;
2234  for (int k1 = 0; k1 < dim; k1++)
2235  {
2236  for (int k2 = 0; (k1 != k2) && (k2 < dim); k2++)
2237  {
2238  for (int i = 0; i < cnt; i++)
2239  {
2240  xtemp(k1*cnt+i) += dx;
2241  xtemp(k2*cnt+i) += dx;
2242  }
2243 
2244  TSpecTemp.NewDataAndSize(tspec_pertmix.GetData() + j*cnt*ncomp, cnt*ncomp);
2245  UpdateTargetSpecification(xtemp, TSpecTemp);
2246 
2247  for (int i = 0; i < cnt; i++)
2248  {
2249  xtemp(k1*cnt+i) -= dx;
2250  xtemp(k2*cnt+i) -= dx;
2251  }
2252  j++;
2253  }
2254  }
2255 
2256  good_tspec_hess = use_flag;
2257 }
2258 
2259 DiscreteAdaptTC::~DiscreteAdaptTC()
2260 {
2261  delete tspec_gf;
2262  delete adapt_eval;
2263  delete tspec_fesv;
2264 #ifdef MFEM_USE_MPI
2265  delete ptspec_fesv;
2266 #endif
2267 }
2268 
2269 void AdaptivityEvaluator::SetSerialMetaInfo(const Mesh &m,
2270  const FiniteElementCollection &fec,
2271  int num_comp)
2272 {
2273  delete fes;
2274  delete mesh;
2275  mesh = new Mesh(m, true);
2276  fes = new FiniteElementSpace(mesh, &fec, num_comp);
2277  dim = fes->GetFE(0)->GetDim();
2278  ncomp = num_comp;
2279 }
2280 
2281 #ifdef MFEM_USE_MPI
2282 void AdaptivityEvaluator::SetParMetaInfo(const ParMesh &m,
2283  const FiniteElementCollection &fec,
2284  int num_comp)
2285 {
2286  delete pfes;
2287  delete pmesh;
2288  pmesh = new ParMesh(m, true);
2289  pfes = new ParFiniteElementSpace(pmesh, &fec, num_comp);
2290  dim = pfes->GetFE(0)->GetDim();
2291  ncomp = num_comp;
2292 }
2293 #endif
2294 
2295 void AdaptivityEvaluator::ClearGeometricFactors()
2296 {
2297 #ifdef MFEM_USE_MPI
2298  if (pmesh) { pmesh->DeleteGeometricFactors(); }
2299 #else
2300  if (mesh) { mesh->DeleteGeometricFactors(); }
2301 #endif
2302 }
2303 
2304 AdaptivityEvaluator::~AdaptivityEvaluator()
2305 {
2306  delete fes;
2307  delete mesh;
2308 #ifdef MFEM_USE_MPI
2309  delete pfes;
2310  delete pmesh;
2311 #endif
2312 }
2313 
2314 void TMOP_Integrator::ReleasePADeviceMemory(bool copy_to_host)
2315 {
2316  if (PA.enabled)
2317  {
2318  PA.H.GetMemory().DeleteDevice(copy_to_host);
2319  PA.H0.GetMemory().DeleteDevice(copy_to_host);
2320  if (!copy_to_host && !PA.Jtr.GetMemory().HostIsValid())
2321  {
2322  PA.Jtr_needs_update = true;
2323  }
2324  PA.Jtr.GetMemory().DeleteDevice(copy_to_host);
2325  }
2326 }
2327 
2328 TMOP_Integrator::~TMOP_Integrator()
2329 {
2330  delete lim_func;
2331  delete adapt_lim_gf;
2332  delete surf_fit_gf;
2333  for (int i = 0; i < ElemDer.Size(); i++)
2334  {
2335  delete ElemDer[i];
2336  delete ElemPertEnergy[i];
2337  }
2338 }
2339 
2340 void TMOP_Integrator::EnableLimiting(const GridFunction &n0,
2341  const GridFunction &dist, Coefficient &w0,
2342  TMOP_LimiterFunction *lfunc)
2343 {
2344  lim_nodes0 = &n0;
2345  lim_coeff = &w0;
2346  lim_dist = &dist;
2347  MFEM_VERIFY(lim_dist->FESpace()->GetVDim() == 1,
2348  "'dist' must be a scalar GridFunction!");
2349 
2350  delete lim_func;
2351  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2352 }
2353 
2354 void TMOP_Integrator::EnableLimiting(const GridFunction &n0, Coefficient &w0,
2355  TMOP_LimiterFunction *lfunc)
2356 {
2357  lim_nodes0 = &n0;
2358  lim_coeff = &w0;
2359  lim_dist = NULL;
2360 
2361  delete lim_func;
2362  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2363 }
2364 
2365 void TMOP_Integrator::EnableAdaptiveLimiting(const GridFunction &z0,
2366  Coefficient &coeff,
2367  AdaptivityEvaluator &ae)
2368 {
2369  adapt_lim_gf0 = &z0;
2370  delete adapt_lim_gf;
2371  adapt_lim_gf = new GridFunction(z0);
2372  adapt_lim_coeff = &coeff;
2373  adapt_lim_eval = &ae;
2374 
2375  adapt_lim_eval->SetSerialMetaInfo(*adapt_lim_gf->FESpace()->GetMesh(),
2376  *adapt_lim_gf->FESpace()->FEColl(), 1);
2377  adapt_lim_eval->SetInitialField
2378  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2379 }
2380 
2381 #ifdef MFEM_USE_MPI
2382 void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
2383  Coefficient &coeff,
2384  AdaptivityEvaluator &ae)
2385 {
2386  adapt_lim_gf0 = &z0;
2387  adapt_lim_pgf0 = &z0;
2388  delete adapt_lim_gf;
2389  adapt_lim_gf = new GridFunction(z0);
2390  adapt_lim_coeff = &coeff;
2391  adapt_lim_eval = &ae;
2392 
2393  adapt_lim_eval->SetParMetaInfo(*z0.ParFESpace()->GetParMesh(),
2394  *z0.ParFESpace()->FEColl(), 1);
2395  adapt_lim_eval->SetInitialField
2396  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2397 }
2398 #endif
2399 
2400 void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &s0,
2401  const Array<bool> &smarker,
2402  Coefficient &coeff,
2403  AdaptivityEvaluator &ae)
2404 {
2405  delete surf_fit_gf;
2406  surf_fit_gf = new GridFunction(s0);
2407  surf_fit_marker = &smarker;
2408  surf_fit_coeff = &coeff;
2409  surf_fit_eval = &ae;
2410 
2411  surf_fit_eval->SetSerialMetaInfo(*s0.FESpace()->GetMesh(),
2412  *s0.FESpace()->FEColl(), 1);
2413  surf_fit_eval->SetInitialField
2414  (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
2415 }
2416 
2417 #ifdef MFEM_USE_MPI
2418 void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
2419  const Array<bool> &smarker,
2420  Coefficient &coeff,
2421  AdaptivityEvaluator &ae)
2422 {
2423  delete surf_fit_gf;
2424  surf_fit_gf = new GridFunction(s0);
2425  surf_fit_marker = &smarker;
2426  surf_fit_coeff = &coeff;
2427  surf_fit_eval = &ae;
2428 
2429  surf_fit_eval->SetParMetaInfo(*s0.ParFESpace()->GetParMesh(),
2430  *s0.ParFESpace()->FEColl(), 1);
2431  surf_fit_eval->SetInitialField
2432  (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
2433 }
2434 #endif
2435 
2436 void TMOP_Integrator::GetSurfaceFittingErrors(double &err_avg, double &err_max)
2437 {
2438  MFEM_VERIFY(surf_fit_gf, "Surface fitting has not been enabled.");
2439 
2440  int loc_cnt = 0;
2441  double loc_max = 0.0, loc_sum = 0.0;
2442  for (int i = 0; i < surf_fit_marker->Size(); i++)
2443  {
2444  if ((*surf_fit_marker)[i] == true)
2445  {
2446  loc_cnt++;
2447  loc_max = std::max(loc_max, std::abs((*surf_fit_gf)(i)));
2448  loc_sum += std::abs((*surf_fit_gf)(i));
2449  }
2450  }
2451  err_avg = loc_sum / loc_cnt;
2452  err_max = loc_max;
2453 
2454 #ifdef MFEM_USE_MPI
2455  if (targetC->Parallel() == false) { return; }
2456  int glob_cnt;
2457  MPI_Comm comm = targetC->GetComm();
2458  MPI_Allreduce(&loc_max, &err_max, 1, MPI_DOUBLE, MPI_MAX, comm);
2459  MPI_Allreduce(&loc_cnt, &glob_cnt, 1, MPI_INT, MPI_SUM, comm);
2460  MPI_Allreduce(&loc_sum, &err_avg, 1, MPI_DOUBLE, MPI_SUM, comm);
2461  err_avg = err_avg / glob_cnt;
2462 #endif
2463 }
2464 
2465 void TMOP_Integrator::UpdateAfterMeshTopologyChange()
2466 {
2467  if (adapt_lim_gf)
2468  {
2469  adapt_lim_gf->Update();
2470  adapt_lim_eval->SetSerialMetaInfo(*adapt_lim_gf->FESpace()->GetMesh(),
2471  *adapt_lim_gf->FESpace()->FEColl(), 1);
2472  adapt_lim_eval->SetInitialField
2473  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2474  }
2475 }
2476 
2477 #ifdef MFEM_USE_MPI
2478 void TMOP_Integrator::ParUpdateAfterMeshTopologyChange()
2479 {
2480  if (adapt_lim_gf)
2481  {
2482  adapt_lim_gf->Update();
2483  adapt_lim_eval->SetParMetaInfo(*adapt_lim_pgf0->ParFESpace()->GetParMesh(),
2484  *adapt_lim_pgf0->ParFESpace()->FEColl(), 1);
2485  adapt_lim_eval->SetInitialField
2486  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2487  }
2488 }
2489 #endif
2490 
2491 double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
2492  ElementTransformation &T,
2493  const Vector &elfun)
2494 {
2495  const int dof = el.GetDof(), dim = el.GetDim();
2496  const int el_id = T.ElementNo;
2497  double energy;
2498 
2499  // No adaptive limiting / surface fitting terms if the function is called
2500  // as part of a FD derivative computation (because we include the exact
2501  // derivatives of these terms in FD computations).
2502  const bool adaptive_limiting = (adapt_lim_gf && fd_call_flag == false);
2503  const bool surface_fit = (surf_fit_gf && fd_call_flag == false);
2504 
2505  DSh.SetSize(dof, dim);
2506  Jrt.SetSize(dim);
2507  Jpr.SetSize(dim);
2508  Jpt.SetSize(dim);
2509  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2510 
2511  const IntegrationRule &ir = EnergyIntegrationRule(el);
2512 
2513  energy = 0.0;
2514  DenseTensor Jtr(dim, dim, ir.GetNPoints());
2515  targetC->ComputeElementTargets(el_id, el, ir, elfun, Jtr);
2516 
2517  // Limited case.
2518  Vector shape, p, p0, d_vals;
2519  DenseMatrix pos0;
2520  if (lim_coeff)
2521  {
2522  shape.SetSize(dof);
2523  p.SetSize(dim);
2524  p0.SetSize(dim);
2525  pos0.SetSize(dof, dim);
2526  Vector pos0V(pos0.Data(), dof * dim);
2527  Array<int> pos_dofs;
2528  lim_nodes0->FESpace()->GetElementVDofs(el_id, pos_dofs);
2529  lim_nodes0->GetSubVector(pos_dofs, pos0V);
2530  if (lim_dist)
2531  {
2532  lim_dist->GetValues(el_id, ir, d_vals);
2533  }
2534  else
2535  {
2536  d_vals.SetSize(ir.GetNPoints()); d_vals = 1.0;
2537  }
2538  }
2539 
2540  // Define ref->physical transformation, when a Coefficient is specified.
2541  IsoparametricTransformation *Tpr = NULL;
2542  if (metric_coeff || lim_coeff || adaptive_limiting || surface_fit)
2543  {
2544  Tpr = new IsoparametricTransformation;
2545  Tpr->SetFE(&el);
2546  Tpr->ElementNo = el_id;
2547  Tpr->ElementType = ElementTransformation::ELEMENT;
2548  Tpr->Attribute = T.Attribute;
2549  Tpr->mesh = T.mesh;
2550  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2551  }
2552  // TODO: computing the coefficients 'metric_coeff' and 'lim_coeff' in physical
2553  // coordinates means that, generally, the gradient and Hessian of the
2554  // TMOP_Integrator will depend on the derivatives of the coefficients.
2555  //
2556  // In some cases the coefficients are independent of any movement of
2557  // the physical coordinates (i.e. changes in 'elfun'), e.g. when the
2558  // coefficient is a ConstantCoefficient or a GridFunctionCoefficient.
2559 
2560  Vector adapt_lim_gf_q, adapt_lim_gf0_q;
2561  if (adaptive_limiting)
2562  {
2563  adapt_lim_gf->GetValues(el_id, ir, adapt_lim_gf_q);
2564  adapt_lim_gf0->GetValues(el_id, ir, adapt_lim_gf0_q);
2565  }
2566 
2567  for (int i = 0; i < ir.GetNPoints(); i++)
2568  {
2569  const IntegrationPoint &ip = ir.IntPoint(i);
2570 
2571  const DenseMatrix &Jtr_i = Jtr(i);
2572  metric->SetTargetJacobian(Jtr_i);
2573  CalcInverse(Jtr_i, Jrt);
2574  const double weight = ip.weight * Jtr_i.Det();
2575 
2576  el.CalcDShape(ip, DSh);
2577  MultAtB(PMatI, DSh, Jpr);
2578  Mult(Jpr, Jrt, Jpt);
2579 
2580  double val = metric_normal * metric->EvalW(Jpt);
2581  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
2582 
2583  if (lim_coeff)
2584  {
2585  el.CalcShape(ip, shape);
2586  PMatI.MultTranspose(shape, p);
2587  pos0.MultTranspose(shape, p0);
2588  val += lim_normal *
2589  lim_func->Eval(p, p0, d_vals(i)) *
2590  lim_coeff->Eval(*Tpr, ip);
2591  }
2592 
2593  // Contribution from the adaptive limiting term.
2594  if (adaptive_limiting)
2595  {
2596  const double diff = adapt_lim_gf_q(i) - adapt_lim_gf0_q(i);
2597  val += adapt_lim_coeff->Eval(*Tpr, ip) * lim_normal * diff * diff;
2598  }
2599 
2600  energy += weight * val;
2601  }
2602 
2603  // Contribution from the surface fitting term.
2604  if (surface_fit)
2605  {
2606  const IntegrationRule &ir_s =
2607  surf_fit_gf->FESpace()->GetFE(el_id)->GetNodes();
2608  Array<int> dofs;
2609  Vector sigma_e;
2610  surf_fit_gf->FESpace()->GetElementDofs(el_id, dofs);
2611  surf_fit_gf->GetSubVector(dofs, sigma_e);
2612  for (int s = 0; s < dofs.Size(); s++)
2613  {
2614  if ((*surf_fit_marker)[dofs[s]] == true)
2615  {
2616  const IntegrationPoint &ip_s = ir_s.IntPoint(s);
2617  Tpr->SetIntPoint(&ip_s);
2618  energy += surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
2619  sigma_e(s) * sigma_e(s);
2620  }
2621  }
2622  }
2623 
2624  delete Tpr;
2625  return energy;
2626 }
2627 
2628 double TMOP_Integrator::GetRefinementElementEnergy(const FiniteElement &el,
2629  ElementTransformation &T,
2630  const Vector &elfun,
2631  const IntegrationRule &irule)
2632 {
2633  int dof = el.GetDof(), dim = el.GetDim(),
2634  NEsplit = elfun.Size() / (dof*dim), el_id = T.ElementNo;
2635  double energy = 0.;
2636 
2637  TargetConstructor *tc = const_cast<TargetConstructor *>(targetC);
2638  DiscreteAdaptTC *dtc = dynamic_cast<DiscreteAdaptTC *>(tc);
2639  // For DiscreteAdaptTC the GridFunctions used to set the targets must be
2640  // mapped onto the fine elements.
2641  if (dtc) { dtc->SetTspecFromIntRule(el_id, irule); }
2642 
2643  for (int e = 0; e < NEsplit; e++)
2644  {
2645  DSh.SetSize(dof, dim);
2646  Jrt.SetSize(dim);
2647  Jpr.SetSize(dim);
2648  Jpt.SetSize(dim);
2649  Vector elfun_child(dof*dim);
2650  for (int i = 0; i < dof; i++)
2651  {
2652  for (int d = 0; d < dim; d++)
2653  {
2654  // elfun is (xe1,xe2,...xen,ye1,ye2...yen) and has nodal coordinates
2655  // for all the children element of the parent element being considered.
2656  // So we must index and get (xek, yek) i.e. nodal coordinates for
2657  // the fine element being considered.
2658  elfun_child(i + d*dof) = elfun(i + e*dof + d*dof*NEsplit);
2659  }
2660  }
2661  PMatI.UseExternalData(elfun_child.GetData(), dof, dim);
2662 
2663  const IntegrationRule &ir = EnergyIntegrationRule(el);
2664 
2665  double el_energy = 0;
2666  DenseTensor Jtr(dim, dim, ir.GetNPoints());
2667  if (dtc)
2668  {
2669  // This is used to index into the tspec vector inside DiscreteAdaptTC.
2670  dtc->SetRefinementSubElement(e);
2671  }
2672  targetC->ComputeElementTargets(el_id, el, ir, elfun_child, Jtr);
2673 
2674  // Define ref->physical transformation, wn a Coefficient is specified.
2675  IsoparametricTransformation *Tpr = NULL;
2676  if (metric_coeff || lim_coeff)
2677  {
2678  Tpr = new IsoparametricTransformation;
2679  Tpr->SetFE(&el);
2680  Tpr->ElementNo = T.ElementNo;
2681  Tpr->ElementType = ElementTransformation::ELEMENT;
2682  Tpr->Attribute = T.Attribute;
2683  Tpr->mesh = T.mesh;
2684  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2685  }
2686 
2687  for (int i = 0; i < ir.GetNPoints(); i++)
2688  {
2689  const IntegrationPoint &ip = ir.IntPoint(i);
2690  const DenseMatrix &Jtr_i = Jtr(i);
2691  h_metric->SetTargetJacobian(Jtr_i);
2692  CalcInverse(Jtr_i, Jrt);
2693  const double weight = ip.weight * Jtr_i.Det();
2694 
2695  el.CalcDShape(ip, DSh);
2696  MultAtB(PMatI, DSh, Jpr);
2697  Mult(Jpr, Jrt, Jpt);
2698 
2699  double val = metric_normal * h_metric->EvalW(Jpt);
2700  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
2701 
2702  el_energy += weight * val;
2703  delete Tpr;
2704  }
2705  energy += el_energy;
2706  }
2707  energy /= NEsplit;
2708 
2709  if (dtc) { dtc->ResetRefinementTspecData(); }
2710 
2711  return energy;
2712 }
2713 
2714 double TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
2715  ElementTransformation &T,
2716  const Vector &elfun)
2717 {
2718  int dof = el.GetDof(), dim = el.GetDim();
2719  double energy = 0.;
2720 
2721  DSh.SetSize(dof, dim);
2722  Jrt.SetSize(dim);
2723  Jpr.SetSize(dim);
2724  Jpt.SetSize(dim);
2725  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2726 
2727  const IntegrationRule &ir = EnergyIntegrationRule(el);
2728 
2729  energy = 0.0;
2730  DenseTensor Jtr(dim, dim, ir.GetNPoints());
2731  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2732 
2733  // Define ref->physical transformation, wn a Coefficient is specified.
2734  IsoparametricTransformation *Tpr = NULL;
2735  if (metric_coeff)
2736  {
2737  Tpr = new IsoparametricTransformation;
2738  Tpr->SetFE(&el);
2739  Tpr->ElementNo = T.ElementNo;
2740  Tpr->ElementType = ElementTransformation::ELEMENT;
2741  Tpr->Attribute = T.Attribute;
2742  Tpr->mesh = T.mesh;
2743  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2744  }
2745 
2746  for (int i = 0; i < ir.GetNPoints(); i++)
2747  {
2748  const IntegrationPoint &ip = ir.IntPoint(i);
2749  const DenseMatrix &Jtr_i = Jtr(i);
2750  h_metric->SetTargetJacobian(Jtr_i);
2751  CalcInverse(Jtr_i, Jrt);
2752  const double weight = ip.weight * Jtr_i.Det();
2753 
2754  el.CalcDShape(ip, DSh);
2755  MultAtB(PMatI, DSh, Jpr);
2756  Mult(Jpr, Jrt, Jpt);
2757 
2758  double val = metric_normal * h_metric->EvalW(Jpt);
2759  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
2760 
2761  energy += weight * val;
2762  }
2763 
2764  delete Tpr;
2765  return energy;
2766 }
2767 
2768 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
2769  ElementTransformation &T,
2770  const Vector &elfun, Vector &elvect)
2771 {
2772  if (!fdflag)
2773  {
2774  AssembleElementVectorExact(el, T, elfun, elvect);
2775  }
2776  else
2777  {
2778  AssembleElementVectorFD(el, T, elfun, elvect);
2779  }
2780 }
2781 
2782 void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,
2783  ElementTransformation &T,
2784  const Vector &elfun,
2785  DenseMatrix &elmat)
2786 {
2787  if (!fdflag)
2788  {
2789  AssembleElementGradExact(el, T, elfun, elmat);
2790  }
2791  else
2792  {
2793  AssembleElementGradFD(el, T, elfun, elmat);
2794  }
2795 }
2796 
2797 void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
2798  ElementTransformation &T,
2799  const Vector &elfun,
2800  Vector &elvect)
2801 {
2802  const int dof = el.GetDof(), dim = el.GetDim();
2803 
2804  DenseMatrix Amat(dim), work1(dim), work2(dim);
2805  DSh.SetSize(dof, dim);
2806  DS.SetSize(dof, dim);
2807  Jrt.SetSize(dim);
2808  Jpt.SetSize(dim);
2809  P.SetSize(dim);
2810  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2811  elvect.SetSize(dof*dim);
2812  PMatO.UseExternalData(elvect.GetData(), dof, dim);
2813 
2814  const IntegrationRule &ir = ActionIntegrationRule(el);
2815  const int nqp = ir.GetNPoints();
2816 
2817  elvect = 0.0;
2818  Vector weights(nqp);
2819  DenseTensor Jtr(dim, dim, nqp);
2820  DenseTensor dJtr(dim, dim, dim*nqp);
2821  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2822 
2823  // Limited case.
2824  DenseMatrix pos0;
2825  Vector shape, p, p0, d_vals, grad;
2826  shape.SetSize(dof);
2827  if (lim_coeff)
2828  {
2829  p.SetSize(dim);
2830  p0.SetSize(dim);
2831  pos0.SetSize(dof, dim);
2832  Vector pos0V(pos0.Data(), dof * dim);
2833  Array<int> pos_dofs;
2834  lim_nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
2835  lim_nodes0->GetSubVector(pos_dofs, pos0V);
2836  if (lim_dist)
2837  {
2838  lim_dist->GetValues(T.ElementNo, ir, d_vals);
2839  }
2840  else
2841  {
2842  d_vals.SetSize(nqp); d_vals = 1.0;
2843  }
2844  }
2845 
2846  // Define ref->physical transformation, when a Coefficient is specified.
2847  IsoparametricTransformation *Tpr = NULL;
2848  if (metric_coeff || lim_coeff || adapt_lim_gf || surf_fit_gf || exact_action)
2849  {
2850  Tpr = new IsoparametricTransformation;
2851  Tpr->SetFE(&el);
2852  Tpr->ElementNo = T.ElementNo;
2853  Tpr->ElementType = ElementTransformation::ELEMENT;
2854  Tpr->Attribute = T.Attribute;
2855  Tpr->mesh = T.mesh;
2856  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
2857  if (exact_action)
2858  {
2859  targetC->ComputeElementTargetsGradient(ir, elfun, *Tpr, dJtr);
2860  }
2861  }
2862 
2863 
2864  Vector d_detW_dx(dim);
2865  Vector d_Winv_dx(dim);
2866 
2867  for (int q = 0; q < nqp; q++)
2868  {
2869  const IntegrationPoint &ip = ir.IntPoint(q);
2870  const DenseMatrix &Jtr_q = Jtr(q);
2871  metric->SetTargetJacobian(Jtr_q);
2872  CalcInverse(Jtr_q, Jrt);
2873  weights(q) = ip.weight * Jtr_q.Det();
2874  double weight_m = weights(q) * metric_normal;
2875 
2876  el.CalcDShape(ip, DSh);
2877  Mult(DSh, Jrt, DS);
2878  MultAtB(PMatI, DS, Jpt);
2879 
2880  metric->EvalP(Jpt, P);
2881 
2882  if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
2883 
2884  P *= weight_m;
2885  AddMultABt(DS, P, PMatO); // w_q det(W) dmu/dx : dA/dx Winv
2886 
2887  if (exact_action)
2888  {
2889  el.CalcShape(ip, shape);
2890  // Derivatives of adaptivity-based targets.
2891  // First term: w_q d*(Det W)/dx * mu(T)
2892  // d(Det W)/dx = det(W)*Tr[Winv*dW/dx]
2893  DenseMatrix dwdx(dim);
2894  for (int d = 0; d < dim; d++)
2895  {
2896  const DenseMatrix &dJtr_q = dJtr(q + d * nqp);
2897  Mult(Jrt, dJtr_q, dwdx );
2898  d_detW_dx(d) = dwdx.Trace();
2899  }
2900  d_detW_dx *= weight_m*metric->EvalW(Jpt); // *[w_q*det(W)]*mu(T)
2901 
2902  // Second term: w_q det(W) dmu/dx : AdWinv/dx
2903  // dWinv/dx = -Winv*dW/dx*Winv
2904  MultAtB(PMatI, DSh, Amat);
2905  for (int d = 0; d < dim; d++)
2906  {
2907  const DenseMatrix &dJtr_q = dJtr(q + d*nqp);
2908  Mult(Jrt, dJtr_q, work1); // Winv*dw/dx
2909  Mult(work1, Jrt, work2); // Winv*dw/dx*Winv
2910  Mult(Amat, work2, work1); // A*Winv*dw/dx*Winv
2911  MultAtB(P, work1, work2); // dmu/dT^T*A*Winv*dw/dx*Winv
2912  d_Winv_dx(d) = work2.Trace(); // Tr[dmu/dT : AWinv*dw/dx*Winv]
2913  }
2914  d_Winv_dx *= -weight_m; // Include (-) factor as well
2915 
2916  d_detW_dx += d_Winv_dx;
2917  AddMultVWt(shape, d_detW_dx, PMatO);
2918  }
2919 
2920  if (lim_coeff)
2921  {
2922  if (!exact_action) { el.CalcShape(ip, shape); }
2923  PMatI.MultTranspose(shape, p);
2924  pos0.MultTranspose(shape, p0);
2925  lim_func->Eval_d1(p, p0, d_vals(q), grad);
2926  grad *= weights(q) * lim_normal * lim_coeff->Eval(*Tpr, ip);
2927  AddMultVWt(shape, grad, PMatO);
2928  }
2929  }
2930 
2931  if (adapt_lim_gf) { AssembleElemVecAdaptLim(el, *Tpr, ir, weights, PMatO); }
2932  if (surf_fit_gf) { AssembleElemVecSurfFit(el, *Tpr, PMatO); }
2933 
2934  delete Tpr;
2935 }
2936 
2937 void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
2938  ElementTransformation &T,
2939  const Vector &elfun,
2940  DenseMatrix &elmat)
2941 {
2942  const int dof = el.GetDof(), dim = el.GetDim();
2943 
2944  DSh.SetSize(dof, dim);
2945  DS.SetSize(dof, dim);
2946  Jrt.SetSize(dim);
2947  Jpt.SetSize(dim);
2948  PMatI.UseExternalData(elfun.GetData(), dof, dim);
2949  elmat.SetSize(dof*dim);
2950 
2951  const IntegrationRule &ir = GradientIntegrationRule(el);
2952  const int nqp = ir.GetNPoints();
2953 
2954  elmat = 0.0;
2955  Vector weights(nqp);
2956  DenseTensor Jtr(dim, dim, nqp);
2957  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
2958 
2959  // Limited case.
2960  DenseMatrix pos0, hess;
2961  Vector shape, p, p0, d_vals;
2962  if (lim_coeff)
2963  {
2964  shape.SetSize(dof);
2965  p.SetSize(dim);
2966  p0.SetSize(dim);
2967  pos0.SetSize(dof, dim);
2968  Vector pos0V(pos0.Data(), dof * dim);
2969  Array<int> pos_dofs;
2970  lim_nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
2971  lim_nodes0->GetSubVector(pos_dofs, pos0V);
2972  if (lim_dist)
2973  {
2974  lim_dist->GetValues(T.ElementNo, ir, d_vals);
2975  }
2976  else
2977  {
2978  d_vals.SetSize(nqp); d_vals = 1.0;
2979  }
2980  }
2981 
2982  // Define ref->physical transformation, when a Coefficient is specified.
2983  IsoparametricTransformation *Tpr = NULL;
2984  if (metric_coeff || lim_coeff || adapt_lim_gf || surf_fit_gf)
2985  {
2986  Tpr = new IsoparametricTransformation;
2987  Tpr->SetFE(&el);
2988  Tpr->ElementNo = T.ElementNo;
2989  Tpr->ElementType = ElementTransformation::ELEMENT;
2990  Tpr->Attribute = T.Attribute;
2991  Tpr->mesh = T.mesh;
2992  Tpr->GetPointMat().Transpose(PMatI);
2993  }
2994 
2995  for (int q = 0; q < nqp; q++)
2996  {
2997  const IntegrationPoint &ip = ir.IntPoint(q);
2998  const DenseMatrix &Jtr_q = Jtr(q);
2999  metric->SetTargetJacobian(Jtr_q);
3000  CalcInverse(Jtr_q, Jrt);
3001  weights(q) = ip.weight * Jtr_q.Det();
3002  double weight_m = weights(q) * metric_normal;
3003 
3004  el.CalcDShape(ip, DSh);
3005  Mult(DSh, Jrt, DS);
3006  MultAtB(PMatI, DS, Jpt);
3007 
3008  if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
3009 
3010  metric->AssembleH(Jpt, DS, weight_m, elmat);
3011 
3012  // TODO: derivatives of adaptivity-based targets.
3013 
3014  if (lim_coeff)
3015  {
3016  el.CalcShape(ip, shape);
3017  PMatI.MultTranspose(shape, p);
3018  pos0.MultTranspose(shape, p0);
3019  weight_m = weights(q) * lim_normal * lim_coeff->Eval(*Tpr, ip);
3020  lim_func->Eval_d2(p, p0, d_vals(q), hess);
3021  for (int i = 0; i < dof; i++)
3022  {
3023  const double w_shape_i = weight_m * shape(i);
3024  for (int j = 0; j < dof; j++)
3025  {
3026  const double w = w_shape_i * shape(j);
3027  for (int d1 = 0; d1 < dim; d1++)
3028  {
3029  for (int d2 = 0; d2 < dim; d2++)
3030  {
3031  elmat(d1*dof + i, d2*dof + j) += w * hess(d1, d2);
3032  }
3033  }
3034  }
3035  }
3036  }
3037  }
3038 
3039  if (adapt_lim_gf) { AssembleElemGradAdaptLim(el, *Tpr, ir, weights, elmat); }
3040  if (surf_fit_gf) { AssembleElemGradSurfFit(el, *Tpr, elmat); }
3041 
3042  delete Tpr;
3043 }
3044 
3045 void TMOP_Integrator::AssembleElemVecAdaptLim(const FiniteElement &el,
3046  IsoparametricTransformation &Tpr,
3047  const IntegrationRule &ir,
3048  const Vector &weights,
3049  DenseMatrix &mat)
3050 {
3051  const int dof = el.GetDof(), dim = el.GetDim(), nqp = weights.Size();
3052  Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3053 
3054  Array<int> dofs;
3055  adapt_lim_gf->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
3056  adapt_lim_gf->GetSubVector(dofs, adapt_lim_gf_e);
3057  adapt_lim_gf->GetValues(Tpr.ElementNo, ir, adapt_lim_gf_q);
3058  adapt_lim_gf0->GetValues(Tpr.ElementNo, ir, adapt_lim_gf0_q);
3059 
3060  // Project the gradient of adapt_lim_gf in the same space.
3061  // The FE coefficients of the gradient go in adapt_lim_gf_grad_e.
3062  DenseMatrix adapt_lim_gf_grad_e(dof, dim);
3063  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3064  el.ProjectGrad(el, Tpr, grad_phys);
3065  Vector grad_ptr(adapt_lim_gf_grad_e.GetData(), dof*dim);
3066  grad_phys.Mult(adapt_lim_gf_e, grad_ptr);
3067 
3068  Vector adapt_lim_gf_grad_q(dim);
3069 
3070  for (int q = 0; q < nqp; q++)
3071  {
3072  const IntegrationPoint &ip = ir.IntPoint(q);
3073  el.CalcShape(ip, shape);
3074  adapt_lim_gf_grad_e.MultTranspose(shape, adapt_lim_gf_grad_q);
3075  adapt_lim_gf_grad_q *= 2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q));
3076  adapt_lim_gf_grad_q *= weights(q) * lim_normal * adapt_lim_coeff->Eval(Tpr, ip);
3077  AddMultVWt(shape, adapt_lim_gf_grad_q, mat);
3078  }
3079 }
3080 
3081 void TMOP_Integrator::AssembleElemGradAdaptLim(const FiniteElement &el,
3082  IsoparametricTransformation &Tpr,
3083  const IntegrationRule &ir,
3084  const Vector &weights,
3085  DenseMatrix &mat)
3086 {
3087  const int dof = el.GetDof(), dim = el.GetDim(), nqp = weights.Size();
3088  Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3089 
3090  Array<int> dofs;
3091  adapt_lim_gf->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
3092  adapt_lim_gf->GetSubVector(dofs, adapt_lim_gf_e);
3093  adapt_lim_gf->GetValues(Tpr.ElementNo, ir, adapt_lim_gf_q);
3094  adapt_lim_gf0->GetValues(Tpr.ElementNo, ir, adapt_lim_gf0_q);
3095 
3096  // Project the gradient of adapt_lim_gf in the same space.
3097  // The FE coefficients of the gradient go in adapt_lim_gf_grad_e.
3098  DenseMatrix adapt_lim_gf_grad_e(dof, dim);
3099  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3100  el.ProjectGrad(el, Tpr, grad_phys);
3101  Vector grad_ptr(adapt_lim_gf_grad_e.GetData(), dof*dim);
3102  grad_phys.Mult(adapt_lim_gf_e, grad_ptr);
3103 
3104  // Project the gradient of each gradient of adapt_lim_gf in the same space.
3105  // The FE coefficients of the second derivatives go in adapt_lim_gf_hess_e.
3106  DenseMatrix adapt_lim_gf_hess_e(dof*dim, dim);
3107  Mult(grad_phys, adapt_lim_gf_grad_e, adapt_lim_gf_hess_e);
3108  // Reshape to be more convenient later (no change in the data).
3109  adapt_lim_gf_hess_e.SetSize(dof, dim*dim);
3110 
3111  Vector adapt_lim_gf_grad_q(dim);
3112  DenseMatrix adapt_lim_gf_hess_q(dim, dim);
3113 
3114  for (int q = 0; q < nqp; q++)
3115  {
3116  const IntegrationPoint &ip = ir.IntPoint(q);
3117  el.CalcShape(ip, shape);
3118 
3119  adapt_lim_gf_grad_e.MultTranspose(shape, adapt_lim_gf_grad_q);
3120  Vector gg_ptr(adapt_lim_gf_hess_q.GetData(), dim*dim);
3121  adapt_lim_gf_hess_e.MultTranspose(shape, gg_ptr);
3122 
3123  const double w = weights(q) * lim_normal * adapt_lim_coeff->Eval(Tpr, ip);
3124  for (int i = 0; i < dof * dim; i++)
3125  {
3126  const int idof = i % dof, idim = i / dof;
3127  for (int j = 0; j <= i; j++)
3128  {
3129  const int jdof = j % dof, jdim = j / dof;
3130  const double entry =
3131  w * ( 2.0 * adapt_lim_gf_grad_q(idim) * shape(idof) *
3132  /* */ adapt_lim_gf_grad_q(jdim) * shape(jdof) +
3133  2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q)) *
3134  adapt_lim_gf_hess_q(idim, jdim) * shape(idof) * shape(jdof));
3135  mat(i, j) += entry;
3136  if (i != j) { mat(j, i) += entry; }
3137  }
3138  }
3139  }
3140 }
3141 
3142 void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x,
3143  IsoparametricTransformation &Tpr,
3144  DenseMatrix &mat)
3145 {
3146  const int el_id = Tpr.ElementNo;
3147  // Check if the element has any DOFs marked for surface fitting.
3148  Array<int> dofs;
3149  surf_fit_gf->FESpace()->GetElementDofs(el_id, dofs);
3150  int count = 0;
3151  for (int s = 0; s < dofs.Size(); s++)
3152  {
3153  count += ((*surf_fit_marker)[dofs[s]]) ? 1 : 0;
3154  }
3155  if (count == 0) { return; }
3156 
3157  const FiniteElement &el_s = *surf_fit_gf->FESpace()->GetFE(el_id);
3158 
3159  const int dof_x = el_x.GetDof(), dim = el_x.GetDim(),
3160  dof_s = el_s.GetDof();
3161 
3162  Vector sigma_e;
3163  surf_fit_gf->GetSubVector(dofs, sigma_e);
3164 
3165  // Project the gradient of sigma in the same space.
3166  // The FE coefficients of the gradient go in surf_fit_grad_e.
3167  DenseMatrix surf_fit_grad_e(dof_s, dim);
3168  Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
3169  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3170  el_s.ProjectGrad(el_s, Tpr, grad_phys);
3171  grad_phys.Mult(sigma_e, grad_ptr);
3172 
3173  Vector shape_x(dof_x), shape_s(dof_s);
3174  const IntegrationRule &ir = el_s.GetNodes();
3175  Vector surf_fit_grad_s(dim);
3176  surf_fit_grad_s = 0.0;
3177 
3178  for (int s = 0; s < dof_s; s++)
3179  {
3180  if ((*surf_fit_marker)[dofs[s]] == false) { continue; }
3181 
3182  const IntegrationPoint &ip = ir.IntPoint(s);
3183  Tpr.SetIntPoint(&ip);
3184  el_x.CalcShape(ip, shape_x);
3185  el_s.CalcShape(ip, shape_s);
3186 
3187  // Note that this gradient is already in physical space.
3188  surf_fit_grad_e.MultTranspose(shape_s, surf_fit_grad_s);
3189  surf_fit_grad_s *= 2.0 * surf_fit_normal *
3190  surf_fit_coeff->Eval(Tpr, ip) * sigma_e(s);
3191 
3192  AddMultVWt(shape_x, surf_fit_grad_s, mat);
3193  }
3194 }
3195 
3196 void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
3197  IsoparametricTransformation &Tpr,
3198  DenseMatrix &mat)
3199 {
3200  const int el_id = Tpr.ElementNo;
3201  // Check if the element has any DOFs marked for surface fitting.
3202  Array<int> dofs;
3203  surf_fit_gf->FESpace()->GetElementDofs(el_id, dofs);
3204  int count = 0;
3205  for (int s = 0; s < dofs.Size(); s++)
3206  {
3207  count += ((*surf_fit_marker)[dofs[s]]) ? 1 : 0;
3208  }
3209  if (count == 0) { return; }
3210 
3211  const FiniteElement &el_s = *surf_fit_gf->FESpace()->GetFE(el_id);
3212 
3213  const int dof_x = el_x.GetDof(), dim = el_x.GetDim(),
3214  dof_s = el_s.GetDof();
3215 
3216  Vector sigma_e;
3217  surf_fit_gf->GetSubVector(dofs, sigma_e);
3218 
3219  DenseMatrix surf_fit_grad_e(dof_s, dim);
3220  Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
3221  DenseMatrix grad_phys;
3222  el_s.ProjectGrad(el_s, Tpr, grad_phys);
3223  grad_phys.Mult(sigma_e, grad_ptr);
3224 
3225  DenseMatrix surf_fit_hess_e(dof_s, dim*dim);
3226  Vector hess_ptr(surf_fit_hess_e.GetData(), dof_s*dim*dim);
3227  surf_fit_hess_e.SetSize(dof_s*dim, dim);
3228  Mult(grad_phys, surf_fit_grad_e, surf_fit_hess_e);
3229  surf_fit_hess_e.SetSize(dof_s, dim * dim);
3230 
3231  const IntegrationRule &ir = el_s.GetNodes();
3232  Vector shape_x(dof_x), shape_s(dof_s);
3233 
3234  Vector surf_fit_grad_s(dim);
3235  DenseMatrix surf_fit_hess_s(dim, dim);
3236 
3237 
3238  for (int s = 0; s < dof_s; s++)
3239  {
3240  if ((*surf_fit_marker)[dofs[s]] == false) { continue; }
3241 
3242  const IntegrationPoint &ip = ir.IntPoint(s);
3243  Tpr.SetIntPoint(&ip);
3244  el_x.CalcShape(ip, shape_x);
3245  el_s.CalcShape(ip, shape_s);
3246 
3247  // These are the sums over k at the dof s (looking at the notes).
3248  surf_fit_grad_e.MultTranspose(shape_s, surf_fit_grad_s);
3249  Vector gg_ptr(surf_fit_hess_s.GetData(), dim * dim);
3250  surf_fit_hess_e.MultTranspose(shape_s, gg_ptr);
3251 
3252  // Loops over the local matrix.
3253  const double w = surf_fit_normal * surf_fit_coeff->Eval(Tpr, ip);
3254  for (int i = 0; i < dof_x * dim; i++)
3255  {
3256  const int idof = i % dof_x, idim = i / dof_x;
3257  for (int j = 0; j <= i; j++)
3258  {
3259  const int jdof = j % dof_x, jdim = j / dof_x;
3260  const double entry =
3261  w * ( 2.0 * surf_fit_grad_s(idim) * shape_x(idof) *
3262  /* */ surf_fit_grad_s(jdim) * shape_x(jdof) +
3263  2.0 * sigma_e(s) * surf_fit_hess_s(idim, jdim) *
3264  /* */ shape_x(idof) * shape_x(jdof));
3265  mat(i, j) += entry;
3266  if (i != j) { mat(j, i) += entry; }
3267  }
3268  }
3269  }
3270 }
3271 
3272 double TMOP_Integrator::GetFDDerivative(const FiniteElement &el,
3273  ElementTransformation &T,
3274  Vector &elfun, const int dofidx,
3275  const int dir, const double e_fx,
3276  bool update_stored)
3277 {
3278  int dof = el.GetDof();
3279  int idx = dir*dof+dofidx;
3280  elfun[idx] += dx;
3281  double e_fxph = GetElementEnergy(el, T, elfun);
3282  elfun[idx] -= dx;
3283  double dfdx = (e_fxph-e_fx)/dx;
3284 
3285  if (update_stored)
3286  {
3287  (*(ElemPertEnergy[T.ElementNo]))(idx) = e_fxph;
3288  (*(ElemDer[T.ElementNo]))(idx) = dfdx;
3289  }
3290 
3291  return dfdx;
3292 }
3293 
3294 void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
3295  ElementTransformation &T,
3296  const Vector &elfun,
3297  Vector &elvect)
3298 {
3299  const int dof = el.GetDof(), dim = el.GetDim(), elnum = T.ElementNo;
3300  if (elnum>=ElemDer.Size())
3301  {
3302  ElemDer.Append(new Vector);
3303  ElemPertEnergy.Append(new Vector);
3304  ElemDer[elnum]->SetSize(dof*dim);
3305  ElemPertEnergy[elnum]->SetSize(dof*dim);
3306  }
3307 
3308  elvect.SetSize(dof*dim);
3309  Vector elfunmod(elfun);
3310 
3311  // In GetElementEnergy(), skip terms that have exact derivative calculations.
3312  fd_call_flag = true;
3313 
3314  // Energy for unperturbed configuration.
3315  const double e_fx = GetElementEnergy(el, T, elfun);
3316 
3317  for (int j = 0; j < dim; j++)
3318  {
3319  for (int i = 0; i < dof; i++)
3320  {
3321  if (discr_tc)
3322  {
3323  discr_tc->UpdateTargetSpecificationAtNode(
3324  el, T, i, j, discr_tc->GetTspecPert1H());
3325  }
3326  elvect(j*dof+i) = GetFDDerivative(el, T, elfunmod, i, j, e_fx, true);
3327  if (discr_tc) { discr_tc->RestoreTargetSpecificationAtNode(T, i); }
3328  }
3329  }
3330  fd_call_flag = false;
3331 
3332  // Contributions from adaptive limiting, surface fitting (exact derivatives).
3333  if (adapt_lim_gf || surf_fit_gf)
3334  {
3335  const IntegrationRule &ir = ActionIntegrationRule(el);
3336  const int nqp = ir.GetNPoints();
3337  DenseTensor Jtr(dim, dim, nqp);
3338  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3339 
3340  IsoparametricTransformation Tpr;
3341  Tpr.SetFE(&el);
3342  Tpr.ElementNo = T.ElementNo;
3343  Tpr.Attribute = T.Attribute;
3344  Tpr.mesh = T.mesh;
3345  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3346  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3347 
3348  Vector weights(nqp);
3349  for (int q = 0; q < nqp; q++)
3350  {
3351  weights(q) = ir.IntPoint(q).weight * Jtr(q).Det();
3352  }
3353 
3354  PMatO.UseExternalData(elvect.GetData(), dof, dim);
3355  if (adapt_lim_gf) { AssembleElemVecAdaptLim(el, Tpr, ir, weights, PMatO); }
3356  if (surf_fit_gf) { AssembleElemVecSurfFit(el, Tpr, PMatO); }
3357  }
3358 }
3359 
3360 void TMOP_Integrator::AssembleElementGradFD(const FiniteElement &el,
3361  ElementTransformation &T,
3362  const Vector &elfun,
3363  DenseMatrix &elmat)
3364 {
3365  const int dof = el.GetDof(), dim = el.GetDim();
3366 
3367  elmat.SetSize(dof*dim);
3368  Vector elfunmod(elfun);
3369 
3370  const Vector &ElemDerLoc = *(ElemDer[T.ElementNo]);
3371  const Vector &ElemPertLoc = *(ElemPertEnergy[T.ElementNo]);
3372 
3373  // In GetElementEnergy(), skip terms that have exact derivative calculations.
3374  fd_call_flag = true;
3375  for (int i = 0; i < dof; i++)
3376  {
3377  for (int j = 0; j < i+1; j++)
3378  {
3379  for (int k1 = 0; k1 < dim; k1++)
3380  {
3381  for (int k2 = 0; k2 < dim; k2++)
3382  {
3383  elfunmod(k2*dof+j) += dx;
3384 
3385  if (discr_tc)
3386  {
3387  discr_tc->UpdateTargetSpecificationAtNode(
3388  el, T, j, k2, discr_tc->GetTspecPert1H());
3389  if (j != i)
3390  {
3391  discr_tc->UpdateTargetSpecificationAtNode(
3392  el, T, i, k1, discr_tc->GetTspecPert1H());
3393  }
3394  else // j==i
3395  {
3396  if (k1 != k2)
3397  {
3398  int idx = k1+k2-1;
3399  discr_tc->UpdateTargetSpecificationAtNode(
3400  el, T, i, idx, discr_tc->GetTspecPertMixH());
3401  }
3402  else // j==i && k1==k2
3403  {
3404  discr_tc->UpdateTargetSpecificationAtNode(
3405  el, T, i, k1, discr_tc->GetTspecPert2H());
3406  }
3407  }
3408  }
3409 
3410  double e_fx = ElemPertLoc(k2*dof+j);
3411  double e_fpxph = GetFDDerivative(el, T, elfunmod, i, k1, e_fx,
3412  false);
3413  elfunmod(k2*dof+j) -= dx;
3414  double e_fpx = ElemDerLoc(k1*dof+i);
3415 
3416  elmat(k1*dof+i, k2*dof+j) = (e_fpxph - e_fpx) / dx;
3417  elmat(k2*dof+j, k1*dof+i) = (e_fpxph - e_fpx) / dx;
3418 
3419  if (discr_tc)
3420  {
3421  discr_tc->RestoreTargetSpecificationAtNode(T, i);
3422  discr_tc->RestoreTargetSpecificationAtNode(T, j);
3423  }
3424  }
3425  }
3426  }
3427  }
3428  fd_call_flag = false;
3429 
3430  // Contributions from adaptive limiting.
3431  if (adapt_lim_gf || surf_fit_gf)
3432  {
3433  const IntegrationRule &ir = GradientIntegrationRule(el);
3434  const int nqp = ir.GetNPoints();
3435  DenseTensor Jtr(dim, dim, nqp);
3436  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3437 
3438  IsoparametricTransformation Tpr;
3439  Tpr.SetFE(&el);
3440  Tpr.ElementNo = T.ElementNo;
3441  Tpr.Attribute = T.Attribute;
3442  Tpr.mesh = T.mesh;
3443  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3444  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3445 
3446  Vector weights(nqp);
3447  for (int q = 0; q < nqp; q++)
3448  {
3449  weights(q) = ir.IntPoint(q).weight * Jtr(q).Det();
3450  }
3451 
3452  if (adapt_lim_gf) { AssembleElemGradAdaptLim(el, Tpr, ir, weights, elmat); }
3453  if (surf_fit_gf) { AssembleElemGradSurfFit(el, Tpr, elmat); }
3454  }
3455 }
3456 
3457 void TMOP_Integrator::UpdateSurfaceFittingWeight(double factor)
3458 {
3459  if (!surf_fit_coeff) { return; }
3460 
3461  if (surf_fit_coeff)
3462  {
3463  auto cf = dynamic_cast<ConstantCoefficient *>(surf_fit_coeff);
3464  MFEM_VERIFY(cf, "Dynamic weight works only with a ConstantCoefficient.");
3465  cf->constant *= factor;
3466  }
3467 }
3468 
3469 double TMOP_Integrator::GetSurfaceFittingWeight()
3470 {
3471  if (surf_fit_coeff)
3472  {
3473  auto cf = dynamic_cast<ConstantCoefficient *>(surf_fit_coeff);
3474  MFEM_VERIFY(cf, "Dynamic weight works only with a ConstantCoefficient.");
3475  return cf->constant;
3476  }
3477  return 0.0;
3478 }
3479 
3480 void TMOP_Integrator::EnableNormalization(const GridFunction &x)
3481 {
3482  ComputeNormalizationEnergies(x, metric_normal, lim_normal, surf_fit_normal);
3483  metric_normal = 1.0 / metric_normal;
3484  lim_normal = 1.0 / lim_normal;
3485  //if (surf_fit_gf) { surf_fit_normal = 1.0 / surf_fit_normal; }
3486  if (surf_fit_gf) { surf_fit_normal = lim_normal; }
3487 }
3488 
3489 #ifdef MFEM_USE_MPI
3490 void TMOP_Integrator::ParEnableNormalization(const ParGridFunction &x)
3491 {
3492  double loc[3];
3493  ComputeNormalizationEnergies(x, loc[0], loc[1], loc[2]);
3494  double rdc[3];
3495  MPI_Allreduce(loc, rdc, 3, MPI_DOUBLE, MPI_SUM, x.ParFESpace()->GetComm());
3496  metric_normal = 1.0 / rdc[0];
3497  lim_normal = 1.0 / rdc[1];
3498  // if (surf_fit_gf) { surf_fit_normal = 1.0 / rdc[2]; }
3499  if (surf_fit_gf) { surf_fit_normal = lim_normal; }
3500 }
3501 #endif
3502 
3503 void TMOP_Integrator::ComputeNormalizationEnergies(const GridFunction &x,
3504  double &metric_energy,
3505  double &lim_energy,
3506  double &surf_fit_gf_energy)
3507 {
3508  Array<int> vdofs;
3509  Vector x_vals;
3510  const FiniteElementSpace* const fes = x.FESpace();
3511 
3512  const int dim = fes->GetMesh()->Dimension();
3513  Jrt.SetSize(dim);
3514  Jpr.SetSize(dim);
3515  Jpt.SetSize(dim);
3516 
3517  metric_energy = 0.0;
3518  lim_energy = 0.0;
3519  surf_fit_gf_energy = 0.0;
3520  for (int i = 0; i < fes->GetNE(); i++)
3521  {
3522  const FiniteElement *fe = fes->GetFE(i);
3523  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
3524  const int nqp = ir.GetNPoints();
3525  DenseTensor Jtr(dim, dim, nqp);
3526  const int dof = fe->GetDof();
3527  DSh.SetSize(dof, dim);
3528 
3529  fes->GetElementVDofs(i, vdofs);
3530  x.GetSubVector(vdofs, x_vals);
3531  PMatI.UseExternalData(x_vals.GetData(), dof, dim);
3532 
3533  targetC->ComputeElementTargets(i, *fe, ir, x_vals, Jtr);
3534 
3535  for (int q = 0; q < nqp; q++)
3536  {
3537  const IntegrationPoint &ip = ir.IntPoint(q);
3538  metric->SetTargetJacobian(Jtr(q));
3539  CalcInverse(Jtr(q), Jrt);
3540  const double weight = ip.weight * Jtr(q).Det();
3541 
3542  fe->CalcDShape(ip, DSh);
3543  MultAtB(PMatI, DSh, Jpr);
3544  Mult(Jpr, Jrt, Jpt);
3545 
3546  metric_energy += weight * metric->EvalW(Jpt);
3547  lim_energy += weight;
3548  }
3549 
3550  // Normalization of the surface fitting term.
3551  if (surf_fit_gf)
3552  {
3553  Array<int> dofs;
3554  Vector sigma_e;
3555  surf_fit_gf->FESpace()->GetElementDofs(i, dofs);
3556  surf_fit_gf->GetSubVector(dofs, sigma_e);
3557  for (int s = 0; s < dofs.Size(); s++)
3558  {
3559  if ((*surf_fit_marker)[dofs[s]] == true)
3560  {
3561  surf_fit_gf_energy += sigma_e(s) * sigma_e(s);
3562  }
3563  }
3564  }
3565  }
3566 
3567  if (targetC->ContainsVolumeInfo() == false)
3568  {
3569  // Special case when the targets don't contain volumetric information.
3570  lim_energy = fes->GetNE();
3571  }
3572 }
3573 
3574 void TMOP_Integrator::ComputeMinJac(const Vector &x,
3575  const FiniteElementSpace &fes)
3576 {
3577  const FiniteElement *fe = fes.GetFE(0);
3578  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
3579  const int NE = fes.GetMesh()->GetNE(), dim = fe->GetDim(),
3580  dof = fe->GetDof(), nsp = ir.GetNPoints();
3581 
3582  Array<int> xdofs(dof * dim);
3583  DenseMatrix dshape(dof, dim), pos(dof, dim);
3584  Vector posV(pos.Data(), dof * dim);
3585  Jpr.SetSize(dim);
3586 
3587  dx = std::numeric_limits<float>::max();
3588 
3589  double detv_sum;
3590  double detv_avg_min = std::numeric_limits<float>::max();
3591  for (int i = 0; i < NE; i++)
3592  {
3593  fes.GetElementVDofs(i, xdofs);
3594  x.GetSubVector(xdofs, posV);
3595  detv_sum = 0.;
3596  for (int j = 0; j < nsp; j++)
3597  {
3598  fes.GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
3599  MultAtB(pos, dshape, Jpr);
3600  detv_sum += std::fabs(Jpr.Det());
3601  }
3602  double detv_avg = pow(detv_sum/nsp, 1./dim);
3603  detv_avg_min = std::min(detv_avg, detv_avg_min);
3604  }
3605  dx = detv_avg_min / dxscale;
3606 }
3607 
3608 void TMOP_Integrator::UpdateAfterMeshPositionChange(const Vector &new_x)
3609 {
3610  if (discr_tc)
3611  {
3612  PA.Jtr_needs_update = true;
3613  }
3614  // Update adapt_lim_gf if adaptive limiting is enabled.
3615  if (adapt_lim_gf) { adapt_lim_eval->ComputeAtNewPosition(new_x, *adapt_lim_gf); }
3616 
3617  // Update surf_fit_gf if surface fitting is enabled.
3618  if (surf_fit_gf)
3619  {
3620  surf_fit_eval->ComputeAtNewPosition(new_x, *surf_fit_gf);
3621  }
3622 }
3623 
3624 void TMOP_Integrator::ComputeFDh(const Vector &x, const FiniteElementSpace &fes)
3625 {
3626  if (!fdflag) { return; }
3627  ComputeMinJac(x, fes);
3628 }
3629 
3630 #ifdef MFEM_USE_MPI
3631 void TMOP_Integrator::ComputeFDh(const Vector &x,
3632  const ParFiniteElementSpace &pfes)
3633 {
3634  if (!fdflag) { return; }
3635  ComputeMinJac(x, pfes);
3636  double min_jac_all;
3637  MPI_Allreduce(&dx, &min_jac_all, 1, MPI_DOUBLE, MPI_MIN, pfes.GetComm());
3638  dx = min_jac_all;
3639 }
3640 #endif
3641 
3642 void TMOP_Integrator::EnableFiniteDifferences(const GridFunction &x)
3643 {
3644  fdflag = true;
3645  const FiniteElementSpace *fes = x.FESpace();
3646  ComputeFDh(x,*fes);
3647  if (discr_tc)
3648  {
3649  discr_tc->UpdateTargetSpecification(x);
3650  discr_tc->UpdateGradientTargetSpecification(x, dx);
3651  discr_tc->UpdateHessianTargetSpecification(x, dx);
3652  }
3653 }
3654 
3655 #ifdef MFEM_USE_MPI
3656 void TMOP_Integrator::EnableFiniteDifferences(const ParGridFunction &x)
3657 {
3658  fdflag = true;
3659  const ParFiniteElementSpace *pfes = x.ParFESpace();
3660  ComputeFDh(x,*pfes);
3661  if (discr_tc)
3662  {
3663  discr_tc->UpdateTargetSpecification(x);
3664  discr_tc->UpdateGradientTargetSpecification(x, dx);
3665  discr_tc->UpdateHessianTargetSpecification(x, dx);
3666  }
3667 }
3668 #endif
3669 
3670 void TMOPComboIntegrator::EnableLimiting(const GridFunction &n0,
3671  const GridFunction &dist,
3672  Coefficient &w0,
3673  TMOP_LimiterFunction *lfunc)
3674 {
3675  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3676 
3677  tmopi[0]->EnableLimiting(n0, dist, w0, lfunc);
3678  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
3679 }
3680 
3681 void TMOPComboIntegrator::EnableLimiting(const GridFunction &n0,
3682  Coefficient &w0,
3683  TMOP_LimiterFunction *lfunc)
3684 {
3685  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3686 
3687  tmopi[0]->EnableLimiting(n0, w0, lfunc);
3688  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
3689 }
3690 
3691 void TMOPComboIntegrator::SetLimitingNodes(const GridFunction &n0)
3692 {
3693  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3694 
3695  tmopi[0]->SetLimitingNodes(n0);
3696  for (int i = 1; i < tmopi.Size(); i++) { tmopi[i]->DisableLimiting(); }
3697 }
3698 
3699 double TMOPComboIntegrator::GetElementEnergy(const FiniteElement &el,
3700  ElementTransformation &T,
3701  const Vector &elfun)
3702 {
3703  double energy= 0.0;
3704  for (int i = 0; i < tmopi.Size(); i++)
3705  {
3706  energy += tmopi[i]->GetElementEnergy(el, T, elfun);
3707  }
3708  return energy;
3709 }
3710 
3711 void TMOPComboIntegrator::AssembleElementVector(const FiniteElement &el,
3712  ElementTransformation &T,
3713  const Vector &elfun,
3714  Vector &elvect)
3715 {
3716  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3717 
3718  tmopi[0]->AssembleElementVector(el, T, elfun, elvect);
3719  for (int i = 1; i < tmopi.Size(); i++)
3720  {
3721  Vector elvect_i;
3722  tmopi[i]->AssembleElementVector(el, T, elfun, elvect_i);
3723  elvect += elvect_i;
3724  }
3725 }
3726 
3727 void TMOPComboIntegrator::AssembleElementGrad(const FiniteElement &el,
3728  ElementTransformation &T,
3729  const Vector &elfun,
3730  DenseMatrix &elmat)
3731 {
3732  MFEM_VERIFY(tmopi.Size() > 0, "No TMOP_Integrators were added.");
3733 
3734  tmopi[0]->AssembleElementGrad(el, T, elfun, elmat);
3735  for (int i = 1; i < tmopi.Size(); i++)
3736  {
3737  DenseMatrix elmat_i;
3738  tmopi[i]->AssembleElementGrad(el, T, elfun, elmat_i);
3739  elmat += elmat_i;
3740  }
3741 }
3742 
3743 double TMOPComboIntegrator::GetRefinementElementEnergy(const FiniteElement &el,
3744  ElementTransformation &T,
3745  const Vector &elfun,
3746  const IntegrationRule &irule)
3747 {
3748  double energy= 0.0;
3749  for (int i = 0; i < tmopi.Size(); i++)
3750  {
3751  energy += tmopi[i]->GetRefinementElementEnergy(el, T, elfun, irule);
3752  }
3753  return energy;
3754 }
3755 
3756 double TMOPComboIntegrator::GetDerefinementElementEnergy(
3757  const FiniteElement &el,
3758  ElementTransformation &T,
3759  const Vector &elfun)
3760 {
3761  double energy= 0.0;
3762  for (int i = 0; i < tmopi.Size(); i++)
3763  {
3764  energy += tmopi[i]->GetDerefinementElementEnergy(el, T, elfun);
3765  }
3766  return energy;
3767 }
3768 
3769 void TMOPComboIntegrator::EnableNormalization(const GridFunction &x)
3770 {
3771  const int cnt = tmopi.Size();
3772  double total_integral = 0.0;
3773  for (int i = 0; i < cnt; i++)
3774  {
3775  tmopi[i]->EnableNormalization(x);
3776  total_integral += 1.0 / tmopi[i]->metric_normal;
3777  }
3778  for (int i = 0; i < cnt; i++)
3779  {
3780  tmopi[i]->metric_normal = 1.0 / total_integral;
3781  }
3782 }
3783 
3784 #ifdef MFEM_USE_MPI
3785 void TMOPComboIntegrator::ParEnableNormalization(const ParGridFunction &x)
3786 {
3787  const int cnt = tmopi.Size();
3788  double total_integral = 0.0;
3789  for (int i = 0; i < cnt; i++)
3790  {
3791  tmopi[i]->ParEnableNormalization(x);
3792  total_integral += 1.0 / tmopi[i]->metric_normal;
3793  }
3794  for (int i = 0; i < cnt; i++)
3795  {
3796  tmopi[i]->metric_normal = 1.0 / total_integral;
3797  }
3798 }
3799 #endif
3800 
3801 void TMOPComboIntegrator::AssemblePA(const FiniteElementSpace &fes)
3802 {
3803  for (int i = 0; i < tmopi.Size(); i++)
3804  {
3805  tmopi[i]->AssemblePA(fes);
3806  }
3807 }
3808 
3809 void TMOPComboIntegrator::AssembleGradPA(const Vector &xe,
3810  const FiniteElementSpace &fes)
3811 {
3812  for (int i = 0; i < tmopi.Size(); i++)
3813  {
3814  tmopi[i]->AssembleGradPA(xe,fes);
3815  }
3816 }
3817 
3818 void TMOPComboIntegrator::AssembleGradDiagonalPA(Vector &de) const
3819 {
3820  for (int i = 0; i < tmopi.Size(); i++)
3821  {
3822  tmopi[i]->AssembleGradDiagonalPA(de);
3823  }
3824 }
3825 
3826 void TMOPComboIntegrator::AddMultPA(const Vector &xe, Vector &ye) const
3827 {
3828  for (int i = 0; i < tmopi.Size(); i++)
3829  {
3830  tmopi[i]->AddMultPA(xe, ye);
3831  }
3832 }
3833 
3834 void TMOPComboIntegrator::AddMultGradPA(const Vector &re, Vector &ce) const
3835 {
3836  for (int i = 0; i < tmopi.Size(); i++)
3837  {
3838  tmopi[i]->AddMultGradPA(re, ce);
3839  }
3840 }
3841 
3842 double TMOPComboIntegrator::GetLocalStateEnergyPA(const Vector &xe) const
3843 {
3844  double energy = 0.0;
3845  for (int i = 0; i < tmopi.Size(); i++)
3846  {
3847  energy += tmopi[i]->GetLocalStateEnergyPA(xe);
3848  }
3849  return energy;
3850 }
3851 
3852 void InterpolateTMOP_QualityMetric(TMOP_QualityMetric &metric,
3853  const TargetConstructor &tc,
3854  const Mesh &mesh, GridFunction &metric_gf)
3855 {
3856  const int NE = mesh.GetNE();
3857  const GridFunction &nodes = *mesh.GetNodes();
3858  const int dim = mesh.Dimension();
3859  DenseMatrix Winv(dim), T(dim), A(dim), dshape, pos;
3860  Array<int> pos_dofs, gf_dofs;
3861  DenseTensor W;
3862  Vector posV;
3863 
3864  for (int i = 0; i < NE; i++)
3865  {
3866  const FiniteElement &fe_pos = *nodes.FESpace()->GetFE(i);
3867  const IntegrationRule &ir = metric_gf.FESpace()->GetFE(i)->GetNodes();
3868  const int nsp = ir.GetNPoints(), dof = fe_pos.GetDof();
3869 
3870  dshape.SetSize(dof, dim);
3871  pos.SetSize(dof, dim);
3872  posV.SetDataAndSize(pos.Data(), dof * dim);
3873 
3874  metric_gf.FESpace()->GetElementDofs(i, gf_dofs);
3875  nodes.FESpace()->GetElementVDofs(i, pos_dofs);
3876  nodes.GetSubVector(pos_dofs, posV);
3877 
3878  W.SetSize(dim, dim, nsp);
3879  tc.ComputeElementTargets(i, fe_pos, ir, posV, W);
3880 
3881  for (int j = 0; j < nsp; j++)
3882  {
3883  const DenseMatrix &Wj = W(j);
3884  metric.SetTargetJacobian(Wj);
3885  CalcInverse(Wj, Winv);
3886 
3887  const IntegrationPoint &ip = ir.IntPoint(j);
3888  fe_pos.CalcDShape(ip, dshape);
3889  MultAtB(pos, dshape, A);
3890  Mult(A, Winv, T);
3891 
3892  metric_gf(gf_dofs[j]) = metric.EvalW(T);
3893  }
3894  }
3895 }
3896 
3897 } // namespace mfem
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
mfem::FiniteElement
Abstract class for all finite elements.
Definition: fe_base.hpp:235
mfem::TMOP_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_Integrator::GetDerefinementElementEnergy
virtual double GetDerefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Definition: tmop.cpp:2714
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:1375
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:138
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_base.hpp:311
mfem::TMOP_QualityMetric::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const =0
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
mfem::AddMultVWt
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:2761
mfem::FiniteElementSpace::GetNDofs
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:560
mfem::TMOP_Integrator::AssembleElemVecSurfFit
void AssembleElemVecSurfFit(const FiniteElement &el_x, IsoparametricTransformation &Tpr, DenseMatrix &mat)
Definition: tmop.cpp:3142
mfem::AdaptivityEvaluator::SetSerialMetaInfo
void SetSerialMetaInfo(const Mesh &m, const FiniteElementCollection &fec, int num_comp)
Definition: tmop.cpp:2269
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:1053
mfem::GridFunction::GetVectorValues
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:703
mfem::TMOP_Integrator::EnableFiniteDifferences
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
Definition: tmop.cpp:3642
mfem::Vector::NewDataAndSize
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:162
mfem::TargetConstructor::target_type
const TargetType target_type
Definition: tmop.hpp:935
mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
mfem::TMOP_Integrator::PMatO
DenseMatrix PMatO
Definition: tmop.hpp:1376
mfem::TMOPComboIntegrator::AddMultPA
virtual void AddMultPA(const Vector &, Vector &) const
Method for partially assembled action.
Definition: tmop.cpp:3826
mfem::TMOP_Integrator::Jtr
DenseTensor Jtr
Definition: tmop.hpp:1406
mfem::DenseTensor::UseExternalData
void UseExternalData(double *ext_data, int i, int j, int k)
Definition: densemat.hpp:847
mfem::TMOP_Metric_022::min_detT
double & min_detT
Definition: tmop.hpp:256
mfem::FiniteElementSpace::Update
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition: fespace.cpp:3336
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:455
mfem::InvariantsEvaluator2D::Assemble_ddI2b
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:384
mfem::IntegrationRules::Get
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
mfem::TMOP_Metric_302::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:674
mfem::TMOP_Integrator::GetSurfaceFittingWeight
double GetSurfaceFittingWeight()
Get the surface fitting weight.
Definition: tmop.cpp:3469
mfem::ParFiniteElementSpace::GetParMesh
ParMesh * GetParMesh() const
Definition: pfespace.hpp:277
mfem::VectorCoefficient
Base class for vector Coefficients that optionally depend on time and space.
Definition: coefficient.hpp:432
mfem::ConstantCoefficient
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
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:1376
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:1328
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetOrientation
virtual void SetSerialDiscreteTargetOrientation(const GridFunction &tspec_)
Definition: tmop.cpp:1474
mfem::TMOPComboIntegrator::GetElementEnergy
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Compute the local energy.
Definition: tmop.cpp:3699
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:5877
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:3480
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:521
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:3852
mfem::TMOP_Integrator::metric
TMOP_QualityMetric * metric
Definition: tmop.hpp:1310
mfem::TMOP_Metric_056::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:308
mfem::DenseMatrix::Det
double Det() const
Definition: densemat.cpp:436
mfem::TMOP_LimiterFunction
Base class for limiting functions to be used in class TMOP_Integrator.
Definition: tmop.hpp:798
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:1832
mfem::Mult
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:472
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:1332
mfem::ParFiniteElementSpace::Update
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3284
mfem::Vector::Norml2
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:792
mfem::DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange
void ParUpdateAfterMeshTopologyChange()
Definition: tmop.cpp:1343
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:1376
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:93
mfem::Vector::GetSubVector
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:534
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:863
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::TMOP_Integrator::GetSurfaceFittingErrors
void GetSurfaceFittingErrors(double &err_avg, double &err_max)
Definition: tmop.cpp:2436
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
mfem::TMOP_Integrator::metric_normal
double metric_normal
Definition: tmop.hpp:1320
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:2365
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:928
mfem::TMOPComboIntegrator::EnableNormalization
void EnableNormalization(const GridFunction &x)
Normalization factor that considers all integrators in the combination.
Definition: tmop.cpp:3769
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_Metric_313::min_detT
double & min_detT
Definition: tmop.hpp:524
mfem::Operator::Mult
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
mfem::TMOP_Metric_085::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
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:3801
mfem::IsoparametricTransformation::GetFE
const FiniteElement * GetFE() const
Get the current element used to compute the transformations.
Definition: eltrans.hpp:389
mfem::TMOP_Metric_009::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:225
mfem::TMOP_Integrator::ir
const IntegrationRule * ir
Definition: tmop.hpp:1414
mfem::TMOP_Metric_058::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:326
mfem::TMOP_Integrator::AssembleElementVectorFD
void AssembleElementVectorFD(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Definition: tmop.cpp:3294
mfem::TargetConstructor::nodes
const GridFunction * nodes
Definition: tmop.hpp:932
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSkew
virtual void SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
Definition: tmop.cpp:1458
mfem::TMOP_Integrator::UpdateSurfaceFittingWeight
void UpdateSurfaceFittingWeight(double factor)
Update the surface fitting weight as surf_fit_coeff *= factor;.
Definition: tmop.cpp:3457
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:2179
mfem::DiscreteAdaptTC::SetRefinementSubElement
void SetRefinementSubElement(int amr_el_)
Definition: tmop.hpp:1291
mfem::InvariantsEvaluator2D::Assemble_TProd
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition: invariants.hpp:432
mfem::InvariantsEvaluator2D::Assemble_ddI1
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:226
mfem::DenseMatrix::GetData
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
mfem::IsoparametricTransformation::OrderJ
virtual int OrderJ() const
Return the order of the elements of the Jacobian of the transformation.
Definition: eltrans.cpp:443
mfem::FiniteElement::ProjectGrad
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:160
mfem::Vector::GetData
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
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:819
mfem::Mesh::GetElementBaseGeometry
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1062
mfem::Vector::UseDevice
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
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:520
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:1617
mfem::TMOP_Integrator::AssembleElementGradExact
void AssembleElementGradExact(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Definition: tmop.cpp:2937
mfem::TMOP_Integrator::ElemDer
Array< Vector * > ElemDer
Definition: tmop.hpp:1363
mfem::DiscreteAdaptTC::GetTspecPertMixH
const Vector & GetTspecPertMixH()
Definition: tmop.hpp:1241
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:3670
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::AnalyticAdaptTC::scalar_tspec
Coefficient * scalar_tspec
Definition: tmop.hpp:1052
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:1466
mfem::Geometry::GetGeomToPerfGeomJac
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:98
mfem::DiscreteAdaptTC::UpdateTargetSpecification
void UpdateTargetSpecification(const Vector &new_x, bool use_flag=false)
Definition: tmop.cpp:1534
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:3834
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_base.hpp:320
mfem::Add
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1916
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:590
mfem::TMOP_QualityMetric
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
mfem::DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec
void FinalizeSerialDiscreteTargetSpec(const GridFunction &tspec_)
Definition: tmop.cpp:1482
mfem::Geometries
Geometry Geometries
Definition: fe.cpp:49
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
mfem::TMOP_Integrator::adapt_lim_coeff
Coefficient * adapt_lim_coeff
Definition: tmop.hpp:1340
mfem::DenseMatrix::MultTranspose
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:188
mfem::AdaptivityEvaluator::~AdaptivityEvaluator
virtual ~AdaptivityEvaluator()
Definition: tmop.cpp:2304
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:2146
mfem::TMOP_Integrator::UpdateAfterMeshPositionChange
void UpdateAfterMeshPositionChange(const Vector &new_x)
Definition: tmop.cpp:3608
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::FiniteElementSpace::GetElementVDofs
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
mfem::TMOPComboIntegrator::SetLimitingNodes
void SetLimitingNodes(const GridFunction &n0)
Update the original/reference nodes used for limiting.
Definition: tmop.cpp:3691
mfem::DiscreteAdaptTC::SetDiscreteTargetBase
void SetDiscreteTargetBase(const GridFunction &tspec_)
Definition: tmop.cpp:1413
mfem::TMOP_Integrator::adapt_lim_eval
AdaptivityEvaluator * adapt_lim_eval
Definition: tmop.hpp:1341
mfem::ParFiniteElementSpace::GetComm
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
mfem::ParFiniteElementSpace::GetFE
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:535
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:2042
mfem::TMOP_Integrator::discr_tc
DiscreteAdaptTC * discr_tc
Definition: tmop.hpp:1350
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:2340
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:685
mfem::DiscreteAdaptTC::UpdateTargetSpecificationAtNode
void UpdateTargetSpecificationAtNode(const FiniteElement &el, ElementTransformation &T, int nodenum, int idir, const Vector &IntData)
Definition: tmop.cpp:1552
mfem::Operator::Height
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
mfem::DiscreteAdaptTC::SetTspecFromIntRule
void SetTspecFromIntRule(int e_id, const IntegrationRule &intrule)
Definition: tmop.cpp:1583
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::ParGridFunction::Update
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:81
mfem::TMOP_Combo_QualityMetric::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:46
mfem::FiniteElementSpace::GetMesh
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
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:529
mfem::DiscreteAdaptTC::ResetRefinementTspecData
void ResetRefinementTspecData()
Definition: tmop.hpp:1276
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:2797
mfem::TMOPComboIntegrator::AssembleGradDiagonalPA
virtual void AssembleGradDiagonalPA(Vector &) const
Method for computing the diagonal of the gradient with partial assembly.
Definition: tmop.cpp:3818
mfem::TMOP_Integrator::targetC
const TargetConstructor * targetC
Definition: tmop.hpp:1311
mfem::TMOP_Integrator::ElemPertEnergy
Array< Vector * > ElemPertEnergy
Definition: tmop.hpp:1364
mfem::TMOP_Integrator::lim_nodes0
const GridFunction * lim_nodes0
Definition: tmop.hpp:1325
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::TMOP_Integrator::adapt_lim_pgf0
const ParGridFunction * adapt_lim_pgf0
Definition: tmop.hpp:1337
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_base.hpp:390
mfem::InvariantsEvaluator3D::SetDerivativeMatrix
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 3, using column-major storage.
Definition: invariants.hpp: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:3809
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:1391
mfem::DenseTensor::GetData
double * GetData(int k)
Definition: densemat.hpp:886
mfem::ElementTransformation::Weight
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
mfem::TMOP_Metric_077::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:491
mfem::DenseTensor::SetSize
void SetSize(int i, int j, int k, MemoryType mt_=MemoryType::PRESERVE)
Definition: densemat.hpp:838
mfem::TMOPComboIntegrator::ParEnableNormalization
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:3785
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:366
mfem::GridFunction::GetValues
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:510
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:2768
mfem::TMOPComboIntegrator::GetLocalStateEnergyPA
virtual double GetLocalStateEnergyPA(const Vector &) const
Compute the local (to the MPI rank) energy with partial assembly.
Definition: tmop.cpp:3842
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::DiscreteAdaptTC::tspec_fesv
FiniteElementSpace * tspec_fesv
Definition: tmop.hpp:1117
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:557
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:121
mfem::TMOP_Integrator::metric_coeff
Coefficient * metric_coeff
Definition: tmop.hpp:1318
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::TMOP_Integrator::AssembleElemGradSurfFit
void AssembleElemGradSurfFit(const FiniteElement &el_x, IsoparametricTransformation &Tpr, DenseMatrix &mat)
Definition: tmop.cpp:3196
mfem::Vector::SetData
void SetData(double *d)
Definition: vector.hpp:149
mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:652
mfem::TargetConstructor::ContainsVolumeInfo
virtual bool ContainsVolumeInfo() const
Checks if the target matrices contain non-trivial size specification.
Definition: tmop.cpp:1128
mfem::TMOP_Integrator::h_metric
TMOP_QualityMetric * h_metric
Definition: tmop.hpp:1309
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:999
mfem::DenseMatrix::GetColumn
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1270
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:442
mfem::DiscreteAdaptTC::RestoreTargetSpecificationAtNode
void RestoreTargetSpecificationAtNode(ElementTransformation &T, int nodenum)
Definition: tmop.cpp:1569
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::FiniteElementSpace::GetElementDofs
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2680
mfem::DiscreteAdaptTC::SetTspecAtIndex
void SetTspecAtIndex(int idx, const GridFunction &tspec_)
Definition: tmop.cpp:1437
mfem::TMOP_Integrator::GetRefinementElementEnergy
virtual double GetRefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, const IntegrationRule &irule)
Computes the mean of the energies of the given element&#39;s children.
Definition: tmop.cpp:2628
mfem::TargetConstructor::ComputeAvgVolume
void ComputeAvgVolume() const
Definition: tmop.cpp:1014
mfem::ad::cos
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
mfem::DenseTensor::SizeI
int SizeI() const
Definition: densemat.hpp:832
mfem::CalcInverse
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2182
mfem::ad::pow
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
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:2328
omega
double omega
Definition: ex25.cpp:141
mfem::TMOP_Integrator::dxscale
double dxscale
Definition: tmop.hpp:1355
mfem::TMOP_Integrator::lim_coeff
Coefficient * lim_coeff
Definition: tmop.hpp:1326
tmop_tools.hpp
mfem::Vector::Min
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1182
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::TMOP_Integrator::surf_fit_eval
AdaptivityEvaluator * surf_fit_eval
Definition: tmop.hpp:1347
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::ad::sin
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
mfem::DiscreteAdaptTC::UpdateHessianTargetSpecification
void UpdateHessianTargetSpecification(const Vector &x, double dx, bool use_flag=false)
Definition: tmop.cpp:2205
mfem::DenseMatrix::Transpose
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1359
mfem::TMOP_Integrator::AssembleElemGradAdaptLim
void AssembleElemGradAdaptLim(const FiniteElement &el, IsoparametricTransformation &Tpr, const IntegrationRule &ir, const Vector &weights, DenseMatrix &m)
Definition: tmop.cpp:3081
mfem::DenseMatrix::Trace
double Trace() const
Trace of a square matrix.
Definition: densemat.cpp:412
mfem::AddMultABt
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2481
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:2782
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
mfem::FiniteElement::GetDof
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
mfem::InvariantsEvaluator3D::Assemble_ddI2
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:948
mfem::TMOP_Metric_211::eps
const double eps
Definition: tmop.hpp:414
mfem::DiscreteAdaptTC::tspec_refine
DenseMatrix tspec_refine
Definition: tmop.hpp:1106
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:686
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:1054
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:1486
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
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:2314
mfem::Operator::NumCols
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition: operator.hpp:75
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::ReleasePADeviceMemory
void ReleasePADeviceMemory(bool copy_to_host=true)
Definition: tmop.cpp:2314
mfem::TargetConstructor::UsesPhysicalCoordinates
bool UsesPhysicalCoordinates() const
Return true if the methods ComputeElementTargets(), ComputeAllElementTargets(), and ComputeElementTar...
Definition: tmop.hpp:1001
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:921
mfem::TMOPComboIntegrator::AssembleElementGrad
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: tmop.cpp:3727
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:1113
mfem::Vector::SetDataAndSize
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
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::GridFunction::Update
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:164
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::DiscreteAdaptTC::tspec_gf
GridFunction * tspec_gf
Definition: tmop.hpp:1119
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_base.hpp:1181
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::UpdateAfterMeshTopologyChange
void UpdateAfterMeshTopologyChange()
Definition: tmop.cpp:2465
mfem::TMOP_Integrator::PMatI
DenseMatrix PMatI
Definition: tmop.hpp:1376
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:833
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:3624
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:454
mfem::TMOP_Integrator::ComputeNormalizationEnergies
void ComputeNormalizationEnergies(const GridFunction &x, double &metric_energy, double &lim_energy, double &surf_fit_gf_energy)
Definition: tmop.cpp:3503
mfem::TMOP_LimiterFunction::Eval
virtual double Eval(const Vector &x, const Vector &x0, double d) const =0
Returns the limiting function, f(x, x0, d).
mfem::TMOP_Integrator::surf_fit_coeff
Coefficient * surf_fit_coeff
Definition: tmop.hpp:1346
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
mfem::DiscreteAdaptTC::adapt_eval
AdaptivityEvaluator * adapt_eval
Definition: tmop.hpp:1136
mfem::DiscreteAdaptTC::ptspec_fesv
ParFiniteElementSpace * ptspec_fesv
Definition: tmop.hpp:1122
mfem::TMOP_Integrator::dx
double dx
Definition: tmop.hpp:1354
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:1413
mfem::TMOPComboIntegrator::tmopi
Array< TMOP_Integrator * > tmopi
Definition: tmop.hpp:1724
mfem::IsoparametricTransformation
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
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:3574
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:348
mfem::TMOP_Integrator::EnergyIntegrationRule
const IntegrationRule & EnergyIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1477
mfem::TMOP_Integrator::GradientIntegrationRule
const IntegrationRule & GradientIntegrationRule(const FiniteElement &el) const
Definition: tmop.hpp:1491
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:3272
mfem::ad::sqrt
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
mfem::TMOP_Integrator::fdflag
bool fdflag
Definition: tmop.hpp:1353
mfem::AdaptivityEvaluator::SetParMetaInfo
void SetParMetaInfo(const ParMesh &m, const FiniteElementCollection &fec, int num_comp)
Parallel version of SetSerialMetaInfo.
Definition: tmop.cpp:2282
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
mfem::TMOP_Integrator::surf_fit_normal
double surf_fit_normal
Definition: tmop.hpp:1348
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:2491
mfem::TMOP_Integrator::AssembleElemVecAdaptLim
void AssembleElemVecAdaptLim(const FiniteElement &el, IsoparametricTransformation &Tpr, const IntegrationRule &ir, const Vector &weights, DenseMatrix &mat)
Definition: tmop.cpp:3045
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:381
mfem::TMOP_Integrator::ParUpdateAfterMeshTopologyChange
void ParUpdateAfterMeshTopologyChange()
Definition: tmop.cpp:2478
mfem::TMOPComboIntegrator::GetDerefinementElementEnergy
virtual double GetDerefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
Definition: tmop.cpp:3756
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:2783
mfem::FiniteElementSpace::GetUpdateOperator
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
Definition: fespace.hpp:873
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:1240
mfem::DiscreteAdaptTC::GetTspecPert1H
const Vector & GetTspecPert1H()
Definition: tmop.hpp:1239
mfem::TMOP_Integrator::AssembleElementGradFD
void AssembleElementGradFD(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Definition: tmop.cpp:3360
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:1407
mfem::TMOP_Metric_007::ie
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:207
mfem::Mesh::ncmesh
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition: mesh.hpp:273
mfem::TMOP_Metric_aspratio3D::EvalW
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
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
t
RefCoord t[3]
Definition: ncmesh_tables.hpp:158
mfem::Coefficient::Eval
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
mfem::TMOP_Integrator::adapt_lim_gf
GridFunction * adapt_lim_gf
Definition: tmop.hpp:1339
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
mfem::TMOP_Integrator::surf_fit_gf
GridFunction * surf_fit_gf
Definition: tmop.hpp:1344
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:2633
mfem::DiscreteAdaptTC::tspec_pgf
ParGridFunction * tspec_pgf
Definition: tmop.hpp:1124
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::FiniteElementSpace::FEColl
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:577
mfem::DiscreteAdaptTC::GetDiscreteTargetSpec
void GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
Get one of the discrete fields from tspec.
Definition: tmop.cpp:1502
mfem::DenseMatrix::Mult
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:160
mfem::TMOP_Integrator::DSh
DenseMatrix DSh
Definition: tmop.hpp:1376
mfem::Mesh::GetNodes
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
mfem::TMOP_Integrator::surf_fit_marker
const Array< bool > * surf_fit_marker
Definition: tmop.hpp:1345
mfem::TMOP_Integrator::DS
DenseMatrix DS
Definition: tmop.hpp:1376
mfem::TMOP_Metric_303::ie
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:487
mfem::DiscreteAdaptTC::SetTspecDataForDerefinement
void SetTspecDataForDerefinement(FiniteElementSpace *fes)
Computes target specification data with respect to the coarse FE space.
Definition: tmop.cpp:1609
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::TMOP_Integrator::EnableSurfaceFitting
void EnableSurfaceFitting(const GridFunction &s0, const Array< bool > &smarker, Coefficient &coeff, AdaptivityEvaluator &ae)
Fitting of certain DOFs to the zero level set of a function.
Definition: tmop.cpp:2400
mfem::Mesh::DeleteGeometricFactors
void DeleteGeometricFactors()
Destroy all GeometricFactors stored by the Mesh.
Definition: mesh.cpp:881
s
RefCoord s[3]
Definition: ncmesh_tables.hpp:158
mfem::TMOP_Metric_211::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp: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:1528
mfem::DiscreteAdaptTC::SetParDiscreteTargetSkew
virtual void SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
Definition: tmop.cpp:1383
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::TMOPComboIntegrator::GetRefinementElementEnergy
virtual double GetRefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, const IntegrationRule &irule)
Definition: tmop.cpp:3743
mfem::TMOP_Integrator::exact_action
bool exact_action
Definition: tmop.hpp:1361
mfem::DiscreteAdaptTC::SetParDiscreteTargetOrientation
virtual void SetParDiscreteTargetOrientation(const ParGridFunction &tspec_)
Definition: tmop.cpp:1399
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:908
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::DiscreteAdaptTC::~DiscreteAdaptTC
virtual ~DiscreteAdaptTC()
Definition: tmop.cpp:2259
mfem::Operator
Abstract operator.
Definition: operator.hpp:24
mfem::IsoparametricTransformation::GetPointMat
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:404
mfem::TMOP_Metric_002::EvalP
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp: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:786
mfem::ParMesh
Class for parallel meshes.
Definition: pmesh.hpp:32
mfem::DiscreteAdaptTC::SetSerialDiscreteTargetSize
virtual void SetSerialDiscreteTargetSize(const GridFunction &tspec_)
Definition: tmop.cpp:1450
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:379
mfem::TMOP_Integrator::lim_func
TMOP_LimiterFunction * lim_func
Definition: tmop.hpp:1330
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::TMOP_Integrator::adapt_lim_gf0
const GridFunction * adapt_lim_gf0
Definition: tmop.hpp:1335
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:3711
mfem::TMOP_Integrator::dim
int dim
Definition: tmop.hpp:1405
mfem::TMOP_Integrator::Jpt
DenseMatrix Jpt
Definition: tmop.hpp:1376
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:438
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:3490
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:1376
mfem::TargetConstructor::GetComm
MPI_Comm GetComm() const
Definition: tmop.hpp:978
mfem::ElementTransformation::mesh
class Mesh * mesh
The Mesh object containing the element.
Definition: eltrans.hpp:84
mfem::TMOP_Metric_352::AssembleH
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
Definition: tmop.cpp:912
mfem::AdaptivityEvaluator::ComputeAtNewPosition
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field)=0
mfem::DiscreteAdaptTC::UpdateAfterMeshTopologyChange
void UpdateAfterMeshTopologyChange()
Update all discrete fields based on tspec and update for AMR.
Definition: tmop.cpp:1516
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:977
mfem::TMOP_Integrator::fd_call_flag
bool fd_call_flag
Definition: tmop.hpp:1358
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:858
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:113
mfem::DiscreteAdaptTC::coarse_tspec_fesv
FiniteElementSpace * coarse_tspec_fesv
Definition: tmop.hpp:1118
mfem::InvariantsEvaluator3D::SetJacobian
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:734
Generated on Mon Mar 21 2022 16:49:50 for MFEM by   doxygen 1.8.5