MFEM  v4.1.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 
29 {
30  ie.SetJacobian(Jpt.GetData());
31  P = ie.Get_dI1();
32 }
33 
35  const DenseMatrix &DS,
36  const double weight,
37  DenseMatrix &A) const
38 {
39  ie.SetJacobian(Jpt.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 
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 
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 double TMOP_Metric_002::EvalW(const DenseMatrix &Jpt) const
164 {
165  ie.SetJacobian(Jpt.GetData());
166  return 0.5 * ie.Get_I1b() - 1.0;
167 }
168 
170 {
171  ie.SetJacobian(Jpt.GetData());
172  P.Set(0.5, ie.Get_dI1b());
173 }
174 
176  const DenseMatrix &DS,
177  const double weight,
178  DenseMatrix &A) const
179 {
180  ie.SetJacobian(Jpt.GetData());
181  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
182  ie.Assemble_ddI1b(0.5*weight, A.GetData());
183 }
184 
185 double TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
186 {
187  // mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
188  ie.SetJacobian(Jpt.GetData());
189  return ie.Get_I1()*(1. + 1./ie.Get_I2()) - 4.0;
190 }
191 
193 {
194  // P = d(I1*(1 + 1/I2)) = (1 + 1/I2) dI1 - I1/I2^2 dI2
195  ie.SetJacobian(Jpt.GetData());
196  const double I2 = ie.Get_I2();
197  Add(1. + 1./I2, ie.Get_dI1(), -ie.Get_I1()/(I2*I2), ie.Get_dI2(), P);
198 }
199 
201  const DenseMatrix &DS,
202  const double weight,
203  DenseMatrix &A) const
204 {
205  // P = d(I1*(1 + 1/I2))
206  // = (1 + 1/I2) dI1 - I1/I2^2 dI2
207  //
208  // dP = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 -
209  // (dI2 x d(I1/I2^2)) - I1/I2^2 ddI2
210  // = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 +
211  // (-1/I2^2) (dI2 x [dI1 - 2 I1/I2 dI2]) - I1/I2^2 ddI2
212  // = (-1/I2^2) (dI1 x dI2 + dI2 x dI1) + (1 + 1/I2) ddI1 +
213  // (2 I1/I2^3) (dI2 x dI2) - I1/I2^2 ddI2
214  ie.SetJacobian(Jpt.GetData());
215  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
216  const double c1 = 1./ie.Get_I2();
217  const double c2 = weight*c1*c1;
218  const double c3 = ie.Get_I1()*c2;
219  ie.Assemble_ddI1(weight*(1. + c1), A.GetData());
220  ie.Assemble_ddI2(-c3, A.GetData());
221  ie.Assemble_TProd(-c2, ie.Get_dI1(), ie.Get_dI2(), A.GetData());
222  ie.Assemble_TProd(2*c1*c3, ie.Get_dI2(), A.GetData());
223 }
224 
225 double TMOP_Metric_009::EvalW(const DenseMatrix &Jpt) const
226 {
227  // mu_9 = det(J)*|J-J^{-t}|^2 = I1b * (I2b^2 + 1) - 4 * I2b
228  // = (I1 - 4)*I2b + I1b
229  ie.SetJacobian(Jpt.GetData());
230  return (ie.Get_I1() - 4.0)*ie.Get_I2b() + ie.Get_I1b();
231 }
232 
234 {
235  // mu_9 = (I1 - 4)*I2b + I1b
236  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
237  ie.SetJacobian(Jpt.GetData());
238  Add(ie.Get_I1() - 4.0, ie.Get_dI2b(), ie.Get_I2b(), ie.Get_dI1(), P);
239  P += ie.Get_dI1b();
240 }
241 
243  const DenseMatrix &DS,
244  const double weight,
245  DenseMatrix &A) const
246 {
247  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
248  // dP = dI2b x dI1 + (I1-4)*ddI2b + dI1 x dI2b + I2b*ddI1 + ddI1b
249  // = (dI1 x dI2b + dI2b x dI1) + (I1-4)*ddI2b + I2b*ddI1 + ddI1b
250  ie.SetJacobian(Jpt.GetData());
251  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
252  ie.Assemble_TProd(weight, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
253  ie.Assemble_ddI2b(weight*(ie.Get_I1()-4.0), A.GetData());
254  ie.Assemble_ddI1(weight*ie.Get_I2b(), A.GetData());
255  ie.Assemble_ddI1b(weight, A.GetData());
256 }
257 
258 double TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const
259 {
260  // mu_22 = (0.5*|J|^2 - det(J)) / (det(J) - tau0)
261  // = (0.5*I1 - I2b) / (I2b - tau0)
262  ie.SetJacobian(Jpt.GetData());
263  const double I2b = ie.Get_I2b();
264  return (0.5*ie.Get_I1() - I2b) / (I2b - tau0);
265 }
266 
268 {
269  // mu_22 = (0.5*I1 - I2b) / (I2b - tau0)
270  // P = 1/(I2b - tau0)*(0.5*dI1 - dI2b) - (0.5*I1 - I2b)/(I2b - tau0)^2*dI2b
271  // = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
272  ie.SetJacobian(Jpt.GetData());
273  const double c1 = 1.0/(ie.Get_I2b() - tau0);
274  Add(c1/2, ie.Get_dI1(), (tau0 - ie.Get_I1()/2)*c1*c1, ie.Get_dI2b(), P);
275 }
276 
278  const DenseMatrix &DS,
279  const double weight,
280  DenseMatrix &A) const
281 {
282  // P = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
283  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b) + 0.5/(I2b - tau0)*ddI1
284  // + (dI2b x dz) + z*ddI2b
285  //
286  // z = (tau0 - 0.5*I1)/(I2b - tau0)^2
287  // dz = -0.5/(I2b - tau0)^2*dI1 - 2*(tau0 - 0.5*I1)/(I2b - tau0)^3*dI2b
288  //
289  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b + dI2b x dI1)
290  // -2*z/(I2b - tau0)*(dI2b x dI2b)
291  // +0.5/(I2b - tau0)*ddI1 + z*ddI2b
292  ie.SetJacobian(Jpt.GetData());
293  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
294  const double c1 = 1.0/(ie.Get_I2b() - tau0);
295  const double c2 = weight*c1/2;
296  const double c3 = c1*c2;
297  const double c4 = (2*tau0 - ie.Get_I1())*c3; // weight*z
298  ie.Assemble_TProd(-c3, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
299  ie.Assemble_TProd(-2*c1*c4, ie.Get_dI2b(), A.GetData());
300  ie.Assemble_ddI1(c2, A.GetData());
301  ie.Assemble_ddI2b(c4, A.GetData());
302 }
303 
304 double TMOP_Metric_050::EvalW(const DenseMatrix &Jpt) const
305 {
306  // mu_50 = 0.5*|J^t J|^2/det(J)^2 - 1
307  // = 0.5*(l1^4 + l2^4)/(l1*l2)^2 - 1
308  // = 0.5*((l1/l2)^2 + (l2/l1)^2) - 1 = 0.5*(l1/l2 - l2/l1)^2
309  // = 0.5*(l1/l2 + l2/l1)^2 - 2 = 0.5*I1b^2 - 2
310  ie.SetJacobian(Jpt.GetData());
311  const double I1b = ie.Get_I1b();
312  return 0.5*I1b*I1b - 2.0;
313 }
314 
316 {
317  // mu_50 = 0.5*I1b^2 - 2
318  // P = I1b*dI1b
319  ie.SetJacobian(Jpt.GetData());
320  P.Set(ie.Get_I1b(), ie.Get_dI1b());
321 }
322 
324  const DenseMatrix &DS,
325  const double weight,
326  DenseMatrix &A) const
327 {
328  // P = I1b*dI1b
329  // dP = dI1b x dI1b + I1b*ddI1b
330  ie.SetJacobian(Jpt.GetData());
331  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
332  ie.Assemble_TProd(weight, ie.Get_dI1b(), A.GetData());
333  ie.Assemble_ddI1b(weight*ie.Get_I1b(), A.GetData());
334 }
335 
336 double TMOP_Metric_055::EvalW(const DenseMatrix &Jpt) const
337 {
338  // mu_55 = (det(J) - 1)^2 = (I2b - 1)^2
339  ie.SetJacobian(Jpt.GetData());
340  const double c1 = ie.Get_I2b() - 1.0;
341  return c1*c1;
342 }
343 
345 {
346  // mu_55 = (I2b - 1)^2
347  // P = 2*(I2b - 1)*dI2b
348  ie.SetJacobian(Jpt.GetData());
349  P.Set(2*(ie.Get_I2b() - 1.0), ie.Get_dI2b());
350 }
351 
353  const DenseMatrix &DS,
354  const double weight,
355  DenseMatrix &A) const
356 {
357  // P = 2*(I2b - 1)*dI2b
358  // dP = 2*(dI2b x dI2b) + 2*(I2b - 1)*ddI2b
359  ie.SetJacobian(Jpt.GetData());
360  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
361  ie.Assemble_TProd(2*weight, ie.Get_dI2b(), A.GetData());
362  ie.Assemble_ddI2b(2*weight*(ie.Get_I2b() - 1.0), A.GetData());
363 }
364 
365 double TMOP_Metric_056::EvalW(const DenseMatrix &Jpt) const
366 {
367  // mu_56 = 0.5*(I2b + 1/I2b) - 1
368  ie.SetJacobian(Jpt.GetData());
369  const double I2b = ie.Get_I2b();
370  return 0.5*(I2b + 1.0/I2b) - 1.0;
371 }
372 
374 {
375  // mu_56 = 0.5*(I2b + 1/I2b) - 1
376  // P = 0.5*(1 - 1/I2b^2)*dI2b
377  ie.SetJacobian(Jpt.GetData());
378  P.Set(0.5 - 0.5/ie.Get_I2(), ie.Get_dI2b());
379 }
380 
382  const DenseMatrix &DS,
383  const double weight,
384  DenseMatrix &A) const
385 {
386  // P = 0.5*(1 - 1/I2b^2)*dI2b
387  // dP = (1/I2b^3)*(dI2b x dI2b) + (0.5 - 0.5/I2)*ddI2b
388  ie.SetJacobian(Jpt.GetData());
389  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
390  ie.Assemble_TProd(weight/(ie.Get_I2()*ie.Get_I2b()),
391  ie.Get_dI2b(), A.GetData());
392  ie.Assemble_ddI2b(weight*(0.5 - 0.5/ie.Get_I2()), A.GetData());
393 }
394 
395 double TMOP_Metric_058::EvalW(const DenseMatrix &Jpt) const
396 {
397  // mu_58 = I1b*(I1b - 2)
398  ie.SetJacobian(Jpt.GetData());
399  const double I1b = ie.Get_I1b();
400  return I1b*(I1b - 1.0);
401 }
402 
404 {
405  // mu_58 = I1b*(I1b - 2)
406  // P = (2*I1b - 2)*dI1b
407  ie.SetJacobian(Jpt.GetData());
408  P.Set(2*ie.Get_I1b() - 2.0, ie.Get_dI1b());
409 }
410 
412  const DenseMatrix &DS,
413  const double weight,
414  DenseMatrix &A) const
415 {
416  // P = (2*I1b - 2)*dI1b
417  // dP = 2*(dI1b x dI1b) + (2*I1b - 2)*ddI1b
418  ie.SetJacobian(Jpt.GetData());
419  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
420  ie.Assemble_TProd(2*weight, ie.Get_dI1b(), A.GetData());
421  ie.Assemble_ddI1b(weight*(2*ie.Get_I1b() - 2.0), A.GetData());
422 }
423 
424 double TMOP_Metric_077::EvalW(const DenseMatrix &Jpt) const
425 {
426  ie.SetJacobian(Jpt.GetData());
427  const double I2 = ie.Get_I2b();
428  return 0.5*(I2*I2 + 1./(I2*I2) - 2.);
429 }
430 
432 {
433  // Using I2b^2 = I2.
434  // dmu77_dJ = 1/2 (1 - 1/I2^2) dI2_dJ.
435  ie.SetJacobian(Jpt.GetData());
436  const double I2 = ie.Get_I2();
437  P.Set(0.5 * (1.0 - 1.0 / (I2 * I2)), ie.Get_dI2());
438 }
439 
441  const DenseMatrix &DS,
442  const double weight,
443  DenseMatrix &A) const
444 {
445  ie.SetJacobian(Jpt.GetData());
446  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
447  const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
448  ie.Assemble_ddI2(weight*0.5*(1.0 - I2inv_sq), A.GetData());
449  ie.Assemble_TProd(weight * I2inv_sq / I2, ie.Get_dI2(), A.GetData());
450 }
451 
452 double TMOP_Metric_211::EvalW(const DenseMatrix &Jpt) const
453 {
454  // mu_211 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
455  // = (I2b - 1)^2 - I2b + sqrt(I2b^2 + eps)
456  ie.SetJacobian(Jpt.GetData());
457  const double I2b = ie.Get_I2b();
458  return (I2b - 1.0)*(I2b - 1.0) - I2b + std::sqrt(I2b*I2b + eps);
459 }
460 
462 {
463  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
464 }
465 
467  const DenseMatrix &DS,
468  const double weight,
469  DenseMatrix &A) const
470 {
471  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
472 }
473 
474 double TMOP_Metric_252::EvalW(const DenseMatrix &Jpt) const
475 {
476  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0).
477  ie.SetJacobian(Jpt.GetData());
478  const double I2b = ie.Get_I2b();
479  return 0.5*(I2b - 1.0)*(I2b - 1.0)/(I2b - tau0);
480 }
481 
483 {
484  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
485  // P = (c - 0.5*c*c ) * dI2b
486  //
487  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
488  ie.SetJacobian(Jpt.GetData());
489  const double I2b = ie.Get_I2b();
490  const double c = (I2b - 1.0)/(I2b - tau0);
491  P.Set(c - 0.5*c*c, ie.Get_dI2b());
492 }
493 
495  const DenseMatrix &DS,
496  const double weight,
497  DenseMatrix &A) const
498 {
499  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
500  //
501  // P = (c - 0.5*c*c ) * dI2b
502  // dP = (1 - c)^2/(I2b - tau0)*(dI2b x dI2b) + (c - 0.5*c*c)*ddI2b
503  ie.SetJacobian(Jpt.GetData());
504  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
505  const double I2b = ie.Get_I2b();
506  const double c0 = 1.0/(I2b - tau0);
507  const double c = c0*(I2b - 1.0);
508  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI2b(), A.GetData());
509  ie.Assemble_ddI2b(weight*(c - 0.5*c*c), A.GetData());
510 }
511 
512 double TMOP_Metric_301::EvalW(const DenseMatrix &Jpt) const
513 {
514  ie.SetJacobian(Jpt.GetData());
515  return std::sqrt(ie.Get_I1b()*ie.Get_I2b())/3. - 1.;
516 }
517 
519 {
520  // W = (1/3)*sqrt(I1b*I2b) - 1
521  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
522  ie.SetJacobian(Jpt.GetData());
523  const double a = 1./(6.*std::sqrt(ie.Get_I1b()*ie.Get_I2b()));
524  Add(a*ie.Get_I2b(), ie.Get_dI1b(), a*ie.Get_I1b(), ie.Get_dI2b(), P);
525 }
526 
528  const DenseMatrix &DS,
529  const double weight,
530  DenseMatrix &A) const
531 {
532  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
533  // dW = (1/6)*[z2*dI1b + z1*dI2b], z1 = sqrt(I1b/I2b), z2 = sqrt(I2b/I1b)
534  // ddW = (1/6)*[dI1b x dz2 + z2*ddI1b + dI2b x dz1 + z1*ddI2b]
535  //
536  // dz1 = (1/2)*sqrt(I2b/I1b) [ (1/I2b)*dI1b + (I1b/(I2b*I2b))*dI2b ]
537  // = (1/2)/sqrt(I1b*I2b) [ dI1b + (I1b/I2b)*dI2b ]
538  // dz2 = (1/2)/sqrt(I1b*I2b) [ (I2b/I1b)*dI1b + dI2b ]
539  //
540  // dI1b x dz2 + dI2b x dz1 =
541  // (1/2)/sqrt(I1b*I2b) dI1b x [ (I2b/I1b)*dI1b + dI2b ] +
542  // (1/2)/sqrt(I1b*I2b) dI2b x [ dI1b + (I1b/I2b)*dI2b ] =
543  // (1/2)/sqrt(I1b*I2b) [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] x
544  // [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] =
545  // (1/2)/sqrt(I1b*I2b) [ 6*dW x 6*dW ] =
546  // (1/2)*(I1b*I2b)^{-3/2} (I2b*dI1b + I1b*dI2b) x (I2b*dI1b + I1b*dI2b)
547  //
548  // z1 = I1b/sqrt(I1b*I2b), z2 = I2b/sqrt(I1b*I2b)
549 
550  ie.SetJacobian(Jpt.GetData());
551  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
552  double d_I1b_I2b_data[9];
553  DenseMatrix d_I1b_I2b(d_I1b_I2b_data, 3, 3);
554  Add(ie.Get_I2b(), ie.Get_dI1b(), ie.Get_I1b(), ie.Get_dI2b(), d_I1b_I2b);
555  const double I1b_I2b = ie.Get_I1b()*ie.Get_I2b();
556  const double a = weight/(6*std::sqrt(I1b_I2b));
557  ie.Assemble_ddI1b(a*ie.Get_I2b(), A.GetData());
558  ie.Assemble_ddI2b(a*ie.Get_I1b(), A.GetData());
559  ie.Assemble_TProd(a/(2*I1b_I2b), d_I1b_I2b_data, A.GetData());
560 }
561 
562 double TMOP_Metric_302::EvalW(const DenseMatrix &Jpt) const
563 {
564  // mu_2 = |J|^2 |J^{-1}|^2 / 9 - 1
565  // = (l1^2 + l2^2 + l3^3)*(l1^{-2} + l2^{-2} + l3^{-2}) / 9 - 1
566  // = I1*(l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/l1^2/l2^2/l3^2/9 - 1
567  // = I1*I2/det(J)^2/9 - 1 = I1b*I2b/9-1
568  ie.SetJacobian(Jpt.GetData());
569  return ie.Get_I1b()*ie.Get_I2b()/9. - 1.;
570 }
571 
573 {
574  // mu_2 = I1b*I2b/9-1
575  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
576  ie.SetJacobian(Jpt.GetData());
577  Add(ie.Get_I1b()/9, ie.Get_dI2b(), ie.Get_I2b()/9, ie.Get_dI1b(), P);
578 }
579 
581  const DenseMatrix &DS,
582  const double weight,
583  DenseMatrix &A) const
584 {
585  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
586  // dP = (dI2b x dI1b)/9 + (I1b/9)*ddI2b + (dI1b x dI2b)/9 + (I2b/9)*ddI1b
587  // = (dI2b x dI1b + dI1b x dI2b)/9 + (I1b/9)*ddI2b + (I2b/9)*ddI1b
588  ie.SetJacobian(Jpt.GetData());
589  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
590  const double c1 = weight/9;
591  ie.Assemble_TProd(c1, ie.Get_dI1b(), ie.Get_dI2b(), A.GetData());
592  ie.Assemble_ddI2b(c1*ie.Get_I1b(), A.GetData());
593  ie.Assemble_ddI1b(c1*ie.Get_I2b(), A.GetData());
594 }
595 
596 double TMOP_Metric_303::EvalW(const DenseMatrix &Jpt) const
597 {
598  ie.SetJacobian(Jpt.GetData());
599  return ie.Get_I1b()/3.0 - 1.0;
600 }
601 
603 {
604  ie.SetJacobian(Jpt.GetData());
605  P.Set(1./3., ie.Get_dI1b());
606 }
607 
609  const DenseMatrix &DS,
610  const double weight,
611  DenseMatrix &A) const
612 {
613  ie.SetJacobian(Jpt.GetData());
614  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
615  ie.Assemble_ddI1b(weight/3., A.GetData());
616 }
617 
618 double TMOP_Metric_315::EvalW(const DenseMatrix &Jpt) const
619 {
620  // mu_315 = mu_15_3D = (det(J) - 1)^2
621  ie.SetJacobian(Jpt.GetData());
622  const double c1 = ie.Get_I3b() - 1.0;
623  return c1*c1;
624 }
625 
627 {
628  // mu_315 = (I3b - 1)^2
629  // P = 2*(I3b - 1)*dI3b
630  ie.SetJacobian(Jpt.GetData());
631  P.Set(2*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
632 }
633 
635  const DenseMatrix &DS,
636  const double weight,
637  DenseMatrix &A) const
638 {
639  // P = 2*(I3b - 1)*dI3b
640  // dP = 2*(dI3b x dI3b) + 2*(I3b - 1)*ddI3b
641  ie.SetJacobian(Jpt.GetData());
642  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
643  ie.Assemble_TProd(2*weight, ie.Get_dI3b(), A.GetData());
644  ie.Assemble_ddI3b(2*weight*(ie.Get_I3b() - 1.0), A.GetData());
645 }
646 
647 double TMOP_Metric_316::EvalW(const DenseMatrix &Jpt) const
648 {
649  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
650  ie.SetJacobian(Jpt.GetData());
651  const double I3b = ie.Get_I3b();
652  return 0.5*(I3b + 1.0/I3b) - 1.0;
653 }
654 
656 {
657  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
658  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
659  ie.SetJacobian(Jpt.GetData());
660  P.Set(0.5 - 0.5/ie.Get_I3(), ie.Get_dI3b());
661 }
662 
664  const DenseMatrix &DS,
665  const double weight,
666  DenseMatrix &A) const
667 {
668  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
669  // dP = (1/I3b^3)*(dI3b x dI3b) + (0.5 - 0.5/I3)*ddI3b
670  ie.SetJacobian(Jpt.GetData());
671  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
672  ie.Assemble_TProd(weight/(ie.Get_I3()*ie.Get_I3b()),
673  ie.Get_dI3b(), A.GetData());
674  ie.Assemble_ddI3b(weight*(0.5 - 0.5/ie.Get_I3()), A.GetData());
675 }
676 
677 double TMOP_Metric_321::EvalW(const DenseMatrix &Jpt) const
678 {
679  // mu_321 = mu_21_3D = |J - J^{-t}|^2
680  // = |J|^2 + |J^{-1}|^2 - 6
681  // = |J|^2 + (l1^{-2} + l2^{-2} + l3^{-2}) - 6
682  // = |J|^2 + (l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/det(J)^2 - 6
683  // = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
684  ie.SetJacobian(Jpt.GetData());
685  return ie.Get_I1() + ie.Get_I2()/ie.Get_I3() - 6.0;
686 }
687 
689 {
690  // mu_321 = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
691  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
692  ie.SetJacobian(Jpt.GetData());
693  const double I3 = ie.Get_I3();
694  Add(1.0/I3, ie.Get_dI2(),
695  -2*ie.Get_I2()/(I3*ie.Get_I3b()), ie.Get_dI3b(), P);
696  P += ie.Get_dI1();
697 }
698 
700  const DenseMatrix &DS,
701  const double weight,
702  DenseMatrix &A) const
703 {
704  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
705  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b) + (1/I3)*ddI2 + (dI3b x dz) + z*ddI3b
706  //
707  // z = -2*I2/I3b^3
708  // dz = (-2/I3b^3)*dI2 + (2*I2)*(3/I3b^4)*dI3b
709  //
710  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b + dI3b x dI2) + (1/I3)*ddI2
711  // + (6*I2/I3b^4)*(dI3b x dI3b) + (-2*I2/I3b^3)*ddI3b
712  ie.SetJacobian(Jpt.GetData());
713  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
714  const double c0 = 1.0/ie.Get_I3b();
715  const double c1 = weight*c0*c0;
716  const double c2 = -2*c0*c1;
717  const double c3 = c2*ie.Get_I2();
718  ie.Assemble_ddI1(weight, A.GetData());
719  ie.Assemble_ddI2(c1, A.GetData());
720  ie.Assemble_ddI3b(c3, A.GetData());
721  ie.Assemble_TProd(c2, ie.Get_dI2(), ie.Get_dI3b(), A.GetData());
722  ie.Assemble_TProd(-3*c0*c3, ie.Get_dI3b(), A.GetData());
723 }
724 
725 double TMOP_Metric_352::EvalW(const DenseMatrix &Jpt) const
726 {
727  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
728  ie.SetJacobian(Jpt.GetData());
729  const double I3b = ie.Get_I3b();
730  return 0.5*(I3b - 1.0)*(I3b - 1.0)/(I3b - tau0);
731 }
732 
734 {
735  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
736  // P = (I3b - 1)/(I3b - tau0)*dI3b + 0.5*(I3b - 1)^2*(-1/(I3b - tau0)^2)*dI3b
737  // = [ (I3b - 1)/(I3b - tau0) - 0.5*(I3b - 1)^2/(I3b - tau0)^2 ] * dI3b
738  // = (c - 0.5*c*c) * dI3b
739  ie.SetJacobian(Jpt.GetData());
740  const double I3b = ie.Get_I3b();
741  const double c = (I3b - 1.0)/(I3b - tau0);
742  P.Set(c - 0.5*c*c, ie.Get_dI3b());
743 }
744 
746  const DenseMatrix &DS,
747  const double weight,
748  DenseMatrix &A) const
749 {
750  // c = (I3b - 1)/(I3b - tau0)
751  //
752  // P = (c - 0.5*c*c) * dI3b
753  // dP = (1 - c)*(dI3b x dc) + (c - 0.5*c*c)*ddI3b
754  //
755  // dc = 1/(I3b - tau0)*dI3b - (I3b - 1)/(I3b - tau)^2*dI3b =
756  // = (1 - c)/(I3b - tau0)*dI3b
757  //
758  // dP = (1 - c)^2/(I3b - tau0)*(dI3b x dI3b) + (c - 0.5*c*c)*ddI3b
759  ie.SetJacobian(Jpt.GetData());
760  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
761  const double I3b = ie.Get_I3b();
762  const double c0 = 1.0/(I3b - tau0);
763  const double c = c0*(I3b - 1.0);
764  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI3b(), A.GetData());
765  ie.Assemble_ddI3b(weight*(c - 0.5*c*c), A.GetData());
766 }
767 
768 
770 {
771  MFEM_VERIFY(nodes, "Nodes are not given!");
772  MFEM_ASSERT(avg_volume == 0.0, "The average volume is already computed!");
773 
774  Mesh *mesh = nodes->FESpace()->GetMesh();
775  const int NE = mesh->GetNE();
777  double volume = 0.0;
778 
779  for (int i = 0; i < NE; i++)
780  {
781  mesh->GetElementTransformation(i, *nodes, &Tr);
782  const IntegrationRule &ir =
783  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr.OrderJ());
784  for (int j = 0; j < ir.GetNPoints(); j++)
785  {
786  const IntegrationPoint &ip = ir.IntPoint(j);
787  Tr.SetIntPoint(&ip);
788  volume += ip.weight * Tr.Weight();
789  }
790  }
791 
792  NCMesh *ncmesh = mesh->ncmesh;
793  if (Parallel() == false)
794  {
795  avg_volume = (ncmesh == NULL) ?
796  volume / NE : volume / ncmesh->GetNumRootElements();
797 
798  }
799 #ifdef MFEM_USE_MPI
800  else
801  {
802  double area_NE[4];
803  area_NE[0] = volume; area_NE[1] = NE;
804  MPI_Allreduce(area_NE, area_NE + 2, 2, MPI_DOUBLE, MPI_SUM, comm);
805  avg_volume = (ncmesh == NULL) ?
806  area_NE[2] / area_NE[3] : area_NE[2] / ncmesh->GetNumRootElements();
807  }
808 #endif
809 }
810 
811 // virtual method
813  const IntegrationRule &ir,
814  const Vector &elfun,
815  DenseTensor &Jtr) const
816 {
817  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
818 
820  nodes->FESpace()->GetFE(e_id) : NULL;
821  const DenseMatrix &Wideal =
823  MFEM_ASSERT(Wideal.Height() == Jtr.SizeI(), "");
824  MFEM_ASSERT(Wideal.Width() == Jtr.SizeJ(), "");
825 
826  switch (target_type)
827  {
829  {
830  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = Wideal; }
831  break;
832  }
834  {
835  if (avg_volume == 0.0) { ComputeAvgVolume(); }
836  DenseMatrix W(Wideal.Height());
837 
838  NCMesh *ncmesh = nodes->FESpace()->GetMesh()->ncmesh;
839  double el_volume = avg_volume;
840  if (ncmesh)
841  {
842  el_volume = avg_volume / ncmesh->GetElementSizeReduction(e_id);
843  }
844 
845  W.Set(std::pow(volume_scale * el_volume / Wideal.Det(),
846  1./W.Height()), Wideal);
847  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = W; }
848  break;
849  }
852  {
853  const int dim = nfe->GetDim(), dof = nfe->GetDof();
854  MFEM_ASSERT(dim == nodes->FESpace()->GetVDim(), "");
855  DenseMatrix dshape(dof, dim), pos(dof, dim);
856  Array<int> xdofs(dof * dim);
857  Vector posV(pos.Data(), dof * dim);
858  double detW;
859 
860  // always initialize detW to suppress a warning:
861  detW = (target_type == IDEAL_SHAPE_GIVEN_SIZE) ? Wideal.Det() : 0.0;
862  nodes->FESpace()->GetElementVDofs(e_id, xdofs);
863  nodes->GetSubVector(xdofs, posV);
864  for (int i = 0; i < ir.GetNPoints(); i++)
865  {
866  nfe->CalcDShape(ir.IntPoint(i), dshape);
867  MultAtB(pos, dshape, Jtr(i));
869  {
870  const double det = Jtr(i).Det();
871  MFEM_VERIFY(det > 0.0, "The given mesh is inverted!");
872  Jtr(i).Set(std::pow(det / detW, 1./dim), Wideal);
873  }
874  }
875  break;
876  }
877  default:
878  MFEM_ABORT("invalid target type!");
879  }
880 }
881 
883  VectorCoefficient *vspec,
884  MatrixCoefficient *mspec)
885 {
886  scalar_tspec = sspec;
887  vector_tspec = vspec;
888  matrix_tspec = mspec;
889 }
890 
892  const IntegrationRule &ir,
893  const Vector &elfun,
894  DenseTensor &Jtr) const
895 {
896  DenseMatrix point_mat;
897  point_mat.UseExternalData(elfun.GetData(), fe.GetDof(), fe.GetDim());
898 
899  switch (target_type)
900  {
901  case GIVEN_FULL:
902  {
903  MFEM_VERIFY(matrix_tspec != NULL,
904  "Target type GIVEN_FULL requires a MatrixCoefficient.");
905 
907  Tpr.SetFE(&fe);
908  Tpr.ElementNo = e_id;
909  Tpr.GetPointMat().Transpose(point_mat);
910 
911  for (int i = 0; i < ir.GetNPoints(); i++)
912  {
913  const IntegrationPoint &ip = ir.IntPoint(i);
914  Tpr.SetIntPoint(&ip);
915  matrix_tspec->Eval(Jtr(i), Tpr, ip);
916  }
917  break;
918  }
919  default:
920  MFEM_ABORT("Incompatible target type for analytic adaptation!");
921  }
922 }
923 
924 #ifdef MFEM_USE_MPI
926 {
927  target_spec.SetSize(tspec.Size());
928  target_spec = tspec;
929  tspec_fes = tspec.FESpace();
930 
931  // Default evaluator is based on CG advection.
932  if (adapt_eval == NULL) { adapt_eval = new AdvectorCG; }
933 
935  *tspec.FESpace()->FEColl(),
936  tspec.FESpace()->GetVDim());
937 
939  (*tspec.FESpace()->GetMesh()->GetNodes(), target_spec);
940 }
941 #endif
942 
944 {
945  target_spec.SetSize(tspec.Size());
946  target_spec = tspec;
947  tspec_fes = tspec.FESpace();
948 
949  // Default evaluator is based on CG advection.
950  if (adapt_eval == NULL) { adapt_eval = new AdvectorCG; }
951 
953  *tspec.FESpace()->FEColl(),
954  tspec.FESpace()->GetVDim());
955 
957  (*tspec.FESpace()->GetMesh()->GetNodes(), target_spec);
958 }
959 
961 {
962  MFEM_VERIFY(target_spec.Size() > 0, "Target specification is not set!");
963 
965 }
966 
968  const IntegrationRule &ir,
969  const Vector &elfun,
970  DenseTensor &Jtr) const
971 {
972  MFEM_VERIFY(tspec_fes, "A call to SetDiscreteTargerSpec() is needed.");
973 
974  switch (target_type)
975  {
977  {
978  const DenseMatrix &Wideal =
980  const int dim = Wideal.Height(),
981  ntspec_dofs = tspec_fes->GetFE(0)->GetDof();
982 
983  Vector shape(ntspec_dofs), tspec_vals(ntspec_dofs);
984  Array<int> dofs;
985  tspec_fes->GetElementDofs(e_id, dofs);
986  target_spec.GetSubVector(dofs, tspec_vals);
987 
988  const double min_size = tspec_vals.Min();
989  MFEM_ASSERT(min_size > 0.0,
990  "Non-positive size propagated in the target definition.");
991 
992  for (int i = 0; i < ir.GetNPoints(); i++)
993  {
994  const IntegrationPoint &ip = ir.IntPoint(i);
995  tspec_fes->GetFE(e_id)->CalcShape(ip, shape);
996  const double size = std::max(shape * tspec_vals, min_size);
997  Jtr(i).Set(std::pow(size / Wideal.Det(), 1.0/dim), Wideal);
998  }
999  break;
1000  }
1001  default:
1002  MFEM_ABORT("Incompatible target type for analytic adaptation!");
1003  }
1004 }
1005 
1007  const FiniteElementCollection &fec,
1008  int num_comp)
1009 {
1010  delete fes;
1011  delete mesh;
1012  mesh = new Mesh(m, true);
1013  fes = new FiniteElementSpace(mesh, &fec, num_comp);
1014 }
1015 
1016 #ifdef MFEM_USE_MPI
1018  const FiniteElementCollection &fec,
1019  int num_comp)
1020 {
1021  delete pfes;
1022  delete pmesh;
1023  pmesh = new ParMesh(m, true);
1024  pfes = new ParFiniteElementSpace(pmesh, &fec, num_comp);
1025 }
1026 #endif
1027 
1029 {
1030  delete fes;
1031  delete mesh;
1032 #ifdef MFEM_USE_MPI
1033  delete pfes;
1034  delete pmesh;
1035 #endif
1036 }
1037 
1039  const GridFunction &dist, Coefficient &w0,
1040  TMOP_LimiterFunction *lfunc)
1041 {
1042  nodes0 = &n0;
1043  coeff0 = &w0;
1044  lim_dist = &dist;
1045 
1046  delete lim_func;
1047  if (lfunc)
1048  {
1049  lim_func = lfunc;
1050  }
1051  else
1052  {
1054  }
1055 }
1057  TMOP_LimiterFunction *lfunc)
1058 {
1059  nodes0 = &n0;
1060  coeff0 = &w0;
1061  lim_dist = NULL;
1062 
1063  delete lim_func;
1064  if (lfunc)
1065  {
1066  lim_func = lfunc;
1067  }
1068  else
1069  {
1071  }
1072 }
1073 
1076  const Vector &elfun)
1077 {
1078  int dof = el.GetDof(), dim = el.GetDim();
1079  double energy;
1080 
1081  DSh.SetSize(dof, dim);
1082  Jrt.SetSize(dim);
1083  Jpr.SetSize(dim);
1084  Jpt.SetSize(dim);
1085  PMatI.UseExternalData(elfun.GetData(), dof, dim);
1086 
1087  const IntegrationRule *ir = IntRule;
1088  if (!ir)
1089  {
1090  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
1091  }
1092 
1093  energy = 0.0;
1094  DenseTensor Jtr(dim, dim, ir->GetNPoints());
1095  targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
1096 
1097  // Limited case.
1098  Vector shape, p, p0, d_vals;
1099  DenseMatrix pos0;
1100  if (coeff0)
1101  {
1102  shape.SetSize(dof);
1103  p.SetSize(dim);
1104  p0.SetSize(dim);
1105  pos0.SetSize(dof, dim);
1106  Vector pos0V(pos0.Data(), dof * dim);
1107  Array<int> pos_dofs;
1108  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
1109  nodes0->GetSubVector(pos_dofs, pos0V);
1110  if (lim_dist)
1111  {
1112  lim_dist->GetValues(T.ElementNo, *ir, d_vals);
1113  }
1114  else
1115  {
1116  d_vals.SetSize(ir->GetNPoints()); d_vals = 1.0;
1117  }
1118  }
1119 
1120  // Define ref->physical transformation, when a Coefficient is specified.
1121  IsoparametricTransformation *Tpr = NULL;
1122  if (coeff1 || coeff0)
1123  {
1124  Tpr = new IsoparametricTransformation;
1125  Tpr->SetFE(&el);
1126  Tpr->ElementNo = T.ElementNo;
1127  Tpr->Attribute = T.Attribute;
1128  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
1129  }
1130  // TODO: computing the coefficients 'coeff1' and 'coeff0' in physical
1131  // coordinates means that, generally, the gradient and Hessian of the
1132  // TMOP_Integrator will depend on the derivatives of the coefficients.
1133  //
1134  // In some cases the coefficients are independent of any movement of
1135  // the physical coordinates (i.e. changes in 'elfun'), e.g. when the
1136  // coefficient is a ConstantCoefficient or a GridFunctionCoefficient.
1137 
1138  for (int i = 0; i < ir->GetNPoints(); i++)
1139  {
1140  const IntegrationPoint &ip = ir->IntPoint(i);
1141  const DenseMatrix &Jtr_i = Jtr(i);
1142  metric->SetTargetJacobian(Jtr_i);
1143  CalcInverse(Jtr_i, Jrt);
1144  const double weight = ip.weight * Jtr_i.Det();
1145 
1146  el.CalcDShape(ip, DSh);
1147  MultAtB(PMatI, DSh, Jpr);
1148  Mult(Jpr, Jrt, Jpt);
1149 
1150  double val = metric_normal * metric->EvalW(Jpt);
1151  if (coeff1) { val *= coeff1->Eval(*Tpr, ip); }
1152 
1153  if (coeff0)
1154  {
1155  el.CalcShape(ip, shape);
1156  PMatI.MultTranspose(shape, p);
1157  pos0.MultTranspose(shape, p0);
1158  val += lim_normal *
1159  lim_func->Eval(p, p0, d_vals(i)) * coeff0->Eval(*Tpr, ip);
1160  }
1161  energy += weight * val;
1162  }
1163  delete Tpr;
1164 
1165  return energy;
1166 }
1167 
1170  const Vector &elfun, Vector &elvect)
1171 {
1172  int dof = el.GetDof(), dim = el.GetDim();
1173 
1174  DSh.SetSize(dof, dim);
1175  DS.SetSize(dof, dim);
1176  Jrt.SetSize(dim);
1177  Jpt.SetSize(dim);
1178  P.SetSize(dim);
1179  PMatI.UseExternalData(elfun.GetData(), dof, dim);
1180  elvect.SetSize(dof*dim);
1181  PMatO.UseExternalData(elvect.GetData(), dof, dim);
1182 
1183  const IntegrationRule *ir = IntRule;
1184  if (!ir)
1185  {
1186  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
1187  }
1188 
1189  elvect = 0.0;
1190  DenseTensor Jtr(dim, dim, ir->GetNPoints());
1191  targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
1192 
1193  // Limited case.
1194  DenseMatrix pos0;
1195  Vector shape, p, p0, d_vals, grad;
1196  if (coeff0)
1197  {
1198  shape.SetSize(dof);
1199  p.SetSize(dim);
1200  p0.SetSize(dim);
1201  pos0.SetSize(dof, dim);
1202  Vector pos0V(pos0.Data(), dof * dim);
1203  Array<int> pos_dofs;
1204  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
1205  nodes0->GetSubVector(pos_dofs, pos0V);
1206  if (lim_dist)
1207  {
1208  lim_dist->GetValues(T.ElementNo, *ir, d_vals);
1209  }
1210  else
1211  {
1212  d_vals.SetSize(ir->GetNPoints()); d_vals = 1.0;
1213  }
1214  }
1215 
1216  // Define ref->physical transformation, when a Coefficient is specified.
1217  IsoparametricTransformation *Tpr = NULL;
1218  if (coeff1 || coeff0)
1219  {
1220  Tpr = new IsoparametricTransformation;
1221  Tpr->SetFE(&el);
1222  Tpr->ElementNo = T.ElementNo;
1223  Tpr->Attribute = T.Attribute;
1224  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
1225  }
1226 
1227  for (int i = 0; i < ir->GetNPoints(); i++)
1228  {
1229  const IntegrationPoint &ip = ir->IntPoint(i);
1230  const DenseMatrix &Jtr_i = Jtr(i);
1231  metric->SetTargetJacobian(Jtr_i);
1232  CalcInverse(Jtr_i, Jrt);
1233  const double weight = ip.weight * Jtr_i.Det();
1234  double weight_m = weight * metric_normal;
1235 
1236  el.CalcDShape(ip, DSh);
1237  Mult(DSh, Jrt, DS);
1238  MultAtB(PMatI, DS, Jpt);
1239 
1240  metric->EvalP(Jpt, P);
1241 
1242  if (coeff1) { weight_m *= coeff1->Eval(*Tpr, ip); }
1243 
1244  P *= weight_m;
1245  AddMultABt(DS, P, PMatO);
1246 
1247  // TODO: derivatives of adaptivity-based targets.
1248 
1249  if (coeff0)
1250  {
1251  el.CalcShape(ip, shape);
1252  PMatI.MultTranspose(shape, p);
1253  pos0.MultTranspose(shape, p0);
1254  lim_func->Eval_d1(p, p0, d_vals(i), grad);
1255  grad *= weight * lim_normal * coeff0->Eval(*Tpr, ip);
1256  AddMultVWt(shape, grad, PMatO);
1257  }
1258  }
1259  delete Tpr;
1260 }
1261 
1264  const Vector &elfun,
1265  DenseMatrix &elmat)
1266 {
1267  int dof = el.GetDof(), dim = el.GetDim();
1268 
1269  DSh.SetSize(dof, dim);
1270  DS.SetSize(dof, dim);
1271  Jrt.SetSize(dim);
1272  Jpt.SetSize(dim);
1273  PMatI.UseExternalData(elfun.GetData(), dof, dim);
1274  elmat.SetSize(dof*dim);
1275 
1276  const IntegrationRule *ir = IntRule;
1277  if (!ir)
1278  {
1279  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
1280  }
1281 
1282  elmat = 0.0;
1283  DenseTensor Jtr(dim, dim, ir->GetNPoints());
1284  targetC->ComputeElementTargets(T.ElementNo, el, *ir, elfun, Jtr);
1285 
1286  // Limited case.
1287  DenseMatrix pos0, grad_grad;
1288  Vector shape, p, p0, d_vals;
1289  if (coeff0)
1290  {
1291  shape.SetSize(dof);
1292  p.SetSize(dim);
1293  p0.SetSize(dim);
1294  pos0.SetSize(dof, dim);
1295  Vector pos0V(pos0.Data(), dof * dim);
1296  Array<int> pos_dofs;
1297  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
1298  nodes0->GetSubVector(pos_dofs, pos0V);
1299  if (lim_dist)
1300  {
1301  lim_dist->GetValues(T.ElementNo, *ir, d_vals);
1302  }
1303  else
1304  {
1305  d_vals.SetSize(ir->GetNPoints()); d_vals = 1.0;
1306  }
1307  }
1308 
1309  // Define ref->physical transformation, when a Coefficient is specified.
1310  IsoparametricTransformation *Tpr = NULL;
1311  if (coeff1 || coeff0)
1312  {
1313  Tpr = new IsoparametricTransformation;
1314  Tpr->SetFE(&el);
1315  Tpr->ElementNo = T.ElementNo;
1316  Tpr->Attribute = T.Attribute;
1317  Tpr->GetPointMat().Transpose(PMatI);
1318  }
1319 
1320  for (int i = 0; i < ir->GetNPoints(); i++)
1321  {
1322  const IntegrationPoint &ip = ir->IntPoint(i);
1323  const DenseMatrix &Jtr_i = Jtr(i);
1324  metric->SetTargetJacobian(Jtr_i);
1325  CalcInverse(Jtr_i, Jrt);
1326  const double weight = ip.weight * Jtr_i.Det();
1327  double weight_m = weight * metric_normal;
1328 
1329  el.CalcDShape(ip, DSh);
1330  Mult(DSh, Jrt, DS);
1331  MultAtB(PMatI, DS, Jpt);
1332 
1333  if (coeff1) { weight_m *= coeff1->Eval(*Tpr, ip); }
1334 
1335  metric->AssembleH(Jpt, DS, weight_m, elmat);
1336 
1337  // TODO: derivatives of adaptivity-based targets.
1338 
1339  if (coeff0)
1340  {
1341  el.CalcShape(ip, shape);
1342  PMatI.MultTranspose(shape, p);
1343  pos0.MultTranspose(shape, p0);
1344  weight_m = weight * lim_normal * coeff0->Eval(*Tpr, ip);
1345  lim_func->Eval_d2(p, p0, d_vals(i), grad_grad);
1346  for (int i = 0; i < dof; i++)
1347  {
1348  const double w_shape_i = weight_m * shape(i);
1349  for (int j = 0; j < dof; j++)
1350  {
1351  const double w = w_shape_i * shape(j);
1352  for (int d1 = 0; d1 < dim; d1++)
1353  {
1354  for (int d2 = 0; d2 < dim; d2++)
1355  {
1356  elmat(d1*dof + i, d2*dof + j) += w * grad_grad(d1, d2);
1357  }
1358  }
1359  }
1360  }
1361  }
1362  }
1363  delete Tpr;
1364 }
1365 
1367 {
1369  metric_normal = 1.0 / metric_normal;
1370  lim_normal = 1.0 / lim_normal;
1371 }
1372 
1373 #ifdef MFEM_USE_MPI
1375 {
1376  double loc[2];
1377  ComputeNormalizationEnergies(x, loc[0], loc[1]);
1378  double rdc[2];
1379  MPI_Allreduce(loc, rdc, 2, MPI_DOUBLE, MPI_SUM, x.ParFESpace()->GetComm());
1380  metric_normal = 1.0 / rdc[0]; lim_normal = 1.0 / rdc[1];
1381 }
1382 #endif
1383 
1385  double &metric_energy,
1386  double &lim_energy)
1387 {
1388  Array<int> vdofs;
1389  Vector x_vals;
1390  const FiniteElementSpace* const fes = x.FESpace();
1391  const FiniteElement *fe = fes->GetFE(0);
1392 
1393  const int dof = fes->GetFE(0)->GetDof(), dim = fes->GetFE(0)->GetDim();
1394 
1395  DSh.SetSize(dof, dim);
1396  Jrt.SetSize(dim);
1397  Jpr.SetSize(dim);
1398  Jpt.SetSize(dim);
1399 
1400  const IntegrationRule *ir = IntRule;
1401  if (!ir)
1402  {
1403  ir = &(IntRules.Get(fe->GetGeomType(), 2*fe->GetOrder() + 3)); // <---
1404  }
1405 
1406  DenseTensor Jtr(dim, dim, ir->GetNPoints());
1407 
1408  metric_energy = 0.0;
1409  lim_energy = 0.0;
1410  for (int i = 0; i < fes->GetNE(); i++)
1411  {
1412  fe = fes->GetFE(i);
1413  fes->GetElementVDofs(i, vdofs);
1414  x.GetSubVector(vdofs, x_vals);
1415  PMatI.UseExternalData(x_vals.GetData(), dof, dim);
1416 
1417  targetC->ComputeElementTargets(i, *fe, *ir, x_vals, Jtr);
1418 
1419  for (int i = 0; i < ir->GetNPoints(); i++)
1420  {
1421  const IntegrationPoint &ip = ir->IntPoint(i);
1422  metric->SetTargetJacobian(Jtr(i));
1423  CalcInverse(Jtr(i), Jrt);
1424  const double weight = ip.weight * Jtr(i).Det();
1425 
1426  fe->CalcDShape(ip, DSh);
1427  MultAtB(PMatI, DSh, Jpr);
1428  Mult(Jpr, Jrt, Jpt);
1429 
1430  metric_energy += weight * metric->EvalW(Jpt);
1431  lim_energy += weight;
1432  }
1433  }
1434 }
1435 
1437  const TargetConstructor &tc,
1438  const Mesh &mesh, GridFunction &metric_gf)
1439 {
1440  const int NE = mesh.GetNE();
1441  const GridFunction &nodes = *mesh.GetNodes();
1442  const int dim = mesh.Dimension();
1443  DenseMatrix Winv(dim), T(dim), A(dim), dshape, pos;
1444  Array<int> pos_dofs, gf_dofs;
1445  DenseTensor W;
1446  Vector posV;
1447 
1448  for (int i = 0; i < NE; i++)
1449  {
1450  const FiniteElement &fe_pos = *nodes.FESpace()->GetFE(i);
1451  const IntegrationRule &ir = metric_gf.FESpace()->GetFE(i)->GetNodes();
1452  const int nsp = ir.GetNPoints(), dof = fe_pos.GetDof();
1453 
1454  dshape.SetSize(dof, dim);
1455  pos.SetSize(dof, dim);
1456  posV.SetDataAndSize(pos.Data(), dof * dim);
1457 
1458  metric_gf.FESpace()->GetElementDofs(i, gf_dofs);
1459  nodes.FESpace()->GetElementVDofs(i, pos_dofs);
1460  nodes.GetSubVector(pos_dofs, posV);
1461 
1462  W.SetSize(dim, dim, nsp);
1463  tc.ComputeElementTargets(i, fe_pos, ir, posV, W);
1464 
1465  for (int j = 0; j < nsp; j++)
1466  {
1467  const DenseMatrix &Wj = W(j);
1468  metric.SetTargetJacobian(Wj);
1469  CalcInverse(Wj, Winv);
1470 
1471  const IntegrationPoint &ip = ir.IntPoint(j);
1472  fe_pos.CalcDShape(ip, dshape);
1473  MultAtB(pos, dshape, A);
1474  Mult(A, Winv, T);
1475 
1476  metric_gf(gf_dofs[j]) = metric.EvalW(T);
1477  }
1478  }
1479 }
1480 
1481 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
Abstract class for Finite Elements.
Definition: fe.hpp:232
MatrixCoefficient * matrix_tspec
Definition: tmop.hpp:652
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:291
virtual void SetParDiscreteTargetSpec(ParGridFunction &tspec)
Definition: tmop.cpp:925
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:277
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
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:258
virtual void SetInitialField(const Vector &init_nodes, const Vector &init_field)=0
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:308
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...
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:2776
void SetSerialMetaInfo(const Mesh &m, const FiniteElementCollection &fec, int num_comp)
Definition: tmop.cpp:1006
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
DenseMatrix & GetPointMat()
Read and write access to the underlying point matrix describing the transformation.
Definition: eltrans.hpp:321
VectorCoefficient * vector_tspec
Definition: tmop.hpp:651
const TargetType target_type
Definition: tmop.hpp:597
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
DenseMatrix PMatO
Definition: tmop.hpp:761
const DenseMatrix * Jtr
Definition: tmop.hpp:26
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:429
void Assemble_ddI2b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:387
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:890
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:572
void Assemble_ddI2b(scalar_t w, scalar_t *A)
DenseMatrix Jrt
Definition: tmop.hpp:761
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 2, using column-major storage.
Definition: invariants.hpp:187
const GridFunction * lim_dist
Definition: tmop.hpp:745
Coefficient * coeff0
Definition: tmop.hpp:743
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:237
void EnableNormalization(const GridFunction &x)
Computes the normalization factors of the metric and limiting integrals using the mesh position given...
Definition: tmop.cpp:1366
const scalar_t * Get_dI3b()
Definition: invariants.hpp:786
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
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:1436
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:172
TMOP_QualityMetric * metric
Definition: tmop.hpp:731
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:254
double Det() const
Definition: densemat.cpp:449
Base class for limiting functions to be used in class TMOP_Integrator.
Definition: tmop.hpp:465
virtual double EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
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.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
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:699
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:626
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:711
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:163
DenseMatrix Jpr
Definition: tmop.hpp:761
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 ...
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:482
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:56
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:472
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:78
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:318
ParFiniteElementSpace * pfes
Definition: tmop.hpp:530
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:562
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
double metric_normal
Definition: tmop.hpp:737
void SetSize(int i, int j, int k)
Definition: densemat.hpp:768
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
virtual void SetSerialDiscreteTargetSpec(GridFunction &tspec)
Definition: tmop.cpp:943
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:494
Coefficient * coeff1
Definition: tmop.hpp:735
Abstract parallel finite element space.
Definition: pfespace.hpp:28
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:186
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:273
const GridFunction * nodes
Definition: tmop.hpp:594
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
Definition: invariants.hpp:435
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:229
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:100
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:166
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:688
Default limiter function in TMOP_Integrator.
Definition: tmop.hpp:486
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:790
void Set(double alpha, const double *A)
Set the matrix to alpha * A, assuming that A has the same dimensions as the matrix and uses column-ma...
Definition: densemat.cpp:533
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:967
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:338
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:365
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:891
Coefficient * scalar_tspec
Definition: tmop.hpp:650
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:518
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:97
void SetTargetJacobian(const DenseMatrix &_Jtr)
Specify the reference-element -&gt; target-element Jacobian matrix for the point of interest.
Definition: tmop.hpp:43
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:68
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:608
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:381
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:311
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:1926
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:408
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:23
Geometry Geometries
Definition: fe.cpp:12638
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:201
virtual ~AdaptivityEvaluator()
Definition: tmop.cpp:1028
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:184
void Assemble_ddI3b(scalar_t w, scalar_t *A)
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:395
MPI_Comm GetComm() const
Definition: pfespace.hpp:239
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:655
void EnableLimiting(const GridFunction &n0, const GridFunction &dist, Coefficient &w0, TMOP_LimiterFunction *lfunc=NULL)
Adds a limiting term to the integrator (general version).
Definition: tmop.cpp:1038
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:452
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
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:323
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:295
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:352
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:225
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:309
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:596
const TargetConstructor * targetC
Definition: tmop.hpp:732
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:461
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...
const scalar_t * Get_dI2b()
Definition: invariants.hpp:219
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:1540
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:411
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:203
const IntegrationRule & GetNodes() const
Definition: fe.hpp:367
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 3, using column-major storage.
Definition: invariants.hpp:743
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:108
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:424
void UpdateTargetSpecification(const Vector &new_x)
Definition: tmop.cpp:960
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:328
int GetNumRootElements()
Return the number of root elements.
Definition: ncmesh.hpp:325
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:442
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: tmop.cpp:1168
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:370
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:200
A class for non-conforming AMR on higher-order hexahedral, prismatic, quadrilateral or triangular mes...
Definition: ncmesh.hpp:100
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:28
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:442
int Dimension() const
Definition: mesh.hpp:744
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1287
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:260
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:192
void ComputeAvgVolume() const
Definition: tmop.cpp:769
int SizeI() const
Definition: densemat.hpp:762
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2197
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:677
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1072
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:96
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1376
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2496
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: tmop.cpp:1262
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:314
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:954
const double eps
Definition: tmop.hpp:308
const scalar_t * Get_dI1()
Definition: invariants.hpp:766
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:431
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:446
const scalar_t * Get_dI1()
Definition: invariants.hpp:207
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:812
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:221
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:44
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:175
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:304
virtual void SetAnalyticTargetSpec(Coefficient *sspec, VectorCoefficient *vspec, MatrixCoefficient *mspec)
Definition: tmop.cpp:882
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:336
const scalar_t * Get_dI1b()
Definition: invariants.hpp:211
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:125
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:474
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:30
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:395
const scalar_t * Get_dI1b()
Definition: invariants.hpp:770
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:185
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:242
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:634
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:403
DenseMatrix PMatI
Definition: tmop.hpp:761
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:411
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:233
int SizeJ() const
Definition: densemat.hpp:763
double a
Definition: lissajous.cpp:41
const scalar_t * Get_dI2b()
Definition: invariants.hpp:778
virtual double Eval(const Vector &x, const Vector &x0, double d) const =0
Returns the limiting function, f(x, x0, d).
Class for integration point with weight.
Definition: intrules.hpp:25
AdaptivityEvaluator * adapt_eval
Definition: tmop.hpp:686
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:347
void Assemble_TProd(scalar_t w, const scalar_t *X, const scalar_t *Y, scalar_t *A)
const scalar_t * Get_dI2()
Definition: invariants.hpp:774
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:580
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:338
void ComputeNormalizationEnergies(const GridFunction &x, double &metric_energy, double &lim_energy)
Definition: tmop.cpp:1384
void SetParMetaInfo(const ParMesh &m, const FiniteElementCollection &fec, int num_comp)
Parallel version of SetSerialMetaInfo.
Definition: tmop.cpp:1017
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:373
int dim
Definition: ex24.cpp:43
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:1074
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:22
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:307
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:796
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.
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:647
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:154
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:527
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:170
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:192
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:128
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:267
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
const scalar_t * Get_dI2()
Definition: invariants.hpp:215
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:663
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:512
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2648
Vector data type.
Definition: vector.hpp:48
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:395
DenseMatrix DSh
Definition: tmop.hpp:761
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:6203
DenseMatrix DS
Definition: tmop.hpp:761
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:379
const GridFunction * nodes0
Definition: tmop.hpp:742
const FiniteElementSpace * tspec_fes
Definition: tmop.hpp:682
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:344
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:466
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:315
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:65
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void Assemble_ddI1b(scalar_t w, scalar_t *A)
Definition: invariants.hpp:830
Base class representing target-matrix construction algorithms for mesh optimization via the target-ma...
Definition: tmop.hpp:570
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:90
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:169
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...
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:725
Class for parallel meshes.
Definition: pmesh.hpp:32
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:376
TMOP_LimiterFunction * lim_func
Definition: tmop.hpp:747
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:618
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:602
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:363
DenseMatrix Jpt
Definition: tmop.hpp:761
void ParEnableNormalization(const ParGridFunction &x)
Definition: tmop.cpp:1374
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:733
DenseMatrix P
Definition: tmop.hpp:761
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:745
virtual void ComputeAtNewPosition(const Vector &new_nodes, Vector &new_field)=0
bool Parallel() const
Definition: tmop.hpp:601
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:725
FiniteElementSpace * fes
Definition: tmop.hpp:525
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:440
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:103
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:740