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