MFEM  v3.3.2
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "tmop.hpp"
13 #include "linearform.hpp"
14 
15 namespace mfem
16 {
17 
18 // Target-matrix optimization paradigm (TMOP) mesh quality metrics.
19 
20 double TMOP_Metric_001::EvalW(const DenseMatrix &Jpt) const
21 {
22  ie.SetJacobian(Jpt.GetData());
23  return ie.Get_I1();
24 }
25 
27 {
28  ie.SetJacobian(Jpt.GetData());
29  P = ie.Get_dI1();
30 }
31 
33  const DenseMatrix &DS,
34  const double weight,
35  DenseMatrix &A) const
36 {
37  ie.SetJacobian(Jpt.GetData());
39  ie.Assemble_ddI1(weight, A.GetData());
40 }
41 
42 double TMOP_Metric_002::EvalW(const DenseMatrix &Jpt) const
43 {
44  ie.SetJacobian(Jpt.GetData());
45  return 0.5 * ie.Get_I1b() - 1.0;
46 }
47 
49 {
50  ie.SetJacobian(Jpt.GetData());
51  P.Set(0.5, ie.Get_dI1b());
52 }
53 
55  const DenseMatrix &DS,
56  const double weight,
57  DenseMatrix &A) const
58 {
59  ie.SetJacobian(Jpt.GetData());
61  ie.Assemble_ddI1b(0.5*weight, A.GetData());
62 }
63 
64 double TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
65 {
66  // mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
67  ie.SetJacobian(Jpt.GetData());
68  return ie.Get_I1()*(1. + 1./ie.Get_I2()) - 4.0;
69 }
70 
72 {
73  // P = d(I1*(1 + 1/I2)) = (1 + 1/I2) dI1 - I1/I2^2 dI2
74  ie.SetJacobian(Jpt.GetData());
75  const double I2 = ie.Get_I2();
76  Add(1. + 1./I2, ie.Get_dI1(), -ie.Get_I1()/(I2*I2), ie.Get_dI2(), P);
77 }
78 
80  const DenseMatrix &DS,
81  const double weight,
82  DenseMatrix &A) const
83 {
84  // P = d(I1*(1 + 1/I2))
85  // = (1 + 1/I2) dI1 - I1/I2^2 dI2
86  //
87  // dP = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 -
88  // (dI2 x d(I1/I2^2)) - I1/I2^2 ddI2
89  // = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 +
90  // (-1/I2^2) (dI2 x [dI1 - 2 I1/I2 dI2]) - I1/I2^2 ddI2
91  // = (-1/I2^2) (dI1 x dI2 + dI2 x dI1) + (1 + 1/I2) ddI1 +
92  // (2 I1/I2^3) (dI2 x dI2) - I1/I2^2 ddI2
93  ie.SetJacobian(Jpt.GetData());
95  const double c1 = 1./ie.Get_I2();
96  const double c2 = weight*c1*c1;
97  const double c3 = ie.Get_I1()*c2;
98  ie.Assemble_ddI1(weight*(1. + c1), A.GetData());
99  ie.Assemble_ddI2(-c3, A.GetData());
100  ie.Assemble_TProd(-c2, ie.Get_dI1(), ie.Get_dI2(), A.GetData());
101  ie.Assemble_TProd(2*c1*c3, ie.Get_dI2(), A.GetData());
102 }
103 
104 double TMOP_Metric_009::EvalW(const DenseMatrix &Jpt) const
105 {
106  // mu_9 = det(J)*|J-J^{-t}|^2 = I1b * (I2b^2 + 1) - 4 * I2b
107  // = (I1 - 4)*I2b + I1b
108  ie.SetJacobian(Jpt.GetData());
109  return (ie.Get_I1() - 4.0)*ie.Get_I2b() + ie.Get_I1b();
110 }
111 
113 {
114  // mu_9 = (I1 - 4)*I2b + I1b
115  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
116  ie.SetJacobian(Jpt.GetData());
117  Add(ie.Get_I1() - 4.0, ie.Get_dI2b(), ie.Get_I2b(), ie.Get_dI1(), P);
118  P += ie.Get_dI1b();
119 }
120 
122  const DenseMatrix &DS,
123  const double weight,
124  DenseMatrix &A) const
125 {
126  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
127  // dP = dI2b x dI1 + (I1-4)*ddI2b + dI1 x dI2b + I2b*ddI1 + ddI1b
128  // = (dI1 x dI2b + dI2b x dI1) + (I1-4)*ddI2b + I2b*ddI1 + ddI1b
129  ie.SetJacobian(Jpt.GetData());
130  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
131  ie.Assemble_TProd(weight, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
132  ie.Assemble_ddI2b(weight*(ie.Get_I1()-4.0), A.GetData());
133  ie.Assemble_ddI1(weight*ie.Get_I2b(), A.GetData());
134  ie.Assemble_ddI1b(weight, A.GetData());
135 }
136 
137 double TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const
138 {
139  // mu_22 = (0.5*|J|^2 - det(J)) / (det(J) - tau0)
140  // = (0.5*I1 - I2b) / (I2b - tau0)
141  ie.SetJacobian(Jpt.GetData());
142  const double I2b = ie.Get_I2b();
143  return (0.5*ie.Get_I1() - I2b) / (I2b - tau0);
144 }
145 
147 {
148  // mu_22 = (0.5*I1 - I2b) / (I2b - tau0)
149  // P = 1/(I2b - tau0)*(0.5*dI1 - dI2b) - (0.5*I1 - I2b)/(I2b - tau0)^2*dI2b
150  // = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
151  ie.SetJacobian(Jpt.GetData());
152  const double c1 = 1.0/(ie.Get_I2b() - tau0);
153  Add(c1/2, ie.Get_dI1(), (tau0 - ie.Get_I1()/2)*c1*c1, ie.Get_dI2b(), P);
154 }
155 
157  const DenseMatrix &DS,
158  const double weight,
159  DenseMatrix &A) const
160 {
161  // P = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
162  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b) + 0.5/(I2b - tau0)*ddI1
163  // + (dI2b x dz) + z*ddI2b
164  //
165  // z = (tau0 - 0.5*I1)/(I2b - tau0)^2
166  // dz = -0.5/(I2b - tau0)^2*dI1 - 2*(tau0 - 0.5*I1)/(I2b - tau0)^3*dI2b
167  //
168  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b + dI2b x dI1)
169  // -2*z/(I2b - tau0)*(dI2b x dI2b)
170  // +0.5/(I2b - tau0)*ddI1 + z*ddI2b
171  ie.SetJacobian(Jpt.GetData());
172  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
173  const double c1 = 1.0/(ie.Get_I2b() - tau0);
174  const double c2 = weight*c1/2;
175  const double c3 = c1*c2;
176  const double c4 = (2*tau0 - ie.Get_I1())*c3; // weight*z
177  ie.Assemble_TProd(-c3, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
178  ie.Assemble_TProd(-2*c1*c4, ie.Get_dI2b(), A.GetData());
179  ie.Assemble_ddI1(c2, A.GetData());
180  ie.Assemble_ddI2b(c4, A.GetData());
181 }
182 
183 double TMOP_Metric_050::EvalW(const DenseMatrix &Jpt) const
184 {
185  // mu_50 = 0.5*|J^t J|^2/det(J)^2 - 1
186  // = 0.5*(l1^4 + l2^4)/(l1*l2)^2 - 1
187  // = 0.5*((l1/l2)^2 + (l2/l1)^2) - 1 = 0.5*(l1/l2 - l2/l1)^2
188  // = 0.5*(l1/l2 + l2/l1)^2 - 2 = 0.5*I1b^2 - 2
189  ie.SetJacobian(Jpt.GetData());
190  const double I1b = ie.Get_I1b();
191  return 0.5*I1b*I1b - 2.0;
192 }
193 
195 {
196  // mu_50 = 0.5*I1b^2 - 2
197  // P = I1b*dI1b
198  ie.SetJacobian(Jpt.GetData());
199  P.Set(ie.Get_I1b(), ie.Get_dI1b());
200 }
201 
203  const DenseMatrix &DS,
204  const double weight,
205  DenseMatrix &A) const
206 {
207  // P = I1b*dI1b
208  // dP = dI1b x dI1b + I1b*ddI1b
209  ie.SetJacobian(Jpt.GetData());
210  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
211  ie.Assemble_TProd(weight, ie.Get_dI1b(), A.GetData());
212  ie.Assemble_ddI1b(weight*ie.Get_I1b(), A.GetData());
213 }
214 
215 double TMOP_Metric_055::EvalW(const DenseMatrix &Jpt) const
216 {
217  // mu_55 = (det(J) - 1)^2 = (I2b - 1)^2
218  ie.SetJacobian(Jpt.GetData());
219  const double c1 = ie.Get_I2b() - 1.0;
220  return c1*c1;
221 }
222 
224 {
225  // mu_55 = (I2b - 1)^2
226  // P = 2*(I2b - 1)*dI2b
227  ie.SetJacobian(Jpt.GetData());
228  P.Set(2*(ie.Get_I2b() - 1.0), ie.Get_dI2b());
229 }
230 
232  const DenseMatrix &DS,
233  const double weight,
234  DenseMatrix &A) const
235 {
236  // P = 2*(I2b - 1)*dI2b
237  // dP = 2*(dI2b x dI2b) + 2*(I2b - 1)*ddI2b
238  ie.SetJacobian(Jpt.GetData());
239  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
240  ie.Assemble_TProd(2*weight, ie.Get_dI2b(), A.GetData());
241  ie.Assemble_ddI2b(2*weight*(ie.Get_I2b() - 1.0), A.GetData());
242 }
243 
244 double TMOP_Metric_056::EvalW(const DenseMatrix &Jpt) const
245 {
246  // mu_56 = 0.5*(I2b + 1/I2b) - 1
247  ie.SetJacobian(Jpt.GetData());
248  const double I2b = ie.Get_I2b();
249  return 0.5*(I2b + 1.0/I2b) - 1.0;
250 }
251 
253 {
254  // mu_56 = 0.5*(I2b + 1/I2b) - 1
255  // P = 0.5*(1 - 1/I2b^2)*dI2b
256  ie.SetJacobian(Jpt.GetData());
257  P.Set(0.5 - 0.5/ie.Get_I2(), ie.Get_dI2b());
258 }
259 
261  const DenseMatrix &DS,
262  const double weight,
263  DenseMatrix &A) const
264 {
265  // P = 0.5*(1 - 1/I2b^2)*dI2b
266  // dP = (1/I2b^3)*(dI2b x dI2b) + (0.5 - 0.5/I2)*ddI2b
267  ie.SetJacobian(Jpt.GetData());
268  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
269  ie.Assemble_TProd(weight/(ie.Get_I2()*ie.Get_I2b()),
270  ie.Get_dI2b(), A.GetData());
271  ie.Assemble_ddI2b(weight*(0.5 - 0.5/ie.Get_I2()), A.GetData());
272 }
273 
274 double TMOP_Metric_058::EvalW(const DenseMatrix &Jpt) const
275 {
276  // mu_58 = I1b*(I1b - 2)
277  ie.SetJacobian(Jpt.GetData());
278  const double I1b = ie.Get_I1b();
279  return I1b*(I1b - 1.0);
280 }
281 
283 {
284  // mu_58 = I1b*(I1b - 2)
285  // P = (2*I1b - 2)*dI1b
286  ie.SetJacobian(Jpt.GetData());
287  P.Set(2*ie.Get_I1b() - 2.0, ie.Get_dI1b());
288 }
289 
291  const DenseMatrix &DS,
292  const double weight,
293  DenseMatrix &A) const
294 {
295  // P = (2*I1b - 2)*dI1b
296  // dP = 2*(dI1b x dI1b) + (2*I1b - 2)*ddI1b
297  ie.SetJacobian(Jpt.GetData());
298  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
299  ie.Assemble_TProd(2*weight, ie.Get_dI1b(), A.GetData());
300  ie.Assemble_ddI1b(weight*(2*ie.Get_I1b() - 2.0), A.GetData());
301 }
302 
303 double TMOP_Metric_077::EvalW(const DenseMatrix &Jpt) const
304 {
305  ie.SetJacobian(Jpt.GetData());
306  const double I2 = ie.Get_I2b();
307  return 0.5*(I2*I2 + 1./(I2*I2) - 2.);
308 }
309 
311 {
312  // Using I2b^2 = I2.
313  // dmu77_dJ = 1/2 (1 - 1/I2^2) dI2_dJ.
314  ie.SetJacobian(Jpt.GetData());
315  const double I2 = ie.Get_I2();
316  P.Set(0.5 * (1.0 - 1.0 / (I2 * I2)), ie.Get_dI2());
317 }
318 
320  const DenseMatrix &DS,
321  const double weight,
322  DenseMatrix &A) const
323 {
324  ie.SetJacobian(Jpt.GetData());
325  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
326  const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
327  ie.Assemble_ddI2(weight*0.5*(1.0 - I2inv_sq), A.GetData());
328  ie.Assemble_TProd(weight * I2inv_sq / I2, ie.Get_dI2(), A.GetData());
329 }
330 
331 double TMOP_Metric_211::EvalW(const DenseMatrix &Jpt) const
332 {
333  // mu_211 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
334  // = (I2b - 1)^2 - I2b + sqrt(I2b^2 + eps)
335  ie.SetJacobian(Jpt.GetData());
336  const double I2b = ie.Get_I2b();
337  return (I2b - 1.0)*(I2b - 1.0) - I2b + std::sqrt(I2b*I2b + eps);
338 }
339 
341 {
342  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
343 }
344 
346  const DenseMatrix &DS,
347  const double weight,
348  DenseMatrix &A) const
349 {
350  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
351 }
352 
353 double TMOP_Metric_252::EvalW(const DenseMatrix &Jpt) const
354 {
355  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0).
356  ie.SetJacobian(Jpt.GetData());
357  const double I2b = ie.Get_I2b();
358  return 0.5*(I2b - 1.0)*(I2b - 1.0)/(I2b - tau0);
359 }
360 
362 {
363  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
364  // P = (c - 0.5*c*c ) * dI2b
365  //
366  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
367  ie.SetJacobian(Jpt.GetData());
368  const double I2b = ie.Get_I2b();
369  const double c = (I2b - 1.0)/(I2b - tau0);
370  P.Set(c - 0.5*c*c, ie.Get_dI2b());
371 }
372 
374  const DenseMatrix &DS,
375  const double weight,
376  DenseMatrix &A) const
377 {
378  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
379  //
380  // P = (c - 0.5*c*c ) * dI2b
381  // dP = (1 - c)^2/(I2b - tau0)*(dI2b x dI2b) + (c - 0.5*c*c)*ddI2b
382  ie.SetJacobian(Jpt.GetData());
383  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
384  const double I2b = ie.Get_I2b();
385  const double c0 = 1.0/(I2b - tau0);
386  const double c = c0*(I2b - 1.0);
387  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI2b(), A.GetData());
388  ie.Assemble_ddI2b(weight*(c - 0.5*c*c), A.GetData());
389 }
390 
391 double TMOP_Metric_301::EvalW(const DenseMatrix &Jpt) const
392 {
393  ie.SetJacobian(Jpt.GetData());
394  return std::sqrt(ie.Get_I1b()*ie.Get_I2b())/3. - 1.;
395 }
396 
398 {
399  // W = (1/3)*sqrt(I1b*I2b) - 1
400  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
401  ie.SetJacobian(Jpt.GetData());
402  const double a = 1./(6.*std::sqrt(ie.Get_I1b()*ie.Get_I2b()));
403  Add(a*ie.Get_I2b(), ie.Get_dI1b(), a*ie.Get_I1b(), ie.Get_dI2b(), P);
404 }
405 
407  const DenseMatrix &DS,
408  const double weight,
409  DenseMatrix &A) const
410 {
411  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
412  // dW = (1/6)*[z2*dI1b + z1*dI2b], z1 = sqrt(I1b/I2b), z2 = sqrt(I2b/I1b)
413  // ddW = (1/6)*[dI1b x dz2 + z2*ddI1b + dI2b x dz1 + z1*ddI2b]
414  //
415  // dz1 = (1/2)*sqrt(I2b/I1b) [ (1/I2b)*dI1b + (I1b/(I2b*I2b))*dI2b ]
416  // = (1/2)/sqrt(I1b*I2b) [ dI1b + (I1b/I2b)*dI2b ]
417  // dz2 = (1/2)/sqrt(I1b*I2b) [ (I2b/I1b)*dI1b + dI2b ]
418  //
419  // dI1b x dz2 + dI2b x dz1 =
420  // (1/2)/sqrt(I1b*I2b) dI1b x [ (I2b/I1b)*dI1b + dI2b ] +
421  // (1/2)/sqrt(I1b*I2b) dI2b x [ dI1b + (I1b/I2b)*dI2b ] =
422  // (1/2)/sqrt(I1b*I2b) [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] x
423  // [sqrt(I2b/I1b)*dI1b + sqrt(I1b/I2b)*dI2b] =
424  // (1/2)/sqrt(I1b*I2b) [ 6*dW x 6*dW ] =
425  // (1/2)*(I1b*I2b)^{-3/2} (I2b*dI1b + I1b*dI2b) x (I2b*dI1b + I1b*dI2b)
426  //
427  // z1 = I1b/sqrt(I1b*I2b), z2 = I2b/sqrt(I1b*I2b)
428 
429  ie.SetJacobian(Jpt.GetData());
430  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
431  double d_I1b_I2b_data[9];
432  DenseMatrix d_I1b_I2b(d_I1b_I2b_data, 3, 3);
433  Add(ie.Get_I2b(), ie.Get_dI1b(), ie.Get_I1b(), ie.Get_dI2b(), d_I1b_I2b);
434  const double I1b_I2b = ie.Get_I1b()*ie.Get_I2b();
435  const double a = weight/(6*std::sqrt(I1b_I2b));
436  ie.Assemble_ddI1b(a*ie.Get_I2b(), A.GetData());
437  ie.Assemble_ddI2b(a*ie.Get_I1b(), A.GetData());
438  ie.Assemble_TProd(a/(2*I1b_I2b), d_I1b_I2b_data, A.GetData());
439 }
440 
441 double TMOP_Metric_302::EvalW(const DenseMatrix &Jpt) const
442 {
443  // mu_2 = |J|^2 |J^{-1}|^2 / 9 - 1
444  // = (l1^2 + l2^2 + l3^3)*(l1^{-2} + l2^{-2} + l3^{-2}) / 9 - 1
445  // = I1*(l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/l1^2/l2^2/l3^2/9 - 1
446  // = I1*I2/det(J)^2/9 - 1 = I1b*I2b/9-1
447  ie.SetJacobian(Jpt.GetData());
448  return ie.Get_I1b()*ie.Get_I2b()/9. - 1.;
449 }
450 
452 {
453  // mu_2 = I1b*I2b/9-1
454  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
455  ie.SetJacobian(Jpt.GetData());
456  Add(ie.Get_I1b()/9, ie.Get_dI2b(), ie.Get_I2b()/9, ie.Get_dI1b(), P);
457 }
458 
460  const DenseMatrix &DS,
461  const double weight,
462  DenseMatrix &A) const
463 {
464  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
465  // dP = (dI2b x dI1b)/9 + (I1b/9)*ddI2b + (dI1b x dI2b)/9 + (I2b/9)*ddI1b
466  // = (dI2b x dI1b + dI1b x dI2b)/9 + (I1b/9)*ddI2b + (I2b/9)*ddI1b
467  ie.SetJacobian(Jpt.GetData());
468  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
469  const double c1 = weight/9;
470  ie.Assemble_TProd(c1, ie.Get_dI1b(), ie.Get_dI2b(), A.GetData());
471  ie.Assemble_ddI2b(c1*ie.Get_I1b(), A.GetData());
472  ie.Assemble_ddI1b(c1*ie.Get_I2b(), A.GetData());
473 }
474 
475 double TMOP_Metric_303::EvalW(const DenseMatrix &Jpt) const
476 {
477  ie.SetJacobian(Jpt.GetData());
478  return ie.Get_I1b()/3.0 - 1.0;
479 }
480 
482 {
483  ie.SetJacobian(Jpt.GetData());
484  P.Set(1./3., ie.Get_dI1b());
485 }
486 
488  const DenseMatrix &DS,
489  const double weight,
490  DenseMatrix &A) const
491 {
492  ie.SetJacobian(Jpt.GetData());
493  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
494  ie.Assemble_ddI1b(weight/3., A.GetData());
495 }
496 
497 double TMOP_Metric_315::EvalW(const DenseMatrix &Jpt) const
498 {
499  // mu_315 = mu_15_3D = (det(J) - 1)^2
500  ie.SetJacobian(Jpt.GetData());
501  const double c1 = ie.Get_I3b() - 1.0;
502  return c1*c1;
503 }
504 
506 {
507  // mu_315 = (I3b - 1)^2
508  // P = 2*(I3b - 1)*dI3b
509  ie.SetJacobian(Jpt.GetData());
510  P.Set(2*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
511 }
512 
514  const DenseMatrix &DS,
515  const double weight,
516  DenseMatrix &A) const
517 {
518  // P = 2*(I3b - 1)*dI3b
519  // dP = 2*(dI3b x dI3b) + 2*(I3b - 1)*ddI3b
520  ie.SetJacobian(Jpt.GetData());
521  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
522  ie.Assemble_TProd(2*weight, ie.Get_dI3b(), A.GetData());
523  ie.Assemble_ddI3b(2*weight*(ie.Get_I3b() - 1.0), A.GetData());
524 }
525 
526 double TMOP_Metric_316::EvalW(const DenseMatrix &Jpt) const
527 {
528  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
529  ie.SetJacobian(Jpt.GetData());
530  const double I3b = ie.Get_I3b();
531  return 0.5*(I3b + 1.0/I3b) - 1.0;
532 }
533 
535 {
536  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
537  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
538  ie.SetJacobian(Jpt.GetData());
539  P.Set(0.5 - 0.5/ie.Get_I3(), ie.Get_dI3b());
540 }
541 
543  const DenseMatrix &DS,
544  const double weight,
545  DenseMatrix &A) const
546 {
547  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
548  // dP = (1/I3b^3)*(dI3b x dI3b) + (0.5 - 0.5/I3)*ddI3b
549  ie.SetJacobian(Jpt.GetData());
550  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
551  ie.Assemble_TProd(weight/(ie.Get_I3()*ie.Get_I3b()),
552  ie.Get_dI3b(), A.GetData());
553  ie.Assemble_ddI3b(weight*(0.5 - 0.5/ie.Get_I3()), A.GetData());
554 }
555 
556 double TMOP_Metric_321::EvalW(const DenseMatrix &Jpt) const
557 {
558  // mu_321 = mu_21_3D = |J - J^{-t}|^2
559  // = |J|^2 + |J^{-1}|^2 - 6
560  // = |J|^2 + (l1^{-2} + l2^{-2} + l3^{-2}) - 6
561  // = |J|^2 + (l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/det(J)^2 - 6
562  // = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
563  ie.SetJacobian(Jpt.GetData());
564  return ie.Get_I1() + ie.Get_I2()/ie.Get_I3() - 6.0;
565 }
566 
568 {
569  // mu_321 = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
570  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
571  ie.SetJacobian(Jpt.GetData());
572  const double I3 = ie.Get_I3();
573  Add(1.0/I3, ie.Get_dI2(),
574  -2*ie.Get_I2()/(I3*ie.Get_I3b()), ie.Get_dI3b(), P);
575  P += ie.Get_dI1();
576 }
577 
579  const DenseMatrix &DS,
580  const double weight,
581  DenseMatrix &A) const
582 {
583  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
584  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b) + (1/I3)*ddI2 + (dI3b x dz) + z*ddI3b
585  //
586  // z = -2*I2/I3b^3
587  // dz = (-2/I3b^3)*dI2 + (2*I2)*(3/I3b^4)*dI3b
588  //
589  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b + dI3b x dI2) + (1/I3)*ddI2
590  // + (6*I2/I3b^4)*(dI3b x dI3b) + (-2*I2/I3b^3)*ddI3b
591  ie.SetJacobian(Jpt.GetData());
592  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
593  const double c0 = 1.0/ie.Get_I3b();
594  const double c1 = weight*c0*c0;
595  const double c2 = -2*c0*c1;
596  const double c3 = c2*ie.Get_I2();
597  ie.Assemble_ddI1(weight, A.GetData());
598  ie.Assemble_ddI2(c1, A.GetData());
599  ie.Assemble_ddI3b(c3, A.GetData());
600  ie.Assemble_TProd(c2, ie.Get_dI2(), ie.Get_dI3b(), A.GetData());
601  ie.Assemble_TProd(-3*c0*c3, ie.Get_dI3b(), A.GetData());
602 }
603 
604 double TMOP_Metric_352::EvalW(const DenseMatrix &Jpt) const
605 {
606  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
607  ie.SetJacobian(Jpt.GetData());
608  const double I3b = ie.Get_I3b();
609  return 0.5*(I3b - 1.0)*(I3b - 1.0)/(I3b - tau0);
610 }
611 
613 {
614  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
615  // P = (I3b - 1)/(I3b - tau0)*dI3b + 0.5*(I3b - 1)^2*(-1/(I3b - tau0)^2)*dI3b
616  // = [ (I3b - 1)/(I3b - tau0) - 0.5*(I3b - 1)^2/(I3b - tau0)^2 ] * dI3b
617  // = (c - 0.5*c*c) * dI3b
618  ie.SetJacobian(Jpt.GetData());
619  const double I3b = ie.Get_I3b();
620  const double c = (I3b - 1.0)/(I3b - tau0);
621  P.Set(c - 0.5*c*c, ie.Get_dI3b());
622 }
623 
625  const DenseMatrix &DS,
626  const double weight,
627  DenseMatrix &A) const
628 {
629  // c = (I3b - 1)/(I3b - tau0)
630  //
631  // P = (c - 0.5*c*c) * dI3b
632  // dP = (1 - c)*(dI3b x dc) + (c - 0.5*c*c)*ddI3b
633  //
634  // dc = 1/(I3b - tau0)*dI3b - (I3b - 1)/(I3b - tau)^2*dI3b =
635  // = (1 - c)/(I3b - tau0)*dI3b
636  //
637  // dP = (1 - c)^2/(I3b - tau0)*(dI3b x dI3b) + (c - 0.5*c*c)*ddI3b
638  ie.SetJacobian(Jpt.GetData());
639  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
640  const double I3b = ie.Get_I3b();
641  const double c0 = 1.0/(I3b - tau0);
642  const double c = c0*(I3b - 1.0);
643  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI3b(), A.GetData());
644  ie.Assemble_ddI3b(weight*(c - 0.5*c*c), A.GetData());
645 }
646 
647 
649 {
650  MFEM_VERIFY(nodes, "Nodes are not given!");
651  MFEM_ASSERT(avg_volume == 0.0, "the average volume is already computed!");
652 
653  Mesh *mesh = nodes->FESpace()->GetMesh();
654  const int NE = mesh->GetNE();
656  double volume = 0.0;
657 
658  for (int i = 0; i < NE; i++)
659  {
660  mesh->GetElementTransformation(i, *nodes, &Tr);
661  const IntegrationRule &ir =
662  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr.OrderJ());
663  for (int j = 0; j < ir.GetNPoints(); j++)
664  {
665  const IntegrationPoint &ip = ir.IntPoint(j);
666  Tr.SetIntPoint(&ip);
667  volume += ip.weight * Tr.Weight();
668  }
669  }
670  if (!Parallel())
671  {
672  avg_volume = volume / NE;
673  }
674 #ifdef MFEM_USE_MPI
675  else
676  {
677  double area_NE[4];
678  area_NE[0] = volume; area_NE[1] = NE;
679  MPI_Allreduce(area_NE, area_NE + 2, 2, MPI_DOUBLE, MPI_SUM, comm);
680  avg_volume = area_NE[2] / area_NE[3];
681  }
682 #endif
683 }
684 
685 // virtual method
687  const IntegrationRule &ir,
688  DenseTensor &Jtr) const
689 {
690  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
691 
693  nodes->FESpace()->GetFE(e_id) : NULL;
694  const DenseMatrix &Wideal =
696  MFEM_ASSERT(Wideal.Height() == Jtr.SizeI(), "");
697  MFEM_ASSERT(Wideal.Width() == Jtr.SizeJ(), "");
698 
699  switch (target_type)
700  {
702  {
703  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = Wideal; }
704  break;
705  }
707  {
708  if (avg_volume == 0.0) { ComputeAvgVolume(); }
709  DenseMatrix W(Wideal.Height());
710  W.Set(std::pow(volume_scale * avg_volume / Wideal.Det(),
711  1./W.Height()), Wideal);
712  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = W; }
713  break;
714  }
717  {
718  const int dim = nfe->GetDim(), dof = nfe->GetDof();
719  MFEM_ASSERT(dim == nodes->FESpace()->GetVDim(), "");
720  DenseMatrix dshape(dof, dim), pos(dof, dim);
721  Array<int> xdofs(dof * dim);
722  Vector posV(pos.Data(), dof * dim);
723  double detW;
724 
725  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { detW = Wideal.Det(); }
726  nodes->FESpace()->GetElementVDofs(e_id, xdofs);
727  nodes->GetSubVector(xdofs, posV);
728  for (int i = 0; i < ir.GetNPoints(); i++)
729  {
730  nfe->CalcDShape(ir.IntPoint(i), dshape);
731  MultAtB(pos, dshape, Jtr(i));
733  {
734  const double det = Jtr(i).Det();
735  MFEM_VERIFY(det > 0.0, "Initial mesh is inverted!");
736  Jtr(i).Set(std::pow(det / detW, 1./dim), Wideal);
737  }
738  }
739  break;
740  }
741  default:
742  MFEM_ABORT("invalid target type!");
743  }
744 }
745 
746 
749  const Vector &elfun)
750 {
751  int dof = el.GetDof(), dim = el.GetDim();
752  double energy;
753 
754  DSh.SetSize(dof, dim);
755  Jrt.SetSize(dim);
756  Jpr.SetSize(dim);
757  Jpt.SetSize(dim);
758  PMatI.UseExternalData(elfun.GetData(), dof, dim);
759 
760  const IntegrationRule *ir = IntRule;
761  if (!ir)
762  {
763  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
764  }
765 
766  energy = 0.0;
767  DenseTensor Jtr(dim, dim, ir->GetNPoints());
768  targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);
769 
770  // Limited case.
771  Vector shape, p, p0;
772  DenseMatrix pos0;
773  if (coeff0)
774  {
775  shape.SetSize(dof);
776  p.SetSize(dim);
777  p0.SetSize(dim);
778  pos0.SetSize(dof, dim);
779  Vector pos0V(pos0.Data(), dof * dim);
780  Array<int> pos_dofs;
781  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
782  nodes0->GetSubVector(pos_dofs, pos0V);
783  }
784 
785  // Define ref->physical transformation, when a Coefficient is specified.
786  IsoparametricTransformation *Tpr = NULL;
787  if (coeff1 || coeff0)
788  {
789  Tpr = new IsoparametricTransformation;
790  Tpr->SetFE(&el);
791  Tpr->ElementNo = T.ElementNo;
792  Tpr->Attribute = T.Attribute;
793  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
794  }
795 
796  for (int i = 0; i < ir->GetNPoints(); i++)
797  {
798  const IntegrationPoint &ip = ir->IntPoint(i);
799  const DenseMatrix &Jtr_i = Jtr(i);
800  metric->SetTargetJacobian(Jtr_i);
801  CalcInverse(Jtr_i, Jrt);
802  const double weight = ip.weight * Jtr_i.Det();
803 
804  el.CalcDShape(ip, DSh);
805  MultAtB(PMatI, DSh, Jpr);
806  Mult(Jpr, Jrt, Jpt);
807 
808  double val = metric->EvalW(Jpt);
809  if (coeff1) { val *= coeff1->Eval(*Tpr, ip); }
810 
811  if (coeff0)
812  {
813  el.CalcShape(ip, shape);
814  PMatI.MultTranspose(shape, p);
815  pos0.MultTranspose(shape, p0);
816  val += 0.5*p.DistanceSquaredTo(p0)*coeff0->Eval(*Tpr, ip);
817  }
818  energy += weight * val;
819  }
820  delete Tpr;
821  return energy;
822 }
823 
826  const Vector &elfun, Vector &elvect)
827 {
828  int dof = el.GetDof(), dim = el.GetDim();
829 
830  DSh.SetSize(dof, dim);
831  DS.SetSize(dof, dim);
832  Jrt.SetSize(dim);
833  Jpt.SetSize(dim);
834  P.SetSize(dim);
835  PMatI.UseExternalData(elfun.GetData(), dof, dim);
836  elvect.SetSize(dof*dim);
837  PMatO.UseExternalData(elvect.GetData(), dof, dim);
838 
839  const IntegrationRule *ir = IntRule;
840  if (!ir)
841  {
842  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
843  }
844 
845  elvect = 0.0;
846  DenseTensor Jtr(dim, dim, ir->GetNPoints());
847  targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);
848 
849  // Limited case.
850  DenseMatrix pos0;
851  Vector shape, p, p0;
852  if (coeff0)
853  {
854  shape.SetSize(dof);
855  p.SetSize(dim);
856  p0.SetSize(dim);
857  pos0.SetSize(dof, dim);
858  Vector pos0V(pos0.Data(), dof * dim);
859  Array<int> pos_dofs;
860  nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
861  nodes0->GetSubVector(pos_dofs, pos0V);
862  }
863 
864  // Define ref->physical transformation, when a Coefficient is specified.
865  IsoparametricTransformation *Tpr = NULL;
866  if (coeff1 || coeff0)
867  {
868  Tpr = new IsoparametricTransformation;
869  Tpr->SetFE(&el);
870  Tpr->ElementNo = T.ElementNo;
871  Tpr->Attribute = T.Attribute;
872  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
873  }
874 
875  for (int i = 0; i < ir->GetNPoints(); i++)
876  {
877  const IntegrationPoint &ip = ir->IntPoint(i);
878  const DenseMatrix &Jtr_i = Jtr(i);
879  metric->SetTargetJacobian(Jtr_i);
880  CalcInverse(Jtr_i, Jrt);
881  const double weight = ip.weight * Jtr_i.Det();
882  double weight_m = weight;
883 
884  el.CalcDShape(ip, DSh);
885  Mult(DSh, Jrt, DS);
886  MultAtB(PMatI, DS, Jpt);
887 
888  metric->EvalP(Jpt, P);
889 
890  if (coeff1) { weight_m *= coeff1->Eval(*Tpr, ip); }
891 
892  P *= weight_m;
893  AddMultABt(DS, P, PMatO);
894 
895  if (coeff0)
896  {
897  el.CalcShape(ip, shape);
898  PMatI.MultTranspose(shape, p);
899  pos0.MultTranspose(shape, p0);
900  weight_m = weight * coeff0->Eval(*Tpr, ip);
901  subtract(weight_m, p, p0, p); // p = weight_m * (p - p0)
902  AddMultVWt(shape, p, PMatO);
903  }
904  }
905  delete Tpr;
906 }
907 
908 
911  const Vector &elfun,
912  DenseMatrix &elmat)
913 {
914  int dof = el.GetDof(), dim = el.GetDim();
915 
916  DSh.SetSize(dof, dim);
917  DS.SetSize(dof, dim);
918  Jrt.SetSize(dim);
919  Jpt.SetSize(dim);
920  PMatI.UseExternalData(elfun.GetData(), dof, dim);
921  elmat.SetSize(dof*dim);
922 
923  const IntegrationRule *ir = IntRule;
924  if (!ir)
925  {
926  ir = &(IntRules.Get(el.GetGeomType(), 2*el.GetOrder() + 3)); // <---
927  }
928 
929  elmat = 0.0;
930  DenseTensor Jtr(dim, dim, ir->GetNPoints());
931  targetC->ComputeElementTargets(T.ElementNo, el, *ir, Jtr);
932 
933  // Define ref->physical transformation, when a Coefficient is specified.
934  IsoparametricTransformation *Tpr = NULL;
935  if (coeff1 || coeff0)
936  {
937  Tpr = new IsoparametricTransformation;
938  Tpr->SetFE(&el);
939  Tpr->ElementNo = T.ElementNo;
940  Tpr->Attribute = T.Attribute;
941  Tpr->GetPointMat().Transpose(PMatI);
942  }
943 
944  Vector shape(coeff0 ? dof : 0);
945 
946  for (int i = 0; i < ir->GetNPoints(); i++)
947  {
948  const IntegrationPoint &ip = ir->IntPoint(i);
949  const DenseMatrix &Jtr_i = Jtr(i);
950  metric->SetTargetJacobian(Jtr_i);
951  CalcInverse(Jtr_i, Jrt);
952  const double weight = ip.weight * Jtr_i.Det();
953  double weight_m = weight;
954 
955  el.CalcDShape(ip, DSh);
956  Mult(DSh, Jrt, DS);
957  MultAtB(PMatI, DS, Jpt);
958 
959  if (coeff1) { weight_m *= coeff1->Eval(*Tpr, ip); }
960 
961  metric->AssembleH(Jpt, DS, weight_m, elmat);
962 
963  if (coeff0)
964  {
965  el.CalcShape(ip, shape);
966  weight_m = weight * coeff0->Eval(*Tpr, ip);
967  for (int i = 0; i < dof; i++)
968  {
969  const double w_shape_i = weight_m * shape(i);
970  for (int j = 0; j <= i; j++)
971  {
972  const double a = w_shape_i * shape(j);
973  for (int d = 0; d < dim; d++)
974  {
975  elmat(i+d*dof, j+d*dof) += a;
976  if (i != j) { elmat(j+d*dof, i+d*dof) += a; }
977  }
978  }
979  }
980  }
981  }
982  delete Tpr;
983 }
984 
985 
987  const TargetConstructor &tc,
988  const Mesh &mesh, GridFunction &metric_gf)
989 {
990  const int NE = mesh.GetNE();
991  const GridFunction &nodes = *mesh.GetNodes();
992  const int dim = mesh.Dimension();
993  DenseMatrix Winv(dim), T(dim), A(dim), dshape, pos;
994  Array<int> pos_dofs, gf_dofs;
995  DenseTensor W;
996  Vector posV;
997 
998  for (int i = 0; i < NE; i++)
999  {
1000  const FiniteElement &fe_pos = *nodes.FESpace()->GetFE(i);
1001  const IntegrationRule &ir = metric_gf.FESpace()->GetFE(i)->GetNodes();
1002  const int nsp = ir.GetNPoints(), dof = fe_pos.GetDof();
1003 
1004  W.SetSize(dim, dim, nsp);
1005  tc.ComputeElementTargets(i, fe_pos, ir, W);
1006 
1007  dshape.SetSize(dof, dim);
1008  pos.SetSize(dof, dim);
1009  posV.SetDataAndSize(pos.Data(), dof * dim);
1010 
1011  metric_gf.FESpace()->GetElementDofs(i, gf_dofs);
1012  nodes.FESpace()->GetElementVDofs(i, pos_dofs);
1013  nodes.GetSubVector(pos_dofs, posV);
1014 
1015  for (int j = 0; j < nsp; j++)
1016  {
1017  const DenseMatrix &Wj = W(j);
1018  metric.SetTargetJacobian(Wj);
1019  CalcInverse(Wj, Winv);
1020 
1021  const IntegrationPoint &ip = ir.IntPoint(j);
1022  fe_pos.CalcDShape(ip, dshape);
1023  MultAtB(pos, dshape, A);
1024  Mult(A, Winv, T);
1025 
1026  metric_gf(gf_dofs[j]) = metric.EvalW(T);
1027  }
1028  }
1029 }
1030 
1031 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
Abstract class for Finite Elements.
Definition: fe.hpp:47
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:232
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:156
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:32
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:137
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:116
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:3764
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
DenseMatrix & GetPointMat()
Read and write access to the underlying point matrix describing the transformation.
Definition: eltrans.hpp:312
const TargetType target_type
Definition: tmop.hpp:435
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
DenseMatrix PMatO
Definition: tmop.hpp:513
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:370
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:861
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:451
void Assemble_ddI2b(scalar_t w, scalar_t *A)
DenseMatrix Jrt
Definition: tmop.hpp:513
void SetDerivativeMatrix(int height, const scalar_t *Deriv)
The Deriv matrix is dof x 2, using column-major storage.
Definition: invariants.hpp:187
Coefficient * coeff0
Definition: tmop.hpp:501
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:178
const scalar_t * Get_dI3b()
Definition: invariants.hpp:786
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
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:986
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:133
TMOP_QualityMetric * metric
Definition: tmop.hpp:493
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:195
double Det() const
Definition: densemat.cpp:435
virtual double EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
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:578
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:505
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:42
DenseMatrix Jpr
Definition: tmop.hpp:513
virtual void ComputeElementTargets(int e_id, const FiniteElement &fe, const IntegrationRule &ir, DenseTensor &Jtr) const
Given an element and quadrature rule, computes ref-&gt;target transformation Jacobians for each quadratu...
Definition: tmop.cpp:686
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:361
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:462
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:79
int GetOrder() const
Returns the order of the finite element.
Definition: fe.hpp:125
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:441
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void SetSize(int i, int j, int k)
Definition: densemat.hpp:662
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:614
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:373
Coefficient * coeff1
Definition: tmop.hpp:497
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:127
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:214
const GridFunction * nodes
Definition: tmop.hpp:432
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:94
double * GetData() const
Definition: vector.hpp:121
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:567
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:519
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:244
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:397
const DenseMatrix & GetGeomToPerfGeomJac(int GeomType) const
Definition: geom.hpp:88
void SetTargetJacobian(const DenseMatrix &_Jtr)
Specify the reference-element -&gt; target-element Jacobian matrix for the point of interest.
Definition: tmop.hpp: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:487
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:260
void Add(const DenseMatrix &A, const DenseMatrix &B, double alpha, DenseMatrix &C)
C = A + alpha*B.
Definition: densemat.cpp:2942
Abstract class for local mesh quality metrics in the target-matrix optimization paradigm (TMOP) by P...
Definition: tmop.hpp:24
Geometry Geometries
Definition: geom.cpp:759
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:225
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:194
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:336
int dim
Definition: ex3.cpp:47
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:534
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:331
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
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:202
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:136
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:231
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:104
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:250
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:475
const TargetConstructor * targetC
Definition: tmop.hpp:494
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:340
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:1072
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:290
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:144
const IntegrationRule & GetNodes() const
Definition: fe.hpp:167
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:303
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:269
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, Vector &elvect)
Perform the local action of the NonlinearFormIntegrator.
Definition: tmop.cpp:824
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:153
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:79
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:26
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:309
int Dimension() const
Definition: mesh.hpp:641
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:71
void ComputeAvgVolume() const
Definition: tmop.cpp:648
int SizeI() const
Definition: densemat.hpp:658
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3156
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:119
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:556
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:92
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:2401
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:3485
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &T, const Vector &elfun, DenseMatrix &elmat)
Assemble the local gradient matrix.
Definition: tmop.cpp:909
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:122
void Assemble_ddI2(scalar_t w, scalar_t *A)
Definition: invariants.hpp:954
const double eps
Definition: tmop.hpp:249
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:310
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:387
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:390
const scalar_t * Get_dI1()
Definition: invariants.hpp:207
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:162
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:54
double DistanceSquaredTo(const double *p) const
Compute the square of the Euclidean distance to another vector.
Definition: vector.hpp:403
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:183
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:215
const scalar_t * Get_dI1b()
Definition: invariants.hpp:211
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:94
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:353
const IntegrationRule * IntRule
Definition: nonlininteg.hpp:29
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:274
const scalar_t * Get_dI1b()
Definition: invariants.hpp:770
int GetElementBaseGeometry(int i=0) const
Definition: mesh.hpp:684
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:64
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:121
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:513
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:282
DenseMatrix PMatI
Definition: tmop.hpp:513
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:352
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:112
int SizeJ() const
Definition: densemat.hpp:659
const scalar_t * Get_dI2b()
Definition: invariants.hpp:778
Class for integration point with weight.
Definition: intrules.hpp:25
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:288
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:459
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:255
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:252
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:747
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:20
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:298
void Assemble_ddI1(scalar_t w, scalar_t *A)
Definition: invariants.hpp:796
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1171
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:526
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:95
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:406
InvariantsEvaluator2D< double > ie
Definition: tmop.hpp:111
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:146
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
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:542
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:391
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:3637
Vector data type.
Definition: vector.hpp:41
DenseMatrix DSh
Definition: tmop.hpp:513
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:5401
DenseMatrix DS
Definition: tmop.hpp:513
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:320
const GridFunction * nodes0
Definition: tmop.hpp:500
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:223
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:345
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:194
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:64
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:410
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:86
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:48
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:634
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:343
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:497
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:481
InvariantsEvaluator3D< double > ie
Definition: tmop.hpp:304
DenseMatrix Jpt
Definition: tmop.hpp:513
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Definition: tmop.cpp:612
DenseMatrix P
Definition: tmop.hpp:513
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:624
bool Parallel() const
Definition: tmop.hpp:439
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt).
Definition: tmop.cpp:604
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:319
void SetJacobian(const scalar_t *Jac)
The Jacobian should use column-major storage.
Definition: invariants.hpp:740