MFEM  v4.6.0
Finite element discretization library
tmop.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "tmop.hpp"
13 #include "linearform.hpp"
14 #include "pgridfunc.hpp"
15 #include "tmop_tools.hpp"
16 #include "../general/forall.hpp"
17 
18 namespace mfem
19 {
20 
21 // Target-matrix optimization paradigm (TMOP) mesh quality metrics.
22 
23 double TMOP_Combo_QualityMetric::EvalWMatrixForm(const DenseMatrix &Jpt) const
24 {
25  double metric = 0.;
26  for (int i = 0; i < tmop_q_arr.Size(); i++)
27  {
28  metric += wt_arr[i]*tmop_q_arr[i]->EvalWMatrixForm(Jpt);
29  }
30  return metric;
31 }
32 
33 double TMOP_Combo_QualityMetric::EvalW(const DenseMatrix &Jpt) const
34 {
35  double metric = 0.;
36  for (int i = 0; i < tmop_q_arr.Size(); i++)
37  {
38  metric += wt_arr[i]*tmop_q_arr[i]->EvalW(Jpt);
39  }
40  return metric;
41 }
42 
43 void TMOP_Combo_QualityMetric::EvalP(const DenseMatrix &Jpt,
44  DenseMatrix &P) const
45 {
46  DenseMatrix Pt(P.Size());
47  P = 0.0;
48  for (int i = 0; i < tmop_q_arr.Size(); i++)
49  {
50  tmop_q_arr[i]->EvalP(Jpt, Pt);
51  P.Add(wt_arr[i], Pt);
52  }
53 }
54 
55 void TMOP_Combo_QualityMetric::AssembleH(const DenseMatrix &Jpt,
56  const DenseMatrix &DS,
57  const double weight,
58  DenseMatrix &A) const
59 {
60  DenseMatrix At(A.Size());
61  for (int i = 0; i < tmop_q_arr.Size(); i++)
62  {
63  At = 0.0;
64  tmop_q_arr[i]->AssembleH(Jpt, DS, weight * wt_arr[i], At);
65  A += At;
66  }
67 }
68 
69 void TMOP_Combo_QualityMetric::
70 ComputeBalancedWeights(const GridFunction &nodes,
71  const TargetConstructor &tc, Vector &weights) const
72 {
73  const int m_cnt = tmop_q_arr.Size();
74  Vector averages;
75  ComputeAvgMetrics(nodes, tc, averages);
76  weights.SetSize(m_cnt);
77 
78  // For [ combo_A_B_C = a m_A + b m_B + c m_C ] we would have:
79  // a = BC / (AB + AC + BC), b = AC / (AB + AC + BC), c = AB / (AB + AC + BC),
80  // where A = avg_m_A, B = avg_m_B, C = avg_m_C.
81  // Nested loop to avoid division, as some avg may be 0.
82  Vector products_no_m(m_cnt); products_no_m = 1.0;
83  for (int m_p = 0; m_p < m_cnt; m_p++)
84  {
85  for (int m_a = 0; m_a < m_cnt; m_a++)
86  {
87  if (m_p != m_a) { products_no_m(m_p) *= averages(m_a); }
88  }
89  }
90  const double pnm_sum = products_no_m.Sum();
91 
92  if (pnm_sum == 0.0) { weights = 1.0 / m_cnt; return; }
93  for (int m = 0; m < m_cnt; m++) { weights(m) = products_no_m(m) / pnm_sum; }
94 
95  MFEM_ASSERT(fabs(weights.Sum() - 1.0) < 1e-14,
96  "Error: sum should be 1 always: " << weights.Sum());
97 }
98 
99 void TMOP_Combo_QualityMetric::ComputeAvgMetrics(const GridFunction &nodes,
100  const TargetConstructor &tc,
101  Vector &averages) const
102 {
103  const int m_cnt = tmop_q_arr.Size(),
104  NE = nodes.FESpace()->GetNE(),
105  dim = nodes.FESpace()->GetMesh()->Dimension();
106 
107  averages.SetSize(m_cnt);
108 
109  // Integrals of all metrics.
110  Array<int> pos_dofs;
111  averages = 0.0;
112  double volume = 0.0;
113  for (int e = 0; e < NE; e++)
114  {
115  const FiniteElement &fe_pos = *nodes.FESpace()->GetFE(e);
116  const IntegrationRule &ir = IntRules.Get(fe_pos.GetGeomType(),
117  2 * fe_pos.GetOrder());
118  const int nsp = ir.GetNPoints(), dof = fe_pos.GetDof();
119 
120  DenseMatrix dshape(dof, dim);
121  DenseMatrix pos(dof, dim);
122  pos.SetSize(dof, dim);
123  Vector posV(pos.Data(), dof * dim);
124 
125  nodes.FESpace()->GetElementVDofs(e, pos_dofs);
126  nodes.GetSubVector(pos_dofs, posV);
127 
128  DenseTensor W(dim, dim, nsp);
129  DenseMatrix Winv(dim), T(dim), A(dim);
130  tc.ComputeElementTargets(e, fe_pos, ir, posV, W);
131 
132  for (int q = 0; q < nsp; q++)
133  {
134  const DenseMatrix &Wj = W(q);
135  CalcInverse(Wj, Winv);
136 
137  const IntegrationPoint &ip = ir.IntPoint(q);
138  fe_pos.CalcDShape(ip, dshape);
139  MultAtB(pos, dshape, A);
140  Mult(A, Winv, T);
141 
142  const double w_detA = ip.weight * A.Det();
143  for (int m = 0; m < m_cnt; m++)
144  {
145  tmop_q_arr[m]->SetTargetJacobian(Wj);
146  averages(m) += tmop_q_arr[m]->EvalW(T) * w_detA;
147  }
148  volume += w_detA;
149  }
150  }
151 
152  // Parallel case.
153 #ifdef MFEM_USE_MPI
154  auto par_nodes = dynamic_cast<const ParGridFunction *>(&nodes);
155  if (par_nodes)
156  {
157  MPI_Allreduce(MPI_IN_PLACE, averages.GetData(), m_cnt,
158  MPI_DOUBLE, MPI_SUM, par_nodes->ParFESpace()->GetComm());
159  MPI_Allreduce(MPI_IN_PLACE, &volume, 1, MPI_DOUBLE, MPI_SUM,
160  par_nodes->ParFESpace()->GetComm());
161  }
162 #endif
163 
164  averages /= volume;
165 }
166 
167 double TMOP_WorstCaseUntangleOptimizer_Metric::EvalW(const DenseMatrix &Jpt)
168 const
169 {
170  double metric_tilde = EvalWBarrier(Jpt);
171  double metric = metric_tilde;
172  if (wctype == WorstCaseType::PMean)
173  {
174  metric = std::pow(metric_tilde, exponent);
175  }
176  else if (wctype == WorstCaseType::Beta)
177  {
178  double beta = max_muT+muT_ep;
179  metric = metric_tilde/(beta-metric_tilde);
180  }
181  return metric;
182 }
183 
184 double TMOP_WorstCaseUntangleOptimizer_Metric::EvalWBarrier(
185  const DenseMatrix &Jpt) const
186 {
187  double denominator = 1.0;
188  if (btype == BarrierType::Shifted)
189  {
190  denominator = 2.0*(Jpt.Det()-std::min(alpha*min_detT-detT_ep, 0.0));
191  }
192  else if (btype == BarrierType::Pseudo)
193  {
194  double detT = Jpt.Det();
195  denominator = detT + std::sqrt(detT*detT + detT_ep*detT_ep);
196  }
197  return tmop_metric.EvalW(Jpt)/denominator;
198 }
199 
200 double TMOP_Metric_001::EvalW(const DenseMatrix &Jpt) const
201 {
202  ie.SetJacobian(Jpt.GetData());
203  return ie.Get_I1();
204 }
205 
206 void TMOP_Metric_001::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
207 {
208  ie.SetJacobian(Jpt.GetData());
209  P = ie.Get_dI1();
210 }
211 
212 void TMOP_Metric_001::AssembleH(const DenseMatrix &Jpt,
213  const DenseMatrix &DS,
214  const double weight,
215  DenseMatrix &A) const
216 {
217  ie.SetJacobian(Jpt.GetData());
218  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
219  ie.Assemble_ddI1(weight, A.GetData());
220 }
221 
222 double TMOP_Metric_skew2D::EvalW(const DenseMatrix &Jpt) const
223 {
224  MFEM_VERIFY(Jtr != NULL,
225  "Requires a target Jacobian, use SetTargetJacobian().");
226 
227  DenseMatrix Jpr(2, 2);
228  Mult(Jpt, *Jtr, Jpr);
229 
230  Vector col1, col2;
231  Jpr.GetColumn(0, col1);
232  Jpr.GetColumn(1, col2);
233  double norm_prod = col1.Norml2() * col2.Norml2();
234  const double cos_Jpr = (col1 * col2) / norm_prod,
235  sin_Jpr = fabs(Jpr.Det()) / norm_prod;
236 
237  Jtr->GetColumn(0, col1);
238  Jtr->GetColumn(1, col2);
239  norm_prod = col1.Norml2() * col2.Norml2();
240  const double cos_Jtr = (col1 * col2) / norm_prod,
241  sin_Jtr = fabs(Jtr->Det()) / norm_prod;
242 
243  return 0.5 * (1.0 - cos_Jpr * cos_Jtr - sin_Jpr * sin_Jtr);
244 }
245 
246 double TMOP_Metric_skew3D::EvalW(const DenseMatrix &Jpt) const
247 {
248  MFEM_VERIFY(Jtr != NULL,
249  "Requires a target Jacobian, use SetTargetJacobian().");
250 
251  DenseMatrix Jpr(3, 3);
252  Mult(Jpt, *Jtr, Jpr);
253 
254  Vector col1, col2, col3;
255  Jpr.GetColumn(0, col1);
256  Jpr.GetColumn(1, col2);
257  Jpr.GetColumn(2, col3);
258  double norm_c1 = col1.Norml2(),
259  norm_c2 = col2.Norml2(),
260  norm_c3 = col3.Norml2();
261  double cos_Jpr_12 = (col1 * col2) / (norm_c1 * norm_c2),
262  cos_Jpr_13 = (col1 * col3) / (norm_c1 * norm_c3),
263  cos_Jpr_23 = (col2 * col3) / (norm_c2 * norm_c3);
264  double sin_Jpr_12 = std::sqrt(1.0 - cos_Jpr_12 * cos_Jpr_12),
265  sin_Jpr_13 = std::sqrt(1.0 - cos_Jpr_13 * cos_Jpr_13),
266  sin_Jpr_23 = std::sqrt(1.0 - cos_Jpr_23 * cos_Jpr_23);
267 
268  Jtr->GetColumn(0, col1);
269  Jtr->GetColumn(1, col2);
270  Jtr->GetColumn(2, col3);
271  norm_c1 = col1.Norml2();
272  norm_c2 = col2.Norml2(),
273  norm_c3 = col3.Norml2();
274  double cos_Jtr_12 = (col1 * col2) / (norm_c1 * norm_c2),
275  cos_Jtr_13 = (col1 * col3) / (norm_c1 * norm_c3),
276  cos_Jtr_23 = (col2 * col3) / (norm_c2 * norm_c3);
277  double sin_Jtr_12 = std::sqrt(1.0 - cos_Jtr_12 * cos_Jtr_12),
278  sin_Jtr_13 = std::sqrt(1.0 - cos_Jtr_13 * cos_Jtr_13),
279  sin_Jtr_23 = std::sqrt(1.0 - cos_Jtr_23 * cos_Jtr_23);
280 
281  return (3.0 - cos_Jpr_12 * cos_Jtr_12 - sin_Jpr_12 * sin_Jtr_12
282  - cos_Jpr_13 * cos_Jtr_13 - sin_Jpr_13 * sin_Jtr_13
283  - cos_Jpr_23 * cos_Jtr_23 - sin_Jpr_23 * sin_Jtr_23) / 6.0;
284 }
285 
286 double TMOP_Metric_aspratio2D::EvalW(const DenseMatrix &Jpt) const
287 {
288  MFEM_VERIFY(Jtr != NULL,
289  "Requires a target Jacobian, use SetTargetJacobian().");
290 
291  DenseMatrix Jpr(2, 2);
292  Mult(Jpt, *Jtr, Jpr);
293 
294  Vector col1, col2;
295  Jpr.GetColumn(0, col1);
296  Jpr.GetColumn(1, col2);
297  const double ratio_Jpr = col2.Norml2() / col1.Norml2();
298 
299  Jtr->GetColumn(0, col1);
300  Jtr->GetColumn(1, col2);
301  const double ratio_Jtr = col2.Norml2() / col1.Norml2();
302 
303  return 0.5 * (ratio_Jpr / ratio_Jtr + ratio_Jtr / ratio_Jpr) - 1.0;
304 }
305 
306 double TMOP_Metric_aspratio3D::EvalW(const DenseMatrix &Jpt) const
307 {
308  MFEM_VERIFY(Jtr != NULL,
309  "Requires a target Jacobian, use SetTargetJacobian().");
310 
311  DenseMatrix Jpr(3, 3);
312  Mult(Jpt, *Jtr, Jpr);
313 
314  Vector col1, col2, col3;
315  Jpr.GetColumn(0, col1);
316  Jpr.GetColumn(1, col2);
317  Jpr.GetColumn(2, col3);
318  double norm_c1 = col1.Norml2(),
319  norm_c2 = col2.Norml2(),
320  norm_c3 = col3.Norml2();
321  double ratio_Jpr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
322  ratio_Jpr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
323  ratio_Jpr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
324 
325  Jtr->GetColumn(0, col1);
326  Jtr->GetColumn(1, col2);
327  Jtr->GetColumn(2, col3);
328  norm_c1 = col1.Norml2();
329  norm_c2 = col2.Norml2();
330  norm_c3 = col3.Norml2();
331  double ratio_Jtr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
332  ratio_Jtr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
333  ratio_Jtr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
334 
335  return ( 0.5 * (ratio_Jpr_1 / ratio_Jtr_1 + ratio_Jtr_1 / ratio_Jpr_1) +
336  0.5 * (ratio_Jpr_2 / ratio_Jtr_2 + ratio_Jtr_2 / ratio_Jpr_2) +
337  0.5 * (ratio_Jpr_3 / ratio_Jtr_3 + ratio_Jtr_3 / ratio_Jpr_3) - 3.0
338  ) / 3.0;
339 }
340 
341 double TMOP_Metric_002::EvalWMatrixForm(const DenseMatrix &Jpt) const
342 {
343  return 0.5 * Jpt.FNorm2() / Jpt.Det() - 1.0;
344 }
345 
346 double TMOP_Metric_002::EvalW(const DenseMatrix &Jpt) const
347 {
348  ie.SetJacobian(Jpt.GetData());
349  return 0.5 * ie.Get_I1b() - 1.0;
350 }
351 
352 void TMOP_Metric_002::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
353 {
354  ie.SetJacobian(Jpt.GetData());
355  P.Set(0.5, ie.Get_dI1b());
356 }
357 
358 void TMOP_Metric_002::AssembleH(const DenseMatrix &Jpt,
359  const DenseMatrix &DS,
360  const double weight,
361  DenseMatrix &A) const
362 {
363  ie.SetJacobian(Jpt.GetData());
364  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
365  ie.Assemble_ddI1b(0.5*weight, A.GetData());
366 }
367 
368 double TMOP_Metric_004::EvalW(const DenseMatrix &Jpt) const
369 {
370  ie.SetJacobian(Jpt.GetData());
371  return ie.Get_I1() - 2.0*ie.Get_I2b();
372 }
373 
374 void TMOP_Metric_004::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
375 {
376  ie.SetJacobian(Jpt.GetData());
377  Add(1.0, ie.Get_dI1(), -2.0, ie.Get_dI2b(), P);
378 }
379 
380 void TMOP_Metric_004::AssembleH(const DenseMatrix &Jpt,
381  const DenseMatrix &DS,
382  const double weight,
383  DenseMatrix &A) const
384 {
385  ie.SetJacobian(Jpt.GetData());
386  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
387 
388  ie.Assemble_ddI1(weight, A.GetData());
389  ie.Assemble_ddI2b(-2.0*weight, A.GetData());
390 }
391 
392 double TMOP_Metric_007::EvalW(const DenseMatrix &Jpt) const
393 {
394  // mu_7 = |J-J^{-t}|^2 = |J|^2 + |J^{-1}|^2 - 4
395  ie.SetJacobian(Jpt.GetData());
396  return ie.Get_I1()*(1. + 1./ie.Get_I2()) - 4.0;
397 }
398 
399 void TMOP_Metric_007::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
400 {
401  // P = d(I1*(1 + 1/I2)) = (1 + 1/I2) dI1 - I1/I2^2 dI2
402  ie.SetJacobian(Jpt.GetData());
403  const double I2 = ie.Get_I2();
404  Add(1. + 1./I2, ie.Get_dI1(), -ie.Get_I1()/(I2*I2), ie.Get_dI2(), P);
405 }
406 
407 void TMOP_Metric_007::AssembleH(const DenseMatrix &Jpt,
408  const DenseMatrix &DS,
409  const double weight,
410  DenseMatrix &A) const
411 {
412  // P = d(I1*(1 + 1/I2))
413  // = (1 + 1/I2) dI1 - I1/I2^2 dI2
414  //
415  // dP = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 -
416  // (dI2 x d(I1/I2^2)) - I1/I2^2 ddI2
417  // = (-1/I2^2) (dI1 x dI2) + (1 + 1/I2) ddI1 +
418  // (-1/I2^2) (dI2 x [dI1 - 2 I1/I2 dI2]) - I1/I2^2 ddI2
419  // = (-1/I2^2) (dI1 x dI2 + dI2 x dI1) + (1 + 1/I2) ddI1 +
420  // (2 I1/I2^3) (dI2 x dI2) - I1/I2^2 ddI2
421  ie.SetJacobian(Jpt.GetData());
422  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
423  const double c1 = 1./ie.Get_I2();
424  const double c2 = weight*c1*c1;
425  const double c3 = ie.Get_I1()*c2;
426  ie.Assemble_ddI1(weight*(1. + c1), A.GetData());
427  ie.Assemble_ddI2(-c3, A.GetData());
428  ie.Assemble_TProd(-c2, ie.Get_dI1(), ie.Get_dI2(), A.GetData());
429  ie.Assemble_TProd(2*c1*c3, ie.Get_dI2(), A.GetData());
430 }
431 
432 double TMOP_Metric_009::EvalW(const DenseMatrix &Jpt) const
433 {
434  // mu_9 = det(J)*|J-J^{-t}|^2 = I1b * (I2b^2 + 1) - 4 * I2b
435  // = (I1 - 4)*I2b + I1b
436  ie.SetJacobian(Jpt.GetData());
437  return (ie.Get_I1() - 4.0)*ie.Get_I2b() + ie.Get_I1b();
438 }
439 
440 void TMOP_Metric_009::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
441 {
442  // mu_9 = (I1 - 4)*I2b + I1b
443  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
444  ie.SetJacobian(Jpt.GetData());
445  Add(ie.Get_I1() - 4.0, ie.Get_dI2b(), ie.Get_I2b(), ie.Get_dI1(), P);
446  P += ie.Get_dI1b();
447 }
448 
449 void TMOP_Metric_009::AssembleH(const DenseMatrix &Jpt,
450  const DenseMatrix &DS,
451  const double weight,
452  DenseMatrix &A) const
453 {
454  // P = (I1 - 4)*dI2b + I2b*dI1 + dI1b
455  // dP = dI2b x dI1 + (I1-4)*ddI2b + dI1 x dI2b + I2b*ddI1 + ddI1b
456  // = (dI1 x dI2b + dI2b x dI1) + (I1-4)*ddI2b + I2b*ddI1 + ddI1b
457  ie.SetJacobian(Jpt.GetData());
458  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
459  ie.Assemble_TProd(weight, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
460  ie.Assemble_ddI2b(weight*(ie.Get_I1()-4.0), A.GetData());
461  ie.Assemble_ddI1(weight*ie.Get_I2b(), A.GetData());
462  ie.Assemble_ddI1b(weight, A.GetData());
463 }
464 
465 // mu_14 = |T-I|^2
466 double TMOP_Metric_014::EvalW(const DenseMatrix &Jpt) const
467 {
468  MFEM_VERIFY(Jtr != NULL,
469  "Requires a target Jacobian, use SetTargetJacobian().");
470 
471  DenseMatrix Id(2,2);
472 
473  Id(0,0) = 1; Id(0,1) = 0;
474  Id(1,0) = 0; Id(1,1) = 1;
475 
476  DenseMatrix Mat(2,2);
477  Mat = Jpt;
478  Mat.Add(-1,Id);
479  return Mat.FNorm2();
480 }
481 
482 double TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const
483 {
484  // mu_22 = (0.5*|J|^2 - det(J)) / (det(J) - tau0)
485  // = (0.5*I1 - I2b) / (I2b - tau0)
486  ie.SetJacobian(Jpt.GetData());
487  const double I2b = ie.Get_I2b();
488 
489  double d = I2b - min_detT;
490  if (d < 0.0 && min_detT == 0.0)
491  {
492  // The mesh has been untangled, but it's still possible to get negative
493  // detJ in FD calculations, as they move the nodes around with some small
494  // increments and can produce negative determinants. Thus we put a small
495  // value in the denominator. Note that here I2b < 0.
496  d = - I2b * 0.1;
497  }
498 
499  return (0.5*ie.Get_I1() - I2b) / d;
500 }
501 
502 void TMOP_Metric_022::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
503 {
504  // mu_22 = (0.5*I1 - I2b) / (I2b - tau0)
505  // P = 1/(I2b - tau0)*(0.5*dI1 - dI2b) - (0.5*I1 - I2b)/(I2b - tau0)^2*dI2b
506  // = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
507  ie.SetJacobian(Jpt.GetData());
508  const double c1 = 1.0/(ie.Get_I2b() - min_detT);
509  Add(c1/2, ie.Get_dI1(), (min_detT - ie.Get_I1()/2)*c1*c1, ie.Get_dI2b(), P);
510 }
511 
512 void TMOP_Metric_022::AssembleH(const DenseMatrix &Jpt,
513  const DenseMatrix &DS,
514  const double weight,
515  DenseMatrix &A) const
516 {
517  // P = 0.5/(I2b - tau0)*dI1 + (tau0 - 0.5*I1)/(I2b - tau0)^2*dI2b
518  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b) + 0.5/(I2b - tau0)*ddI1
519  // + (dI2b x dz) + z*ddI2b
520  //
521  // z = (tau0 - 0.5*I1)/(I2b - tau0)^2
522  // dz = -0.5/(I2b - tau0)^2*dI1 - 2*(tau0 - 0.5*I1)/(I2b - tau0)^3*dI2b
523  //
524  // dP = -0.5/(I2b - tau0)^2*(dI1 x dI2b + dI2b x dI1)
525  // -2*z/(I2b - tau0)*(dI2b x dI2b)
526  // +0.5/(I2b - tau0)*ddI1 + z*ddI2b
527  ie.SetJacobian(Jpt.GetData());
528  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
529  const double c1 = 1.0/(ie.Get_I2b() - min_detT);
530  const double c2 = weight*c1/2;
531  const double c3 = c1*c2;
532  const double c4 = (2*min_detT - ie.Get_I1())*c3; // weight*z
533  ie.Assemble_TProd(-c3, ie.Get_dI1(), ie.Get_dI2b(), A.GetData());
534  ie.Assemble_TProd(-2*c1*c4, ie.Get_dI2b(), A.GetData());
535  ie.Assemble_ddI1(c2, A.GetData());
536  ie.Assemble_ddI2b(c4, A.GetData());
537 }
538 
539 double TMOP_Metric_050::EvalWMatrixForm(const DenseMatrix &Jpt) const
540 {
541  // mu_50 = 0.5 |J^t J|^2 / det(J)^2 - 1.
542  DenseMatrix JtJ(2);
543  MultAAt(Jpt, JtJ);
544  JtJ.Transpose();
545  double det = Jpt.Det();
546 
547  return 0.5 * JtJ.FNorm2()/(det*det) - 1.0;
548 }
549 
550 double TMOP_Metric_050::EvalW(const DenseMatrix &Jpt) const
551 {
552  // mu_50 = 0.5*|J^t J|^2/det(J)^2 - 1
553  // = 0.5*(l1^4 + l2^4)/(l1*l2)^2 - 1
554  // = 0.5*((l1/l2)^2 + (l2/l1)^2) - 1 = 0.5*(l1/l2 - l2/l1)^2
555  // = 0.5*(l1/l2 + l2/l1)^2 - 2 = 0.5*I1b^2 - 2.
556  ie.SetJacobian(Jpt.GetData());
557  const double I1b = ie.Get_I1b();
558  return 0.5*I1b*I1b - 2.0;
559 }
560 
561 void TMOP_Metric_050::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
562 {
563  // mu_50 = 0.5*I1b^2 - 2
564  // P = I1b*dI1b
565  ie.SetJacobian(Jpt.GetData());
566  P.Set(ie.Get_I1b(), ie.Get_dI1b());
567 }
568 
569 void TMOP_Metric_050::AssembleH(const DenseMatrix &Jpt,
570  const DenseMatrix &DS,
571  const double weight,
572  DenseMatrix &A) const
573 {
574  // P = I1b*dI1b
575  // dP = dI1b x dI1b + I1b*ddI1b
576  ie.SetJacobian(Jpt.GetData());
577  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
578  ie.Assemble_TProd(weight, ie.Get_dI1b(), A.GetData());
579  ie.Assemble_ddI1b(weight*ie.Get_I1b(), A.GetData());
580 }
581 
582 double TMOP_Metric_055::EvalW(const DenseMatrix &Jpt) const
583 {
584  // mu_55 = (det(J) - 1)^2 = (I2b - 1)^2
585  ie.SetJacobian(Jpt.GetData());
586  const double c1 = ie.Get_I2b() - 1.0;
587  return c1*c1;
588 }
589 
590 void TMOP_Metric_055::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
591 {
592  // mu_55 = (I2b - 1)^2
593  // P = 2*(I2b - 1)*dI2b
594  ie.SetJacobian(Jpt.GetData());
595  P.Set(2*(ie.Get_I2b() - 1.0), ie.Get_dI2b());
596 }
597 
598 void TMOP_Metric_055::AssembleH(const DenseMatrix &Jpt,
599  const DenseMatrix &DS,
600  const double weight,
601  DenseMatrix &A) const
602 {
603  // P = 2*(I2b - 1)*dI2b
604  // dP = 2*(dI2b x dI2b) + 2*(I2b - 1)*ddI2b
605  ie.SetJacobian(Jpt.GetData());
606  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
607  ie.Assemble_TProd(2*weight, ie.Get_dI2b(), A.GetData());
608  ie.Assemble_ddI2b(2*weight*(ie.Get_I2b() - 1.0), A.GetData());
609 }
610 
611 double TMOP_Metric_056::EvalWMatrixForm(const DenseMatrix &Jpt) const
612 {
613  // mu_56 = 0.5 (det(J) + 1 / det(J)) - 1.
614  const double d = Jpt.Det();
615  return 0.5 * (d + 1.0 / d) - 1.0;
616 }
617 
618 double TMOP_Metric_056::EvalW(const DenseMatrix &Jpt) const
619 {
620  // mu_56 = 0.5*(I2b + 1/I2b) - 1.
621  ie.SetJacobian(Jpt.GetData());
622  const double I2b = ie.Get_I2b();
623  return 0.5*(I2b + 1.0/I2b) - 1.0;
624 }
625 
626 void TMOP_Metric_056::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
627 {
628  // mu_56 = 0.5*(I2b + 1/I2b) - 1.
629  // P = 0.5*(1 - 1/I2b^2)*dI2b.
630  ie.SetJacobian(Jpt.GetData());
631  P.Set(0.5 - 0.5/ie.Get_I2(), ie.Get_dI2b());
632 }
633 
634 void TMOP_Metric_056::AssembleH(const DenseMatrix &Jpt,
635  const DenseMatrix &DS,
636  const double weight,
637  DenseMatrix &A) const
638 {
639  // P = 0.5*(1 - 1/I2b^2)*dI2b.
640  // dP = (1/I2b^3)*(dI2b x dI2b) + (0.5 - 0.5/I2)*ddI2b.
641  ie.SetJacobian(Jpt.GetData());
642  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
643  ie.Assemble_TProd(weight/(ie.Get_I2()*ie.Get_I2b()),
644  ie.Get_dI2b(), A.GetData());
645  ie.Assemble_ddI2b(weight*(0.5 - 0.5/ie.Get_I2()), A.GetData());
646 }
647 
648 double TMOP_Metric_058::EvalWMatrixForm(const DenseMatrix &Jpt) const
649 {
650  // mu_58 = |J^t J|^2 / det(J)^2 - 2|J|^2 / det(J) + 2.
651  DenseMatrix JtJ(2);
652  MultAAt(Jpt, JtJ);
653  JtJ.Transpose();
654  double det = Jpt.Det();
655 
656  return JtJ.FNorm2()/(det*det) - 2*Jpt.FNorm2()/det + 2.0;
657 }
658 
659 double TMOP_Metric_058::EvalW(const DenseMatrix &Jpt) const
660 {
661  // mu_58 = I1b*(I1b - 2)
662  ie.SetJacobian(Jpt.GetData());
663  const double I1b = ie.Get_I1b();
664  return I1b*(I1b - 2.0);
665 }
666 
667 void TMOP_Metric_058::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
668 {
669  // mu_58 = I1b*(I1b - 2)
670  // P = (2*I1b - 2)*dI1b
671  ie.SetJacobian(Jpt.GetData());
672  P.Set(2*ie.Get_I1b() - 2.0, ie.Get_dI1b());
673 }
674 
675 void TMOP_Metric_058::AssembleH(const DenseMatrix &Jpt,
676  const DenseMatrix &DS,
677  const double weight,
678  DenseMatrix &A) const
679 {
680  // P = (2*I1b - 2)*dI1b
681  // dP = 2*(dI1b x dI1b) + (2*I1b - 2)*ddI1b
682  ie.SetJacobian(Jpt.GetData());
683  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
684  ie.Assemble_TProd(2*weight, ie.Get_dI1b(), A.GetData());
685  ie.Assemble_ddI1b(weight*(2*ie.Get_I1b() - 2.0), A.GetData());
686 }
687 
688 double TMOP_Metric_077::EvalWMatrixForm(const DenseMatrix &Jpt) const
689 {
690  // mu_77 = 0.5 (det(J)^2 + 1 / det(J)^2) - 1.
691  const double d = Jpt.Det();
692  return 0.5 * (d*d + 1.0/(d*d)) - 1.0;
693 }
694 double TMOP_Metric_077::EvalW(const DenseMatrix &Jpt) const
695 {
696  // mu_77 = 0.5 (I2 + 1 / I2) - 1.0.
697  ie.SetJacobian(Jpt.GetData());
698  const double I2 = ie.Get_I2();
699  return 0.5*(I2 + 1.0/I2) - 1.0;
700 }
701 
702 void TMOP_Metric_077::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
703 {
704  // mu_77 = 0.5 (I2 + 1 / I2) - 1.0.
705  // P = 1/2 (1 - 1/I2^2) dI2_dJ.
706  ie.SetJacobian(Jpt.GetData());
707  const double I2 = ie.Get_I2();
708  P.Set(0.5 * (1.0 - 1.0 / (I2 * I2)), ie.Get_dI2());
709 }
710 
711 void TMOP_Metric_077::AssembleH(const DenseMatrix &Jpt,
712  const DenseMatrix &DS,
713  const double weight,
714  DenseMatrix &A) const
715 {
716  ie.SetJacobian(Jpt.GetData());
717  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
718  const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
719  ie.Assemble_ddI2(weight*0.5*(1.0 - I2inv_sq), A.GetData());
720  ie.Assemble_TProd(weight * I2inv_sq / I2, ie.Get_dI2(), A.GetData());
721 }
722 
723 // mu_85 = |T-T'|^2, where T'= |T|*I/sqrt(2)
724 double TMOP_Metric_085::EvalW(const DenseMatrix &Jpt) const
725 {
726  MFEM_VERIFY(Jtr != NULL,
727  "Requires a target Jacobian, use SetTargetJacobian().");
728 
729  DenseMatrix Id(2,2);
730  DenseMatrix Mat(2,2);
731  Mat = Jpt;
732 
733  Id(0,0) = 1; Id(0,1) = 0;
734  Id(1,0) = 0; Id(1,1) = 1;
735  Id *= Mat.FNorm()/pow(2,0.5);
736 
737  Mat.Add(-1.,Id);
738  return Mat.FNorm2();
739 }
740 
741 // mu_98 = 1/(tau)|T-I|^2
742 double TMOP_Metric_098::EvalW(const DenseMatrix &Jpt) const
743 {
744  MFEM_VERIFY(Jtr != NULL,
745  "Requires a target Jacobian, use SetTargetJacobian().");
746 
747  DenseMatrix Id(2,2);
748 
749  Id(0,0) = 1; Id(0,1) = 0;
750  Id(1,0) = 0; Id(1,1) = 1;
751 
752  DenseMatrix Mat(2,2);
753  Mat = Jpt;
754  Mat.Add(-1,Id);
755  return Mat.FNorm2()/Jtr->Det();
756 }
757 
758 double TMOP_Metric_211::EvalW(const DenseMatrix &Jpt) const
759 {
760  // mu_211 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
761  // = (I2b - 1)^2 - I2b + sqrt(I2b^2 + eps)
762  ie.SetJacobian(Jpt.GetData());
763  const double I2b = ie.Get_I2b();
764  return (I2b - 1.0)*(I2b - 1.0) - I2b + std::sqrt(I2b*I2b + eps);
765 }
766 
767 void TMOP_Metric_211::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
768 {
769  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
770 }
771 
772 void TMOP_Metric_211::AssembleH(const DenseMatrix &Jpt,
773  const DenseMatrix &DS,
774  const double weight,
775  DenseMatrix &A) const
776 {
777  MFEM_ABORT("Metric not implemented yet. Use metric mu_55 instead.");
778 }
779 
780 double TMOP_Metric_252::EvalW(const DenseMatrix &Jpt) const
781 {
782  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0).
783  ie.SetJacobian(Jpt.GetData());
784  const double I2b = ie.Get_I2b();
785  return 0.5*(I2b - 1.0)*(I2b - 1.0)/(I2b - tau0);
786 }
787 
788 void TMOP_Metric_252::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
789 {
790  // mu_252 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
791  // P = (c - 0.5*c*c ) * dI2b
792  //
793  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
794  ie.SetJacobian(Jpt.GetData());
795  const double I2b = ie.Get_I2b();
796  const double c = (I2b - 1.0)/(I2b - tau0);
797  P.Set(c - 0.5*c*c, ie.Get_dI2b());
798 }
799 
800 void TMOP_Metric_252::AssembleH(const DenseMatrix &Jpt,
801  const DenseMatrix &DS,
802  const double weight,
803  DenseMatrix &A) const
804 {
805  // c = (I2b - 1)/(I2b - tau0), see TMOP_Metric_352 for details
806  //
807  // P = (c - 0.5*c*c ) * dI2b
808  // dP = (1 - c)^2/(I2b - tau0)*(dI2b x dI2b) + (c - 0.5*c*c)*ddI2b
809  ie.SetJacobian(Jpt.GetData());
810  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
811  const double I2b = ie.Get_I2b();
812  const double c0 = 1.0/(I2b - tau0);
813  const double c = c0*(I2b - 1.0);
814  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI2b(), A.GetData());
815  ie.Assemble_ddI2b(weight*(c - 0.5*c*c), A.GetData());
816 }
817 
818 double TMOP_Metric_301::EvalWMatrixForm(const DenseMatrix &Jpt) const
819 {
820  // mu_301 = 1/3 |J| |J^-1| - 1.
821  ie.SetJacobian(Jpt.GetData());
822  DenseMatrix inv(3);
823  CalcInverse(Jpt, inv);
824  return Jpt.FNorm() * inv.FNorm() / 3.0 - 1.0;
825 }
826 
827 double TMOP_Metric_301::EvalW(const DenseMatrix &Jpt) const
828 {
829  // mu_301 = 1/3 sqrt(I1b * I2b) - 1
830  ie.SetJacobian(Jpt.GetData());
831  return std::sqrt(ie.Get_I1b()*ie.Get_I2b())/3. - 1.;
832 }
833 
834 void TMOP_Metric_301::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
835 {
836  // W = (1/3)*sqrt(I1b*I2b) - 1
837  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
838  ie.SetJacobian(Jpt.GetData());
839  const double a = 1./(6.*std::sqrt(ie.Get_I1b()*ie.Get_I2b()));
840  Add(a*ie.Get_I2b(), ie.Get_dI1b(), a*ie.Get_I1b(), ie.Get_dI2b(), P);
841 }
842 
843 void TMOP_Metric_301::AssembleH(const DenseMatrix &Jpt,
844  const DenseMatrix &DS,
845  const double weight,
846  DenseMatrix &A) const
847 {
848  // dW = (1/6)/sqrt(I1b*I2b)*[I2b*dI1b + I1b*dI2b]
849  // dW = (1/6)*[z2*dI1b + z1*dI2b], z1 = sqrt(I1b/I2b), z2 = sqrt(I2b/I1b)
850  // ddW = (1/6)*[dI1b x dz2 + z2*ddI1b + dI2b x dz1 + z1*ddI2b]
851  //
852  // dz1 = (1/2)*sqrt(I2b/I1b) [ (1/I2b)*dI1b - (I1b/(I2b*I2b))*dI2b ]
853  // = (1/2)/sqrt(I1b*I2b) [ dI1b - (I1b/I2b)*dI2b ]
854  // dz2 = (1/2)/sqrt(I1b*I2b) [ dI2b - (I2b/I1b)*dI1b ]
855  //
856  // dI1b x dz2 + dI2b x dz1 =
857  // (1/2)/sqrt(I1b*I2b) dI1b x [ dI2b - (I2b/I1b)*dI1b ] +
858  // (1/2)/sqrt(I1b*I2b) dI2b x [ dI1b - (I1b/I2b)*dI2b ] =
859  // (1/2)/sqrt(I1b*I2b) [sqrt(I1b/I2b)*dI2b - sqrt(I2b/I1b)*dI1b] x
860  // [sqrt(I2b/I1b)*dI1b - sqrt(I1b/I2b)*dI2b] =
861  // (1/2)*(I1b*I2b)^{-3/2} (I1b*dI2b - I2b*dI1b) x (I2b*dI1b - I1b*dI2b)
862  // and the last two parentheses are the same up to a sign.
863  //
864  // z1 = I1b/sqrt(I1b*I2b), z2 = I2b/sqrt(I1b*I2b)
865 
866  ie.SetJacobian(Jpt.GetData());
867  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
868  double X_data[9];
869  DenseMatrix X(X_data, 3, 3);
870  Add(- ie.Get_I2b(), ie.Get_dI1b(), ie.Get_I1b(), ie.Get_dI2b(), X);
871  const double I1b_I2b = ie.Get_I1b()*ie.Get_I2b();
872  const double a = weight/(6*std::sqrt(I1b_I2b));
873  ie.Assemble_ddI1b(a*ie.Get_I2b(), A.GetData());
874  ie.Assemble_ddI2b(a*ie.Get_I1b(), A.GetData());
875  ie.Assemble_TProd(-a/(2*I1b_I2b), X_data, A.GetData());
876 }
877 
878 double TMOP_Metric_302::EvalWMatrixForm(const DenseMatrix &Jpt) const
879 {
880  // mu_301 = |J|^2 |J^{-1}|^2 / 9 - 1.
881  ie.SetJacobian(Jpt.GetData());
882  DenseMatrix inv(3);
883  CalcInverse(Jpt, inv);
884  return Jpt.FNorm2() * inv.FNorm2() / 9.0 - 1.0;
885 }
886 
887 double TMOP_Metric_302::EvalW(const DenseMatrix &Jpt) const
888 {
889  // mu_2 = |J|^2 |J^{-1}|^2 / 9 - 1
890  // = (l1^2 + l2^2 + l3^3)*(l1^{-2} + l2^{-2} + l3^{-2}) / 9 - 1
891  // = I1*(l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/l1^2/l2^2/l3^2/9 - 1
892  // = I1*I2/det(J)^2/9 - 1 = I1b*I2b/9-1
893  ie.SetJacobian(Jpt.GetData());
894  return ie.Get_I1b()*ie.Get_I2b()/9. - 1.;
895 }
896 
897 void TMOP_Metric_302::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
898 {
899  // mu_2 = I1b*I2b/9-1
900  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
901  ie.SetJacobian(Jpt.GetData());
902  Add(ie.Get_I1b()/9, ie.Get_dI2b(), ie.Get_I2b()/9, ie.Get_dI1b(), P);
903 }
904 
905 void TMOP_Metric_302::AssembleH(const DenseMatrix &Jpt,
906  const DenseMatrix &DS,
907  const double weight,
908  DenseMatrix &A) const
909 {
910  // P = (I1b/9)*dI2b + (I2b/9)*dI1b
911  // dP = (dI2b x dI1b)/9 + (I1b/9)*ddI2b + (dI1b x dI2b)/9 + (I2b/9)*ddI1b
912  // = (dI2b x dI1b + dI1b x dI2b)/9 + (I1b/9)*ddI2b + (I2b/9)*ddI1b
913  ie.SetJacobian(Jpt.GetData());
914  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
915  const double c1 = weight/9;
916  ie.Assemble_TProd(c1, ie.Get_dI1b(), ie.Get_dI2b(), A.GetData());
917  ie.Assemble_ddI2b(c1*ie.Get_I1b(), A.GetData());
918  ie.Assemble_ddI1b(c1*ie.Get_I2b(), A.GetData());
919 }
920 
921 double TMOP_Metric_303::EvalWMatrixForm(const DenseMatrix &Jpt) const
922 {
923  // mu_303 = |J|^2 / 3 / det(J)^(2/3) - 1.
924  ie.SetJacobian(Jpt.GetData());
925  return Jpt.FNorm2() / 3.0 / pow(Jpt.Det(), 2.0 / 3.0) - 1.0;
926 }
927 
928 double TMOP_Metric_303::EvalW(const DenseMatrix &Jpt) const
929 {
930  // mu_303 = |J|^2 / 3 / det(J)^(2/3) - 1 = I1b/3 - 1.
931  ie.SetJacobian(Jpt.GetData());
932  return ie.Get_I1b()/3.0 - 1.0;
933 }
934 
935 void TMOP_Metric_303::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
936 {
937  // mu_304 = I1b/3 - 1.
938  // P = dI1b/3.
939  ie.SetJacobian(Jpt.GetData());
940  P.Set(1./3., ie.Get_dI1b());
941 }
942 
943 void TMOP_Metric_303::AssembleH(const DenseMatrix &Jpt,
944  const DenseMatrix &DS,
945  const double weight,
946  DenseMatrix &A) const
947 {
948  // P = dI1b/3.
949  // dP = ddI1b/3.
950  ie.SetJacobian(Jpt.GetData());
951  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
952  ie.Assemble_ddI1b(weight/3., A.GetData());
953 }
954 
955 double TMOP_Metric_304::EvalWMatrixForm(const DenseMatrix &Jpt) const
956 {
957  // mu_304 = |J|^3 / 3^(3/2) / det(J) - 1
958  const double fnorm = Jpt.FNorm();
959  return fnorm * fnorm * fnorm / pow(3.0, 1.5) / Jpt.Det() - 1.0;
960 }
961 
962 double TMOP_Metric_304::EvalW(const DenseMatrix &Jpt) const
963 {
964  // mu_304 = (I1b/3)^3/2 - 1.
965  ie.SetJacobian(Jpt.GetData());
966  return pow(ie.Get_I1b()/3.0, 1.5) - 1.0;
967 }
968 
969 void TMOP_Metric_304::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
970 {
971  // mu_304 = (I1b/3)^3/2 - 1.
972  // P = 3/2 * (I1b/3)^1/2 * dI1b / 3 = 1/2 * (I1b/3)^1/2 * dI1b.
973  ie.SetJacobian(Jpt.GetData());
974  P.Set(0.5 * sqrt(ie.Get_I1b()/3.0), ie.Get_dI1b());
975 }
976 
977 void TMOP_Metric_304::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
978  const double weight, DenseMatrix &A) const
979 {
980  // P = 1/2 * (I1b/3)^1/2 * dI1b.
981  // dP = 1/12 * (I1b/3)^(-1/2) * (dI1b x dI1b) + 1/2 * (I1b/3)^1/2 * ddI1b.
982  ie.SetJacobian(Jpt.GetData());
983  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
984  ie.Assemble_TProd(weight / 12.0 / sqrt(ie.Get_I1b()/3.0),
985  ie.Get_dI1b(), A.GetData());
986  ie.Assemble_ddI1b(weight / 2.0 * sqrt(ie.Get_I1b()/3.0), A.GetData());
987 }
988 
989 double TMOP_Metric_311::EvalW(const DenseMatrix &Jpt) const
990 {
991  // mu_311 = (det(J) - 1)^2 - det(J) + (det(J)^2 + eps)^{1/2}
992  // = (I3b - 1)^2 - I3b + sqrt(I3b^2 + eps)
993  ie.SetJacobian(Jpt.GetData());
994  const double I3b = ie.Get_I3b();
995  return (I3b - 1.0)*(I3b - 1.0) - I3b + std::sqrt(I3b*I3b + eps);
996 }
997 
998 void TMOP_Metric_311::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
999 {
1000  ie.SetJacobian(Jpt.GetData());
1001  const double I3b = ie.Get_I3b();
1002  const double c = 2*I3b-3+(I3b)/(std::pow((I3b*I3b+eps),0.5));
1003  P.Set(c, ie.Get_dI3b());
1004 }
1005 
1006 void TMOP_Metric_311::AssembleH(const DenseMatrix &Jpt,
1007  const DenseMatrix &DS,
1008  const double weight,
1009  DenseMatrix &A) const
1010 {
1011  ie.SetJacobian(Jpt.GetData());
1012  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1013  const double I3b = ie.Get_I3b();
1014  const double c0 = I3b*I3b+eps;
1015  const double c1 = 2 + 1/(pow(c0,0.5)) - I3b*I3b/(pow(c0,1.5));
1016  const double c2 = 2*I3b - 3 + I3b/(pow(c0,0.5));
1017  ie.Assemble_TProd(weight*c1, ie.Get_dI3b(), A.GetData());
1018  ie.Assemble_ddI3b(c2*weight, A.GetData());
1019 }
1020 
1021 double TMOP_Metric_313::EvalW(const DenseMatrix &Jpt) const
1022 {
1023  ie.SetJacobian(Jpt.GetData());
1024 
1025  const double I3b = ie.Get_I3b();
1026  double d = I3b - min_detT;
1027  if (d < 0.0 && min_detT == 0.0)
1028  {
1029  // The mesh has been untangled, but it's still possible to get negative
1030  // detJ in FD calculations, as they move the nodes around with some small
1031  // increments and can produce negative determinants. Thus we put a small
1032  // value in the denominator. Note that here I3b < 0.
1033  d = - I3b * 0.1;
1034  }
1035 
1036  const double c = std::pow(d, -2.0/3.0);
1037 
1038  return ie.Get_I1() * c / 3.0;
1039 }
1040 
1041 void TMOP_Metric_313::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1042 {
1043  MFEM_ABORT("Metric not implemented yet.");
1044 }
1045 
1046 void TMOP_Metric_313::AssembleH(const DenseMatrix &Jpt,
1047  const DenseMatrix &DS,
1048  const double weight,
1049  DenseMatrix &A) const
1050 {
1051  MFEM_ABORT("Metric not implemented yet.");
1052 }
1053 
1054 double TMOP_Metric_315::EvalW(const DenseMatrix &Jpt) const
1055 {
1056  // mu_315 = mu_15_3D = (det(J) - 1)^2
1057  ie.SetJacobian(Jpt.GetData());
1058  const double c1 = ie.Get_I3b() - 1.0;
1059  return c1*c1;
1060 }
1061 
1062 void TMOP_Metric_315::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1063 {
1064  // mu_315 = (I3b - 1)^2
1065  // P = 2*(I3b - 1)*dI3b
1066  ie.SetJacobian(Jpt.GetData());
1067  P.Set(2*(ie.Get_I3b() - 1.0), ie.Get_dI3b());
1068 }
1069 
1070 void TMOP_Metric_315::AssembleH(const DenseMatrix &Jpt,
1071  const DenseMatrix &DS,
1072  const double weight,
1073  DenseMatrix &A) const
1074 {
1075  // P = 2*(I3b - 1)*dI3b
1076  // dP = 2*(dI3b x dI3b) + 2*(I3b - 1)*ddI3b
1077  ie.SetJacobian(Jpt.GetData());
1078  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1079  ie.Assemble_TProd(2*weight, ie.Get_dI3b(), A.GetData());
1080  ie.Assemble_ddI3b(2*weight*(ie.Get_I3b() - 1.0), A.GetData());
1081 }
1082 
1083 double TMOP_Metric_316::EvalWMatrixForm(const DenseMatrix &Jpt) const
1084 {
1085  // mu_316 = 0.5 (det(J) + 1/det(J)) - 1.
1086  return 0.5 * (Jpt.Det() + 1.0 / Jpt.Det()) - 1.0;
1087 }
1088 
1089 double TMOP_Metric_316::EvalW(const DenseMatrix &Jpt) const
1090 {
1091  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
1092  ie.SetJacobian(Jpt.GetData());
1093  const double I3b = ie.Get_I3b();
1094  return 0.5*(I3b + 1.0/I3b) - 1.0;
1095 }
1096 
1097 void TMOP_Metric_316::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1098 {
1099  // mu_316 = mu_16_3D = 0.5*(I3b + 1/I3b) - 1
1100  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
1101  ie.SetJacobian(Jpt.GetData());
1102  P.Set(0.5 - 0.5/ie.Get_I3(), ie.Get_dI3b());
1103 }
1104 
1105 void TMOP_Metric_316::AssembleH(const DenseMatrix &Jpt,
1106  const DenseMatrix &DS,
1107  const double weight,
1108  DenseMatrix &A) const
1109 {
1110  // P = 0.5*(1 - 1/I3b^2)*dI3b = (0.5 - 0.5/I3)*dI3b
1111  // dP = (1/I3b^3)*(dI3b x dI3b) + (0.5 - 0.5/I3)*ddI3b
1112  ie.SetJacobian(Jpt.GetData());
1113  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1114  ie.Assemble_TProd(weight/(ie.Get_I3()*ie.Get_I3b()),
1115  ie.Get_dI3b(), A.GetData());
1116  ie.Assemble_ddI3b(weight*(0.5 - 0.5/ie.Get_I3()), A.GetData());
1117 }
1118 
1119 double TMOP_Metric_318::EvalWMatrixForm(const DenseMatrix &Jpt) const
1120 {
1121  // mu_318 = 0.5 (det(J)^2 + 1/det(J)^2) - 1.
1122  double d = Jpt.Det();
1123  return 0.5 * (d*d + 1.0 / (d*d)) - 1.0;
1124 }
1125 
1126 double TMOP_Metric_318::EvalW(const DenseMatrix &Jpt) const
1127 {
1128  // mu_318 = mu_77_3D = 0.5 * (I3 + 1/I3) - 1.
1129  ie.SetJacobian(Jpt.GetData());
1130  const double I3 = ie.Get_I3();
1131  return 0.5*(I3 + 1.0/I3) - 1.0;
1132 }
1133 
1134 void TMOP_Metric_318::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1135 {
1136  // mu_318 = mu_77_3D = 0.5*(I3 + 1/I3) - 1.
1137  // P = 0.5*(1 - 1/I3^2)*dI3 = (0.5 - 0.5/I3^2)*dI3.
1138  ie.SetJacobian(Jpt.GetData());
1139  P.Set(0.5 - 0.5/(ie.Get_I3()*ie.Get_I3()), ie.Get_dI3());
1140 }
1141 
1142 void TMOP_Metric_318::AssembleH(const DenseMatrix &Jpt,
1143  const DenseMatrix &DS,
1144  const double weight,
1145  DenseMatrix &A) const
1146 {
1147  // P = (0.5 - 0.5/I3^2)*dI3.
1148  // dP = (1/I3^3)*(dI3 x dI3) +(0.5 - 0.5/I3^2)*ddI3
1149  ie.SetJacobian(Jpt.GetData());
1150  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1151  const double i3 = ie.Get_I3();
1152  ie.Assemble_TProd(weight/(i3 * i3 * i3), ie.Get_dI3(), A.GetData());
1153  ie.Assemble_ddI3(weight*(0.5 - 0.5 / (i3 * i3)), A.GetData());
1154 }
1155 
1156 double TMOP_Metric_321::EvalWMatrixForm(const DenseMatrix &Jpt) const
1157 {
1158  // mu_321 = |J - J^-t|^2.
1159  ie.SetJacobian(Jpt.GetData());
1160  DenseMatrix invt(3);
1161  CalcInverseTranspose(Jpt, invt);
1162  invt.Add(-1.0, Jpt);
1163  return invt.FNorm2();
1164 }
1165 
1166 double TMOP_Metric_321::EvalW(const DenseMatrix &Jpt) const
1167 {
1168  // mu_321 = mu_21_3D = |J - J^{-t}|^2
1169  // = |J|^2 + |J^{-1}|^2 - 6
1170  // = |J|^2 + (l1^{-2} + l2^{-2} + l3^{-2}) - 6
1171  // = |J|^2 + (l2^2*l3^2 + l1^2*l3^2 + l1^2*l2^2)/det(J)^2 - 6
1172  // = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
1173  ie.SetJacobian(Jpt.GetData());
1174  return ie.Get_I1() + ie.Get_I2()/ie.Get_I3() - 6.0;
1175 }
1176 
1177 void TMOP_Metric_321::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1178 {
1179  // mu_321 = I1 + I2/I3b^2 - 6 = I1 + I2/I3 - 6
1180  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
1181  ie.SetJacobian(Jpt.GetData());
1182  const double I3 = ie.Get_I3();
1183  Add(1.0/I3, ie.Get_dI2(),
1184  -2*ie.Get_I2()/(I3*ie.Get_I3b()), ie.Get_dI3b(), P);
1185  P += ie.Get_dI1();
1186 }
1187 
1188 void TMOP_Metric_321::AssembleH(const DenseMatrix &Jpt,
1189  const DenseMatrix &DS,
1190  const double weight,
1191  DenseMatrix &A) const
1192 {
1193  // P = dI1 + (1/I3)*dI2 - (2*I2/I3b^3)*dI3b
1194  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b) + (1/I3)*ddI2 + (dI3b x dz) + z*ddI3b
1195  //
1196  // z = -2*I2/I3b^3
1197  // dz = (-2/I3b^3)*dI2 + (2*I2)*(3/I3b^4)*dI3b
1198  //
1199  // dP = ddI1 + (-2/I3b^3)*(dI2 x dI3b + dI3b x dI2) + (1/I3)*ddI2
1200  // + (6*I2/I3b^4)*(dI3b x dI3b) + (-2*I2/I3b^3)*ddI3b
1201  ie.SetJacobian(Jpt.GetData());
1202  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1203  const double c0 = 1.0/ie.Get_I3b();
1204  const double c1 = weight*c0*c0;
1205  const double c2 = -2*c0*c1;
1206  const double c3 = c2*ie.Get_I2();
1207  ie.Assemble_ddI1(weight, A.GetData());
1208  ie.Assemble_ddI2(c1, A.GetData());
1209  ie.Assemble_ddI3b(c3, A.GetData());
1210  ie.Assemble_TProd(c2, ie.Get_dI2(), ie.Get_dI3b(), A.GetData());
1211  ie.Assemble_TProd(-3*c0*c3, ie.Get_dI3b(), A.GetData());
1212 }
1213 
1214 double TMOP_Metric_322::EvalWMatrixForm(const DenseMatrix &Jpt) const
1215 {
1216  // mu_322 = 1 / (6 det(J)) |J - adj(J)^t|^2
1217  DenseMatrix adj_J_t(3);
1218  CalcAdjugateTranspose(Jpt, adj_J_t);
1219  adj_J_t *= -1.0;
1220  adj_J_t.Add(1.0, Jpt);
1221  return 1.0 / 6.0 / Jpt.Det() * adj_J_t.FNorm2();
1222 }
1223 
1224 double TMOP_Metric_322::EvalW(const DenseMatrix &Jpt) const
1225 {
1226  // mu_322 = 1 / (6 det(J)) |J - adj(J)^t|^2
1227  // = 1 / (6 det(J)) |J|^2 + 1/6 det(J) |J^{-1}|^2 - 1
1228  // = I1b / (I3b^-1/3) / 6 + I2b (I3b^1/3) / 6 - 1
1229  ie.SetJacobian(Jpt.GetData());
1230 
1231  return ie.Get_I1b() / pow(ie.Get_I3b(), 1.0/3.0) / 6.0 +
1232  ie.Get_I2b() * pow(ie.Get_I3b(), 1.0/3.0) / 6.0 - 1.0;
1233 }
1234 
1235 void TMOP_Metric_322::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1236 {
1237  // mu_322 = I1b (I3b^-1/3) / 6 + I2b (I3b^1/3) / 6 - 1
1238  // P = 1/6 (I3b^-1/3) dI1b - 1/18 I1b (I3b^-4/3) dI3b
1239  // + 1/6 (I3b^1/3) dI2b + 1/18 I2b (I3b^-2/3) dI3b
1240  ie.SetJacobian(Jpt.GetData());
1241  P.Set(1.0/6.0 * pow(ie.Get_I3b(), -1.0/3.0),
1242  ie.Get_dI1b());
1243  P.Add(-1.0/18.0 * ie.Get_I1b() * pow(ie.Get_I3b(), -4.0/3.0),
1244  ie.Get_dI3b());
1245  P.Add(1.0/6.0 * pow(ie.Get_I3b(), 1.0/3.0),
1246  ie.Get_dI2b());
1247  P.Add(1.0/18.0 * ie.Get_I2b() * pow(ie.Get_I3b(), -2.0/3.0),
1248  ie.Get_dI3b());
1249 }
1250 
1251 void TMOP_Metric_322::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
1252  const double weight, DenseMatrix &A) const
1253 {
1254  // P = 1/6 (I3b^-1/3) dI1b - 1/18 I1b (I3b^-4/3) dI3b
1255  // + 1/6 (I3b^1/3) dI2b + 1/18 I2b (I3b^-2/3) dI3b
1256  // dP = 1/6 (I3b^-1/3) ddI1b - 1/18 (I3b^-4/3) (dI1b x dI3b)
1257  // - 1/18 I1b (I3b^-4/3) ddI3b
1258  // - 1/18 (I3b^-4/3) (dI3b x dI1b)
1259  // + 2/27 I1b (I3b^-7/3) (dI3b x dI3b)
1260  // + 1/6 (I3b^1/3) ddI2b + 1/18 (I3b^-2/3) (dI2b x dI3b)
1261  // + 1/18 I2b (I3b^-2/3) ddI3b
1262  // + 1/18 (I3b^-2/3) (dI3b x dI2b)
1263  // - 1/27 I2b (I3b^-5/3) (dI3b x dI3b)
1264  ie.SetJacobian(Jpt.GetData());
1265  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1266  const double p13 = weight * pow(ie.Get_I3b(), 1.0/3.0),
1267  m13 = weight * pow(ie.Get_I3b(), -1.0/3.0),
1268  m23 = weight * pow(ie.Get_I3b(), -2.0/3.0),
1269  m43 = weight * pow(ie.Get_I3b(), -4.0/3.0),
1270  m53 = weight * pow(ie.Get_I3b(), -5.0/3.0),
1271  m73 = weight * pow(ie.Get_I3b(), -7.0/3.0);
1272  ie.Assemble_ddI1b(1.0/6.0 * m13, A.GetData());
1273  // Combines - 1/18 (I3b^-4/3) (dI1b x dI3b) - 1/18 (I3b^-4/3) (dI3b x dI1b).
1274  ie.Assemble_TProd(-1.0/18.0 * m43,
1275  ie.Get_dI1b(), ie.Get_dI3b(), A.GetData());
1276  ie.Assemble_ddI3b(-1.0/18.0 * ie.Get_I1b() * m43, A.GetData());
1277  ie.Assemble_TProd(2.0/27.0 * ie.Get_I1b() * m73,
1278  ie.Get_dI3b(), A.GetData());
1279  ie.Assemble_ddI2b(1.0/6.0 * p13, A.GetData());
1280  // Combines + 1/18 (I3b^-2/3) (dI2b x dI3b) + 1/18 (I3b^-2/3) (dI3b x dI2b).
1281  ie.Assemble_TProd(1.0/18.0 * m23,
1282  ie.Get_dI2b(), ie.Get_dI3b(), A.GetData());
1283  ie.Assemble_ddI3b(1.0/18.0 * ie.Get_I2b() * m23, A.GetData());
1284  ie.Assemble_TProd(-1.0/27.0 * ie.Get_I2b() * m53,
1285  ie.Get_dI3b(), A.GetData());
1286 }
1287 
1288 double TMOP_Metric_323::EvalWMatrixForm(const DenseMatrix &Jpt) const
1289 {
1290  // mu_323 = |J|^3 - 3 sqrt(3) ln(det(J)) - 3 sqrt(3).
1291  double fnorm = Jpt.FNorm();
1292  return fnorm * fnorm * fnorm - 3.0 * sqrt(3.0) * (log(Jpt.Det()) + 1.0);
1293 }
1294 
1295 double TMOP_Metric_323::EvalW(const DenseMatrix &Jpt) const
1296 {
1297  // mu_323 = I1^3/2 - 3 sqrt(3) ln(I3b) - 3 sqrt(3).
1298  ie.SetJacobian(Jpt.GetData());
1299  return pow(ie.Get_I1(), 1.5) - 3.0 * sqrt(3.0) * (log(ie.Get_I3b()) + 1.0);
1300 }
1301 
1302 void TMOP_Metric_323::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1303 {
1304  // mu_323 = I1^3/2 - 3 sqrt(3) ln(I3b) - 3 sqrt(3).
1305  // P = 3/2 (I1^1/2) dI1 - 3 sqrt(3) (I3b^-1) dI3b.
1306  ie.SetJacobian(Jpt.GetData());
1307  P.Set(1.5 * sqrt(ie.Get_I1()), ie.Get_dI1());
1308  P.Add(- 3.0 * sqrt(3.0) / ie.Get_I3b(), ie.Get_dI3b());
1309 }
1310 
1311 void TMOP_Metric_323::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
1312  const double weight, DenseMatrix &A) const
1313 {
1314  // P = 3/2 (I1^1/2) dI1 - 3 sqrt(3) (I3b^-1) dI3b
1315  // dP = 3/2 (I1^1/2) ddI1 + 3/4 (I1^-1/2) (dI1 x dI1)
1316  // - 3 sqrt(3) (I3b^-1) ddI3b + 3 sqrt(3) (I3b^-2) (dI3b x dI3b)
1317  ie.SetJacobian(Jpt.GetData());
1318  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1319  ie.Assemble_ddI1(weight * 1.5 * sqrt(ie.Get_I1()), A.GetData());
1320  ie.Assemble_TProd(weight * 0.75 / sqrt(ie.Get_I1()),
1321  ie.Get_dI1(), A.GetData());
1322  ie.Assemble_ddI3b(- weight * 3.0 * sqrt(3.0) / ie.Get_I3b(), A.GetData());
1323  ie.Assemble_TProd(weight * 3.0 * sqrt(3.0) / ie.Get_I3b() / ie.Get_I3b(),
1324  ie.Get_dI3b(), A.GetData());
1325 }
1326 
1327 double TMOP_Metric_352::EvalW(const DenseMatrix &Jpt) const
1328 {
1329  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
1330  ie.SetJacobian(Jpt.GetData());
1331  const double I3b = ie.Get_I3b();
1332  return 0.5*(I3b - 1.0)*(I3b - 1.0)/(I3b - tau0);
1333 }
1334 
1335 void TMOP_Metric_352::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1336 {
1337  // mu_352 = 0.5*(det(J) - 1)^2 / (det(J) - tau0)
1338  // P = (I3b - 1)/(I3b - tau0)*dI3b + 0.5*(I3b - 1)^2*(-1/(I3b - tau0)^2)*dI3b
1339  // = [ (I3b - 1)/(I3b - tau0) - 0.5*(I3b - 1)^2/(I3b - tau0)^2 ] * dI3b
1340  // = (c - 0.5*c*c) * dI3b
1341  ie.SetJacobian(Jpt.GetData());
1342  const double I3b = ie.Get_I3b();
1343  const double c = (I3b - 1.0)/(I3b - tau0);
1344  P.Set(c - 0.5*c*c, ie.Get_dI3b());
1345 }
1346 
1347 void TMOP_Metric_352::AssembleH(const DenseMatrix &Jpt,
1348  const DenseMatrix &DS,
1349  const double weight,
1350  DenseMatrix &A) const
1351 {
1352  // c = (I3b - 1)/(I3b - tau0)
1353  //
1354  // P = (c - 0.5*c*c) * dI3b
1355  // dP = (1 - c)*(dI3b x dc) + (c - 0.5*c*c)*ddI3b
1356  //
1357  // dc = 1/(I3b - tau0)*dI3b - (I3b - 1)/(I3b - tau)^2*dI3b =
1358  // = (1 - c)/(I3b - tau0)*dI3b
1359  //
1360  // dP = (1 - c)^2/(I3b - tau0)*(dI3b x dI3b) + (c - 0.5*c*c)*ddI3b
1361  ie.SetJacobian(Jpt.GetData());
1362  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1363  const double I3b = ie.Get_I3b();
1364  const double c0 = 1.0/(I3b - tau0);
1365  const double c = c0*(I3b - 1.0);
1366  ie.Assemble_TProd(weight*c0*(1.0 - c)*(1.0 - c), ie.Get_dI3b(), A.GetData());
1367  ie.Assemble_ddI3b(weight*(c - 0.5*c*c), A.GetData());
1368 }
1369 
1370 
1371 
1372 
1373 double TMOP_Metric_360::EvalWMatrixForm(const DenseMatrix &Jpt) const
1374 {
1375  // mu_360 = |J|^3 / 3^(3/2) - det(J)
1376  const double fnorm = Jpt.FNorm();
1377  return fnorm * fnorm * fnorm / pow(3.0, 1.5) - Jpt.Det();
1378 }
1379 
1380 double TMOP_Metric_360::EvalW(const DenseMatrix &Jpt) const
1381 {
1382  // mu_360 = (I1/3)^(3/2) - I3b.
1383  ie.SetJacobian(Jpt.GetData());
1384  return pow(ie.Get_I1()/3.0, 1.5) - ie.Get_I3b();
1385 }
1386 
1387 void TMOP_Metric_360::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
1388 {
1389  // mu_360 = (I1/3)^(3/2) - I3b.
1390  // P = 3/2 * (I1/3)^1/2 * dI1 / 3 - dI3b
1391  // = 1/2 * (I1/3)^1/2 * dI1 - dI3b.
1392  ie.SetJacobian(Jpt.GetData());
1393  Add(0.5 * sqrt(ie.Get_I1()/3.0), ie.Get_dI1(), -1.0, ie.Get_dI3b(), P);
1394 }
1395 
1396 void TMOP_Metric_360::AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
1397  const double weight, DenseMatrix &A) const
1398 {
1399  // P = 1/2 * (I1/3)^1/2 * dI1 - dI3b.
1400  // dP = 1/12 * (I1/3)^(-1/2) * (dI1 x dI1) + 1/2 * (I1/3)^1/2 * ddI1 - ddI3b
1401  ie.SetJacobian(Jpt.GetData());
1402  ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
1403  ie.Assemble_TProd(weight / 12.0 / sqrt(ie.Get_I1()/3.0),
1404  ie.Get_dI1(), A.GetData());
1405  ie.Assemble_ddI1(weight / 2.0 * sqrt(ie.Get_I1()/3.0), A.GetData());
1406  ie.Assemble_ddI3b(-weight, A.GetData());
1407 }
1408 
1409 double TMOP_AMetric_011::EvalW(const DenseMatrix &Jpt) const
1410 {
1411  MFEM_VERIFY(Jtr != NULL,
1412  "Requires a target Jacobian, use SetTargetJacobian().");
1413 
1414  int dim = Jpt.Size();
1415 
1416  DenseMatrix Jpr(dim, dim);
1417  Mult(Jpt, *Jtr, Jpr);
1418 
1419  double alpha = Jpr.Det(),
1420  omega = Jtr->Det();
1421 
1422  DenseMatrix AdjAt(dim), WtW(dim), WRK(dim), Jtrt(dim);
1423  CalcAdjugateTranspose(Jpr, AdjAt);
1424  Jtrt.Transpose(*Jtr);
1425  MultAAt(Jtrt, WtW);
1426  WtW *= 1./omega;
1427  Mult(AdjAt, WtW, WRK);
1428 
1429  WRK -= Jpr;
1430  WRK *= -1.;
1431 
1432  return (0.25/alpha)*WRK.FNorm2();
1433 }
1434 
1435 double TMOP_AMetric_014a::EvalW(const DenseMatrix &Jpt) const
1436 {
1437  MFEM_VERIFY(Jtr != NULL,
1438  "Requires a target Jacobian, use SetTargetJacobian().");
1439 
1440  int dim = Jpt.Size();
1441 
1442  DenseMatrix Jpr(dim, dim);
1443  Mult(Jpt, *Jtr, Jpr);
1444 
1445  double sqalpha = pow(Jpr.Det(), 0.5),
1446  sqomega = pow(Jtr->Det(), 0.5);
1447 
1448  return 0.5*pow(sqalpha/sqomega - sqomega/sqalpha, 2.);
1449 }
1450 
1451 double TMOP_AMetric_036::EvalW(const DenseMatrix &Jpt) const
1452 {
1453  MFEM_VERIFY(Jtr != NULL,
1454  "Requires a target Jacobian, use SetTargetJacobian().");
1455 
1456  int dim = Jpt.Size();
1457 
1458  DenseMatrix Jpr(dim, dim);
1459  Mult(Jpt, *Jtr, Jpr); // T*W = A
1460 
1461  double alpha = Jpr.Det(); // det(A)
1462  Jpr -= *Jtr; // A-W
1463 
1464  return (1./alpha)*(Jpr.FNorm2()); //(1/alpha)*(|A-W|^2)
1465 }
1466 
1467 double TMOP_AMetric_107a::EvalW(const DenseMatrix &Jpt) const
1468 {
1469  MFEM_VERIFY(Jtr != NULL,
1470  "Requires a target Jacobian, use SetTargetJacobian().");
1471 
1472  int dim = Jpt.Size();
1473 
1474  DenseMatrix Jpr(dim, dim);
1475  Mult(Jpt, *Jtr, Jpr);
1476 
1477  double alpha = Jpr.Det(),
1478  aw = Jpr.FNorm()/Jtr->FNorm();
1479 
1480  DenseMatrix W = *Jtr;
1481  W *= aw;
1482  Jpr -= W;
1483 
1484  return (0.5/alpha)*Jpr.FNorm2();
1485 }
1486 
1487 
1488 void TargetConstructor::ComputeAvgVolume() const
1489 {
1490  MFEM_VERIFY(nodes, "Nodes are not given!");
1491  MFEM_ASSERT(avg_volume == 0.0, "The average volume is already computed!");
1492 
1493  Mesh *mesh = nodes->FESpace()->GetMesh();
1494  const int NE = mesh->GetNE();
1495  IsoparametricTransformation Tr;
1496  double volume = 0.0;
1497 
1498  for (int i = 0; i < NE; i++)
1499  {
1500  mesh->GetElementTransformation(i, *nodes, &Tr);
1501  const IntegrationRule &ir =
1502  IntRules.Get(mesh->GetElementBaseGeometry(i), Tr.OrderJ());
1503  for (int j = 0; j < ir.GetNPoints(); j++)
1504  {
1505  const IntegrationPoint &ip = ir.IntPoint(j);
1506  Tr.SetIntPoint(&ip);
1507  volume += ip.weight * Tr.Weight();
1508  }
1509  }
1510 
1511  NCMesh *ncmesh = mesh->ncmesh;
1512  if (Parallel() == false)
1513  {
1514  avg_volume = (ncmesh == NULL) ?
1515  volume / NE : volume / ncmesh->GetNumRootElements();
1516 
1517  }
1518 #ifdef MFEM_USE_MPI
1519  else
1520  {
1521  double area_NE[4];
1522  area_NE[0] = volume; area_NE[1] = NE;
1523  MPI_Allreduce(area_NE, area_NE + 2, 2, MPI_DOUBLE, MPI_SUM, comm);
1524  avg_volume = (ncmesh == NULL) ?
1525  area_NE[2] / area_NE[3] : area_NE[2] / ncmesh->GetNumRootElements();
1526  }
1527 #endif
1528 }
1529 
1530 void TargetConstructor::ComputeAllElementTargets_Fallback(
1531  const FiniteElementSpace &fes,
1532  const IntegrationRule &ir,
1533  const Vector &xe,
1534  DenseTensor &Jtr) const
1535 {
1536  // Fallback to the 1-element method, ComputeElementTargets()
1537 
1538  // When UsesPhysicalCoordinates() == true, we assume 'xe' uses
1539  // ElementDofOrdering::LEXICOGRAPHIC iff 'fe' is a TensorFiniteElement.
1540 
1541  const Mesh *mesh = fes.GetMesh();
1542  const int NE = mesh->GetNE();
1543  // Quick return for empty processors:
1544  if (NE == 0) { return; }
1545  const int dim = mesh->Dimension();
1546  MFEM_VERIFY(mesh->GetNumGeometries(dim) <= 1,
1547  "mixed meshes are not supported");
1548  MFEM_VERIFY(!fes.IsVariableOrder(), "variable orders are not supported");
1549  const FiniteElement &fe = *fes.GetFE(0);
1550  const int sdim = fes.GetVDim();
1551  const int nvdofs = sdim*fe.GetDof();
1552  MFEM_VERIFY(!UsesPhysicalCoordinates() ||
1553  xe.Size() == NE*nvdofs, "invalid input Vector 'xe'!");
1554  const int NQ = ir.GetNPoints();
1555  const Array<int> *dof_map = nullptr;
1556  if (UsesPhysicalCoordinates())
1557  {
1558  const TensorBasisElement *tfe =
1559  dynamic_cast<const TensorBasisElement *>(&fe);
1560  if (tfe)
1561  {
1562  dof_map = &tfe->GetDofMap();
1563  if (dof_map->Size() == 0) { dof_map = nullptr; }
1564  }
1565  }
1566 
1567  Vector elfun_lex, elfun_nat;
1568  DenseTensor J;
1569  xe.HostRead();
1570  Jtr.HostWrite();
1571  if (UsesPhysicalCoordinates() && dof_map != nullptr)
1572  {
1573  elfun_nat.SetSize(nvdofs);
1574  }
1575  for (int e = 0; e < NE; e++)
1576  {
1577  if (UsesPhysicalCoordinates())
1578  {
1579  if (!dof_map)
1580  {
1581  elfun_nat.SetDataAndSize(xe.GetData()+e*nvdofs, nvdofs);
1582  }
1583  else
1584  {
1585  elfun_lex.SetDataAndSize(xe.GetData()+e*nvdofs, nvdofs);
1586  const int ndofs = fe.GetDof();
1587  for (int d = 0; d < sdim; d++)
1588  {
1589  for (int i_lex = 0; i_lex < ndofs; i_lex++)
1590  {
1591  elfun_nat[(*dof_map)[i_lex]+d*ndofs] =
1592  elfun_lex[i_lex+d*ndofs];
1593  }
1594  }
1595  }
1596  }
1597  J.UseExternalData(Jtr(e*NQ).Data(), sdim, dim, NQ);
1598  ComputeElementTargets(e, fe, ir, elfun_nat, J);
1599  }
1600 }
1601 
1602 bool TargetConstructor::ContainsVolumeInfo() const
1603 {
1604  switch (target_type)
1605  {
1606  case IDEAL_SHAPE_UNIT_SIZE: return false;
1607  case IDEAL_SHAPE_EQUAL_SIZE:
1608  case IDEAL_SHAPE_GIVEN_SIZE:
1609  case GIVEN_SHAPE_AND_SIZE:
1610  case GIVEN_FULL: return true;
1611  default: MFEM_ABORT("TargetType not added to ContainsVolumeInfo.");
1612  }
1613  return false;
1614 }
1615 
1616 void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
1617  const IntegrationRule &ir,
1618  const Vector &elfun,
1619  DenseTensor &Jtr) const
1620 {
1621  MFEM_CONTRACT_VAR(elfun);
1622  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1623 
1624  const FiniteElement *nfe = (target_type != IDEAL_SHAPE_UNIT_SIZE) ?
1625  nodes->FESpace()->GetFE(e_id) : NULL;
1626  const DenseMatrix &Wideal =
1627  Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
1628  MFEM_ASSERT(Wideal.Height() == Jtr.SizeI(), "");
1629  MFEM_ASSERT(Wideal.Width() == Jtr.SizeJ(), "");
1630 
1631  switch (target_type)
1632  {
1633  case IDEAL_SHAPE_UNIT_SIZE:
1634  {
1635  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = Wideal; }
1636  break;
1637  }
1638  case IDEAL_SHAPE_EQUAL_SIZE:
1639  {
1640  if (avg_volume == 0.0) { ComputeAvgVolume(); }
1641  DenseMatrix W(Wideal.Height());
1642 
1643  NCMesh *ncmesh = nodes->FESpace()->GetMesh()->ncmesh;
1644  double el_volume = avg_volume;
1645  if (ncmesh)
1646  {
1647  el_volume = avg_volume / ncmesh->GetElementSizeReduction(e_id);
1648  }
1649 
1650  W.Set(std::pow(volume_scale * el_volume / Wideal.Det(),
1651  1./W.Height()), Wideal);
1652  for (int i = 0; i < ir.GetNPoints(); i++) { Jtr(i) = W; }
1653  break;
1654  }
1655  case IDEAL_SHAPE_GIVEN_SIZE:
1656  case GIVEN_SHAPE_AND_SIZE:
1657  {
1658  const int dim = nfe->GetDim(), dof = nfe->GetDof();
1659  MFEM_ASSERT(dim == nodes->FESpace()->GetVDim(), "");
1660  DenseMatrix dshape(dof, dim), pos(dof, dim);
1661  Array<int> xdofs(dof * dim);
1662  Vector posV(pos.Data(), dof * dim);
1663  double detW;
1664 
1665  // always initialize detW to suppress a warning:
1666  detW = (target_type == IDEAL_SHAPE_GIVEN_SIZE) ? Wideal.Det() : 0.0;
1667  nodes->FESpace()->GetElementVDofs(e_id, xdofs);
1668  nodes->GetSubVector(xdofs, posV);
1669  for (int i = 0; i < ir.GetNPoints(); i++)
1670  {
1671  nfe->CalcDShape(ir.IntPoint(i), dshape);
1672  MultAtB(pos, dshape, Jtr(i));
1673  if (target_type == IDEAL_SHAPE_GIVEN_SIZE)
1674  {
1675  const double det = Jtr(i).Det();
1676  MFEM_VERIFY(det > 0.0, "The given mesh is inverted!");
1677  Jtr(i).Set(std::pow(det / detW, 1./dim), Wideal);
1678  }
1679  }
1680  break;
1681  }
1682  default:
1683  MFEM_ABORT("invalid target type!");
1684  }
1685 }
1686 
1687 void TargetConstructor::ComputeElementTargetsGradient(
1688  const IntegrationRule &ir,
1689  const Vector &elfun,
1690  IsoparametricTransformation &Tpr,
1691  DenseTensor &dJtr) const
1692 {
1693  MFEM_CONTRACT_VAR(elfun);
1694  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
1695 
1696  // TODO: Compute derivative for targets with GIVEN_SHAPE or/and GIVEN_SIZE
1697  for (int i = 0; i < Tpr.GetFE()->GetDim()*ir.GetNPoints(); i++)
1698  { dJtr(i) = 0.; }
1699 }
1700 
1701 void AnalyticAdaptTC::SetAnalyticTargetSpec(Coefficient *sspec,
1702  VectorCoefficient *vspec,
1703  TMOPMatrixCoefficient *mspec)
1704 {
1705  scalar_tspec = sspec;
1706  vector_tspec = vspec;
1707  matrix_tspec = mspec;
1708 }
1709 
1710 void AnalyticAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
1711  const IntegrationRule &ir,
1712  const Vector &elfun,
1713  DenseTensor &Jtr) const
1714 {
1715  DenseMatrix point_mat;
1716  point_mat.UseExternalData(elfun.GetData(), fe.GetDof(), fe.GetDim());
1717 
1718  switch (target_type)
1719  {
1720  case GIVEN_FULL:
1721  {
1722  MFEM_VERIFY(matrix_tspec != NULL,
1723  "Target type GIVEN_FULL requires a MatrixCoefficient.");
1724 
1725  IsoparametricTransformation Tpr;
1726  Tpr.SetFE(&fe);
1727  Tpr.ElementNo = e_id;
1728  Tpr.ElementType = ElementTransformation::ELEMENT;
1729  Tpr.GetPointMat().Transpose(point_mat);
1730 
1731  for (int i = 0; i < ir.GetNPoints(); i++)
1732  {
1733  const IntegrationPoint &ip = ir.IntPoint(i);
1734  Tpr.SetIntPoint(&ip);
1735  matrix_tspec->Eval(Jtr(i), Tpr, ip);
1736  }
1737  break;
1738  }
1739  default:
1740  MFEM_ABORT("Incompatible target type for analytic adaptation!");
1741  }
1742 }
1743 
1744 void AnalyticAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
1745  const Vector &elfun,
1746  IsoparametricTransformation &Tpr,
1747  DenseTensor &dJtr) const
1748 {
1749  const FiniteElement *fe = Tpr.GetFE();
1750  DenseMatrix point_mat;
1751  point_mat.UseExternalData(elfun.GetData(), fe->GetDof(), fe->GetDim());
1752 
1753  switch (target_type)
1754  {
1755  case GIVEN_FULL:
1756  {
1757  MFEM_VERIFY(matrix_tspec != NULL,
1758  "Target type GIVEN_FULL requires a TMOPMatrixCoefficient.");
1759 
1760  for (int d = 0; d < fe->GetDim(); d++)
1761  {
1762  for (int i = 0; i < ir.GetNPoints(); i++)
1763  {
1764  const IntegrationPoint &ip = ir.IntPoint(i);
1765  Tpr.SetIntPoint(&ip);
1766  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
1767  matrix_tspec->EvalGrad(dJtr_i, Tpr, ip, d);
1768  }
1769  }
1770  break;
1771  }
1772  default:
1773  MFEM_ABORT("Incompatible target type for analytic adaptation!");
1774  }
1775 }
1776 
1777 namespace internal
1778 {
1779 
1780 // mfem::forall-based copy kernel -- used by protected methods below.
1781 // Needed as a workaround for the nvcc restriction that methods with mfem::forall
1782 // in them must to be public.
1783 static inline void device_copy(double *d_dest, const double *d_src, int size)
1784 {
1785  mfem::forall(size, [=] MFEM_HOST_DEVICE (int i) { d_dest[i] = d_src[i]; });
1786 }
1787 
1788 } // namespace internal
1789 
1790 #ifdef MFEM_USE_MPI
1791 void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction &t)
1792 {
1793  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1794  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1795 
1796  ParFiniteElementSpace *ptspec_fes = t.ParFESpace();
1797 
1798  tspec_sav = tspec;
1799 
1800  delete tspec_fesv;
1801  tspec_fesv = new FiniteElementSpace(ptspec_fes->GetMesh(),
1802  ptspec_fes->FEColl(), ncomp);
1803 
1804  delete ptspec_fesv;
1805  ptspec_fesv = new ParFiniteElementSpace(ptspec_fes->GetParMesh(),
1806  ptspec_fes->FEColl(), ncomp);
1807 
1808  delete tspec_pgf;
1809  tspec_pgf = new ParGridFunction(ptspec_fesv, tspec);
1810  tspec_gf = tspec_pgf;
1811 
1812  adapt_eval->SetParMetaInfo(*ptspec_fes->GetParMesh(), *ptspec_fesv);
1813  adapt_eval->SetInitialField(*ptspec_fes->GetMesh()->GetNodes(), tspec);
1814 }
1815 
1816 void DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange()
1817 {
1818  ptspec_fesv->Update();
1819  if (tspec_fesv)
1820  {
1821  delete tspec_fesv;
1822  tspec_fesv = new FiniteElementSpace(ptspec_fesv->GetMesh(),
1823  ptspec_fesv->FEColl(), ncomp);
1824  }
1825  tspec_pgf->Update();
1826  tspec_gf = tspec_pgf;
1827  tspec.SetDataAndSize(tspec_pgf->GetData(), tspec_pgf->Size());
1828  tspec_sav = tspec;
1829 
1830  adapt_eval->SetParMetaInfo(*ptspec_fesv->GetParMesh(), *ptspec_fesv);
1831  adapt_eval->SetInitialField(*ptspec_fesv->GetMesh()->GetNodes(), tspec);
1832 }
1833 
1834 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const ParGridFunction &tspec_)
1835 {
1836  const int vdim = tspec_.FESpace()->GetVDim(),
1837  ndof = tspec_.FESpace()->GetNDofs();
1838  MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
1839 
1840  const auto tspec__d = tspec_.Read();
1841  auto tspec_d = tspec.ReadWrite();
1842  const int offset = idx*ndof;
1843  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1844  FinalizeParDiscreteTargetSpec(tspec_);
1845 }
1846 
1847 void DiscreteAdaptTC::SetParDiscreteTargetSize(const ParGridFunction &tspec_)
1848 {
1849  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1850  "Discrete target size should be ordered byNodes.");
1851  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1852  sizeidx = ncomp;
1853  SetDiscreteTargetBase(tspec_);
1854  FinalizeParDiscreteTargetSpec(tspec_);
1855 }
1856 
1857 void DiscreteAdaptTC::SetParDiscreteTargetSkew(const ParGridFunction &tspec_)
1858 {
1859  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1860  "Discrete target skewness should be ordered byNodes.");
1861  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1862  skewidx = ncomp;
1863  SetDiscreteTargetBase(tspec_);
1864  FinalizeParDiscreteTargetSpec(tspec_);
1865 }
1866 
1867 void DiscreteAdaptTC::SetParDiscreteTargetAspectRatio(const ParGridFunction &ar)
1868 {
1869  MFEM_VERIFY(ar.FESpace()->GetOrdering() == Ordering::byNODES,
1870  "Discrete target aspect ratio should be ordered byNodes.");
1871  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
1872  aspectratioidx = ncomp;
1873  SetDiscreteTargetBase(ar);
1874  FinalizeParDiscreteTargetSpec(ar);
1875 }
1876 
1877 void DiscreteAdaptTC::SetParDiscreteTargetOrientation(const ParGridFunction &o)
1878 {
1879  MFEM_VERIFY(o.FESpace()->GetOrdering() == Ordering::byNODES,
1880  "Discrete target orientation should be ordered byNodes.");
1881  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
1882  orientationidx = ncomp;
1883  SetDiscreteTargetBase(o);
1884  FinalizeParDiscreteTargetSpec(o);
1885 }
1886 
1887 void DiscreteAdaptTC::SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
1888 {
1889  SetParDiscreteTargetSize(tspec_);
1890 }
1891 #endif // MFEM_USE_MPI
1892 
1893 void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
1894 {
1895  const int vdim = tspec_.FESpace()->GetVDim(),
1896  ndof = tspec_.FESpace()->GetNDofs();
1897 
1898  ncomp += vdim;
1899 
1900  // need to append data to tspec
1901  // make a copy of tspec->tspec_temp, increase its size, and
1902  // copy data from tspec_temp -> tspec, then add new entries
1903  Vector tspec_temp = tspec;
1904  tspec.UseDevice(true);
1905  tspec_sav.UseDevice(true);
1906  tspec.SetSize(ncomp*ndof);
1907 
1908  const auto tspec_temp_d = tspec_temp.Read();
1909  auto tspec_d = tspec.ReadWrite();
1910  internal::device_copy(tspec_d, tspec_temp_d, tspec_temp.Size());
1911 
1912  const auto tspec__d = tspec_.Read();
1913  const int offset = (ncomp-vdim)*ndof;
1914  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1915 }
1916 
1917 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
1918 {
1919  const int vdim = tspec_.FESpace()->GetVDim(),
1920  ndof = tspec_.FESpace()->GetNDofs();
1921  MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTargetSpec.");
1922 
1923  const auto tspec__d = tspec_.Read();
1924  auto tspec_d = tspec.ReadWrite();
1925  const int offset = idx*ndof;
1926  internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1927  FinalizeSerialDiscreteTargetSpec(tspec_);
1928 }
1929 
1930 void DiscreteAdaptTC::SetSerialDiscreteTargetSize(const GridFunction &tspec_)
1931 {
1932  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1933  "Discrete target size should be ordered byNodes.");
1934  if (sizeidx > -1) { SetTspecAtIndex(sizeidx, tspec_); return; }
1935  sizeidx = ncomp;
1936  SetDiscreteTargetBase(tspec_);
1937  FinalizeSerialDiscreteTargetSpec(tspec_);
1938 }
1939 
1940 void DiscreteAdaptTC::SetSerialDiscreteTargetSkew(const GridFunction &tspec_)
1941 {
1942  MFEM_VERIFY(tspec_.FESpace()->GetOrdering() == Ordering::byNODES,
1943  "Discrete target skewness should be ordered byNodes.");
1944  if (skewidx > -1) { SetTspecAtIndex(skewidx, tspec_); return; }
1945  skewidx = ncomp;
1946  SetDiscreteTargetBase(tspec_);
1947  FinalizeSerialDiscreteTargetSpec(tspec_);
1948 }
1949 
1950 void DiscreteAdaptTC::SetSerialDiscreteTargetAspectRatio(const GridFunction &ar)
1951 {
1952  MFEM_VERIFY(ar.FESpace()->GetOrdering() == Ordering::byNODES,
1953  "Discrete target aspect ratio should be ordered byNodes.");
1954  if (aspectratioidx > -1) { SetTspecAtIndex(aspectratioidx, ar); return; }
1955  aspectratioidx = ncomp;
1956  SetDiscreteTargetBase(ar);
1957  FinalizeSerialDiscreteTargetSpec(ar);
1958 }
1959 
1960 void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(const GridFunction &o)
1961 {
1962  MFEM_VERIFY(o.FESpace()->GetOrdering() == Ordering::byNODES,
1963  "Discrete target orientation should be ordered byNodes.");
1964  if (orientationidx > -1) { SetTspecAtIndex(orientationidx, o); return; }
1965  orientationidx = ncomp;
1966  SetDiscreteTargetBase(o);
1967  FinalizeSerialDiscreteTargetSpec(o);
1968 }
1969 
1970 void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec(const GridFunction &t)
1971 {
1972  MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
1973  MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
1974 
1975  const FiniteElementSpace *tspec_fes = t.FESpace();
1976 
1977  tspec_sav = tspec;
1978 
1979  delete tspec_fesv;
1980  tspec_fesv = new FiniteElementSpace(tspec_fes->GetMesh(),
1981  tspec_fes->FEColl(), ncomp,
1982  Ordering::byNODES);
1983 
1984  delete tspec_gf;
1985  tspec_gf = new GridFunction(tspec_fesv, tspec);
1986 
1987  adapt_eval->SetSerialMetaInfo(*tspec_fes->GetMesh(), *tspec_fesv);
1988  adapt_eval->SetInitialField(*tspec_fes->GetMesh()->GetNodes(), tspec);
1989 }
1990 
1991 void DiscreteAdaptTC::GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
1992 {
1993  if (idx < 0) { return; }
1994  const int ndof = tspec_.FESpace()->GetNDofs(),
1995  vdim = tspec_.FESpace()->GetVDim();
1996  MFEM_VERIFY(ndof == tspec.Size()/ncomp,
1997  "Inconsistency in GetSerialDiscreteTargetSpec.");
1998 
1999  for (int i = 0; i < ndof*vdim; i++)
2000  {
2001  tspec_(i) = tspec(i + idx*ndof);
2002  }
2003 }
2004 
2005 void DiscreteAdaptTC::UpdateAfterMeshTopologyChange()
2006 {
2007  tspec_fesv->Update();
2008  tspec_gf->Update();
2009  tspec.SetDataAndSize(tspec_gf->GetData(), tspec_gf->Size());
2010  tspec_sav = tspec;
2011 
2012  adapt_eval->SetSerialMetaInfo(*tspec_fesv->GetMesh(), *tspec_fesv);
2013  adapt_eval->SetInitialField(*tspec_fesv->GetMesh()->GetNodes(), tspec);
2014 }
2015 
2016 void DiscreteAdaptTC::SetSerialDiscreteTargetSpec(const GridFunction &tspec_)
2017 {
2018  SetSerialDiscreteTargetSize(tspec_);
2019 }
2020 
2021 
2022 void DiscreteAdaptTC::UpdateTargetSpecification(const Vector &new_x,
2023  bool reuse_flag,
2024  int new_x_ordering)
2025 {
2026  if (reuse_flag && good_tspec) { return; }
2027 
2028  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
2029  adapt_eval->ComputeAtNewPosition(new_x, tspec, new_x_ordering);
2030  tspec_sav = tspec;
2031 
2032  good_tspec = reuse_flag;
2033 }
2034 
2035 void DiscreteAdaptTC::UpdateTargetSpecification(Vector &new_x,
2036  Vector &IntData,
2037  int new_x_ordering)
2038 {
2039  adapt_eval->ComputeAtNewPosition(new_x, IntData, new_x_ordering);
2040 }
2041 
2042 void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
2043  ElementTransformation &T,
2044  int dofidx, int dir,
2045  const Vector &IntData)
2046 {
2047  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
2048 
2049  Array<int> dofs;
2050  tspec_fesv->GetElementDofs(T.ElementNo, dofs);
2051  const int cnt = tspec.Size()/ncomp; // dofs per scalar-field
2052 
2053  for (int i = 0; i < ncomp; i++)
2054  {
2055  tspec(dofs[dofidx]+i*cnt) = IntData(dofs[dofidx] + i*cnt + dir*cnt*ncomp);
2056  }
2057 }
2058 
2059 void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
2060  int dofidx)
2061 {
2062  MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
2063 
2064  Array<int> dofs;
2065  tspec_fesv->GetElementDofs(T.ElementNo, dofs);
2066  const int cnt = tspec.Size()/ncomp;
2067  for (int i = 0; i < ncomp; i++)
2068  {
2069  tspec(dofs[dofidx] + i*cnt) = tspec_sav(dofs[dofidx] + i*cnt);
2070  }
2071 }
2072 
2073 void DiscreteAdaptTC::SetTspecFromIntRule(int e_id,
2074  const IntegrationRule &intrule)
2075 {
2076  switch (target_type)
2077  {
2078  case IDEAL_SHAPE_GIVEN_SIZE:
2079  case GIVEN_SHAPE_AND_SIZE:
2080  {
2081  const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
2082  ntspec_dofs = ndofs*ncomp;
2083 
2084  Vector tspec_vals(ntspec_dofs);
2085 
2086  Array<int> dofs;
2087  tspec_fesv->GetElementVDofs(e_id, dofs);
2088  tspec.GetSubVector(dofs, tspec_vals);
2089  DenseMatrix tr;
2090  tspec_gf->GetVectorValues(e_id, intrule, tspec_refine, tr);
2091  tspec_refine.Transpose();
2092  break;
2093  }
2094  default:
2095  MFEM_ABORT("Incompatible target type for discrete adaptation!");
2096  }
2097 }
2098 
2099 void DiscreteAdaptTC::SetTspecDataForDerefinement(FiniteElementSpace *fes)
2100 {
2101  coarse_tspec_fesv = fes;
2102  const Operator *c_op = fes->GetUpdateOperator();
2103  tspec_derefine.SetSize(c_op->Height());
2104  c_op->Mult(tspec, tspec_derefine);
2105 }
2106 
2107 void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
2108  const IntegrationRule &ir,
2109  const Vector &elfun,
2110  DenseTensor &Jtr) const
2111 {
2112  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
2113  const int dim = fe.GetDim(),
2114  nqp = ir.GetNPoints();
2115  Jtrcomp.SetSize(dim, dim, 4*nqp);
2116 
2117  FiniteElementSpace *src_fes = tspec_fesv;
2118 
2119  switch (target_type)
2120  {
2121  case IDEAL_SHAPE_GIVEN_SIZE:
2122  case GIVEN_SHAPE_AND_SIZE:
2123  {
2124  const DenseMatrix &Wideal =
2125  Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
2126  const int ndofs = tspec_fesv->GetFE(e_id)->GetDof(),
2127  ntspec_dofs = ndofs*ncomp;
2128 
2129  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
2130  par_vals_c1, par_vals_c2, par_vals_c3;
2131 
2132  Array<int> dofs;
2133  DenseMatrix D_rho(dim), Q_phi(dim), R_theta(dim);
2134  tspec_fesv->GetElementVDofs(e_id, dofs);
2135  tspec.UseDevice(true);
2136  tspec.GetSubVector(dofs, tspec_vals);
2137  if (tspec_refine.NumCols() > 0) // Refinement
2138  {
2139  MFEM_VERIFY(amr_el >= 0, " Target being constructed for an AMR element.");
2140  for (int i = 0; i < ncomp; i++)
2141  {
2142  for (int j = 0; j < ndofs; j++)
2143  {
2144  tspec_vals(j + i*ndofs) = tspec_refine(j + amr_el*ndofs, i);
2145  }
2146  }
2147  }
2148  else if (tspec_derefine.Size() > 0) // Derefinement
2149  {
2150  dofs.SetSize(0);
2151  coarse_tspec_fesv->GetElementVDofs(e_id, dofs);
2152  tspec_derefine.GetSubVector(dofs, tspec_vals);
2153  src_fes = coarse_tspec_fesv;
2154  }
2155 
2156  for (int q = 0; q < nqp; q++)
2157  {
2158  const IntegrationPoint &ip = ir.IntPoint(q);
2159  src_fes->GetFE(e_id)->CalcShape(ip, shape);
2160  Jtr(q) = Wideal; // Initialize to identity
2161  for (int d = 0; d < 4; d++)
2162  {
2163  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(d + 4*q), dim, dim);
2164  Jtrcomp_q = Wideal; // Initialize to identity
2165  }
2166 
2167  if (sizeidx != -1) // Set size
2168  {
2169  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
2170  double min_size = par_vals.Min();
2171  if (lim_min_size > 0.) { min_size = lim_min_size; }
2172  MFEM_VERIFY(min_size > 0.0,
2173  "Non-positive size propagated in the target definition.");
2174 
2175  double size = fmax(shape * par_vals, min_size);
2176  NCMesh *ncmesh = tspec_fesv->GetMesh()->ncmesh;
2177  if (ncmesh)
2178  {
2179  size /= ncmesh->GetElementSizeReduction(e_id);
2180  }
2181  Jtr(q).Set(std::pow(size, 1.0/dim), Jtr(q));
2182  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(0 + 4*q), dim, dim);
2183  Jtrcomp_q = Jtr(q);
2184  } // Done size
2185 
2186  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
2187 
2188  if (aspectratioidx != -1) // Set aspect ratio
2189  {
2190  if (dim == 2)
2191  {
2192  par_vals.SetDataAndSize(tspec_vals.GetData()+
2193  aspectratioidx*ndofs, ndofs);
2194  const double min_size = par_vals.Min();
2195  MFEM_VERIFY(min_size > 0.0,
2196  "Non-positive aspect-ratio propagated in the target definition.");
2197 
2198  const double aspectratio = shape * par_vals;
2199  D_rho = 0.;
2200  D_rho(0,0) = 1./pow(aspectratio,0.5);
2201  D_rho(1,1) = pow(aspectratio,0.5);
2202  }
2203  else
2204  {
2205  par_vals.SetDataAndSize(tspec_vals.GetData()+
2206  aspectratioidx*ndofs, ndofs*3);
2207  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
2208  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
2209  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
2210 
2211  const double rho1 = shape * par_vals_c1;
2212  const double rho2 = shape * par_vals_c2;
2213  const double rho3 = shape * par_vals_c3;
2214  D_rho = 0.;
2215  D_rho(0,0) = pow(rho1,2./3.);
2216  D_rho(1,1) = pow(rho2,2./3.);
2217  D_rho(2,2) = pow(rho3,2./3.);
2218  }
2219  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(1 + 4*q), dim, dim);
2220  Jtrcomp_q = D_rho;
2221  DenseMatrix Temp = Jtr(q);
2222  Mult(D_rho, Temp, Jtr(q));
2223  } // Done aspect ratio
2224 
2225  if (skewidx != -1) // Set skew
2226  {
2227  if (dim == 2)
2228  {
2229  par_vals.SetDataAndSize(tspec_vals.GetData()+
2230  skewidx*ndofs, ndofs);
2231 
2232  const double skew = shape * par_vals;
2233 
2234  Q_phi = 0.;
2235  Q_phi(0,0) = 1.;
2236  Q_phi(0,1) = cos(skew);
2237  Q_phi(1,1) = sin(skew);
2238  }
2239  else
2240  {
2241  par_vals.SetDataAndSize(tspec_vals.GetData()+
2242  skewidx*ndofs, ndofs*3);
2243  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
2244  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
2245  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
2246 
2247  const double phi12 = shape * par_vals_c1;
2248  const double phi13 = shape * par_vals_c2;
2249  const double chi = shape * par_vals_c3;
2250 
2251  Q_phi = 0.;
2252  Q_phi(0,0) = 1.;
2253  Q_phi(0,1) = cos(phi12);
2254  Q_phi(0,2) = cos(phi13);
2255 
2256  Q_phi(1,1) = sin(phi12);
2257  Q_phi(1,2) = sin(phi13)*cos(chi);
2258 
2259  Q_phi(2,2) = sin(phi13)*sin(chi);
2260  }
2261  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*q), dim, dim);
2262  Jtrcomp_q = Q_phi;
2263  DenseMatrix Temp = Jtr(q);
2264  Mult(Q_phi, Temp, Jtr(q));
2265  } // Done skew
2266 
2267  if (orientationidx != -1) // Set orientation
2268  {
2269  if (dim == 2)
2270  {
2271  par_vals.SetDataAndSize(tspec_vals.GetData()+
2272  orientationidx*ndofs, ndofs);
2273 
2274  const double theta = shape * par_vals;
2275  R_theta(0,0) = cos(theta);
2276  R_theta(0,1) = -sin(theta);
2277  R_theta(1,0) = sin(theta);
2278  R_theta(1,1) = cos(theta);
2279  }
2280  else
2281  {
2282  par_vals.SetDataAndSize(tspec_vals.GetData()+
2283  orientationidx*ndofs, ndofs*3);
2284  par_vals_c1.SetDataAndSize(par_vals.GetData(), ndofs);
2285  par_vals_c2.SetDataAndSize(par_vals.GetData()+ndofs, ndofs);
2286  par_vals_c3.SetDataAndSize(par_vals.GetData()+2*ndofs, ndofs);
2287 
2288  const double theta = shape * par_vals_c1;
2289  const double psi = shape * par_vals_c2;
2290  const double beta = shape * par_vals_c3;
2291 
2292  double ct = cos(theta), st = sin(theta),
2293  cp = cos(psi), sp = sin(psi),
2294  cb = cos(beta), sb = sin(beta);
2295 
2296  R_theta = 0.;
2297  R_theta(0,0) = ct*sp;
2298  R_theta(1,0) = st*sp;
2299  R_theta(2,0) = cp;
2300 
2301  R_theta(0,1) = -st*cb + ct*cp*sb;
2302  R_theta(1,1) = ct*cb + st*cp*sb;
2303  R_theta(2,1) = -sp*sb;
2304 
2305  R_theta(0,0) = -st*sb - ct*cp*cb;
2306  R_theta(1,0) = ct*sb - st*cp*cb;
2307  R_theta(2,0) = sp*cb;
2308  }
2309  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(3 + 4*q), dim, dim);
2310  Jtrcomp_q = R_theta;
2311  DenseMatrix Temp = Jtr(q);
2312  Mult(R_theta, Temp, Jtr(q));
2313  } // Done orientation
2314  }
2315  break;
2316  }
2317  default:
2318  MFEM_ABORT("Incompatible target type for discrete adaptation!");
2319  }
2320 }
2321 
2322 void DiscreteAdaptTC::ComputeElementTargetsGradient(const IntegrationRule &ir,
2323  const Vector &elfun,
2324  IsoparametricTransformation &Tpr,
2325  DenseTensor &dJtr) const
2326 {
2327  MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
2328 
2329  MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
2330 
2331  dJtr = 0.;
2332  const int e_id = Tpr.ElementNo;
2333  const FiniteElement *fe = Tpr.GetFE();
2334 
2335  switch (target_type)
2336  {
2337  case IDEAL_SHAPE_GIVEN_SIZE:
2338  case GIVEN_SHAPE_AND_SIZE:
2339  {
2340  const DenseMatrix &Wideal =
2341  Geometries.GetGeomToPerfGeomJac(fe->GetGeomType());
2342  const int dim = Wideal.Height(),
2343  ndofs = fe->GetDof(),
2344  ntspec_dofs = ndofs*ncomp;
2345 
2346  Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
2347  par_vals_c1(ndofs), par_vals_c2(ndofs), par_vals_c3(ndofs);
2348 
2349  Array<int> dofs;
2350  DenseMatrix dD_rho(dim), dQ_phi(dim), dR_theta(dim);
2351  DenseMatrix dQ_phi13(dim), dQ_phichi(dim); // dQ_phi is used for dQ/dphi12 in 3D
2352  DenseMatrix dR_psi(dim), dR_beta(dim);
2353  tspec_fesv->GetElementVDofs(e_id, dofs);
2354  tspec.GetSubVector(dofs, tspec_vals);
2355 
2356  DenseMatrix grad_e_c1(ndofs, dim),
2357  grad_e_c2(ndofs, dim),
2358  grad_e_c3(ndofs, dim);
2359  Vector grad_ptr_c1(grad_e_c1.GetData(), ndofs*dim),
2360  grad_ptr_c2(grad_e_c2.GetData(), ndofs*dim),
2361  grad_ptr_c3(grad_e_c3.GetData(), ndofs*dim);
2362 
2363  DenseMatrix grad_phys; // This will be (dof x dim, dof).
2364  fe->ProjectGrad(*fe, Tpr, grad_phys);
2365 
2366  for (int i = 0; i < ir.GetNPoints(); i++)
2367  {
2368  const IntegrationPoint &ip = ir.IntPoint(i);
2369  DenseMatrix Jtrcomp_s(Jtrcomp.GetData(0 + 4*i), dim, dim); // size
2370  DenseMatrix Jtrcomp_d(Jtrcomp.GetData(1 + 4*i), dim, dim); // aspect-ratio
2371  DenseMatrix Jtrcomp_q(Jtrcomp.GetData(2 + 4*i), dim, dim); // skew
2372  DenseMatrix Jtrcomp_r(Jtrcomp.GetData(3 + 4*i), dim, dim); // orientation
2373  DenseMatrix work1(dim), work2(dim), work3(dim);
2374 
2375  if (sizeidx != -1) // Set size
2376  {
2377  par_vals.SetDataAndSize(tspec_vals.GetData()+sizeidx*ndofs, ndofs);
2378 
2379  grad_phys.Mult(par_vals, grad_ptr_c1);
2380  Vector grad_q(dim);
2381  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2382  grad_e_c1.MultTranspose(shape, grad_q);
2383 
2384  const double min_size = par_vals.Min();
2385  MFEM_VERIFY(min_size > 0.0,
2386  "Non-positive size propagated in the target definition.");
2387  const double size = std::max(shape * par_vals, min_size);
2388  double dz_dsize = (1./dim)*pow(size, 1./dim - 1.);
2389 
2390  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2391  Mult(Jtrcomp_r, work1, work2); // R*Q*D
2392 
2393  for (int d = 0; d < dim; d++)
2394  {
2395  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2396  work1 = Wideal;
2397  work1.Set(dz_dsize, work1); // dz/dsize
2398  work1 *= grad_q(d); // dz/dsize*dsize/dx
2399  AddMult(work1, work2, dJtr_i); // dz/dx*R*Q*D
2400  }
2401  } // Done size
2402 
2403  if (target_type == IDEAL_SHAPE_GIVEN_SIZE) { continue; }
2404 
2405  if (aspectratioidx != -1) // Set aspect ratio
2406  {
2407  if (dim == 2)
2408  {
2409  par_vals.SetDataAndSize(tspec_vals.GetData()+
2410  aspectratioidx*ndofs, ndofs);
2411 
2412  grad_phys.Mult(par_vals, grad_ptr_c1);
2413  Vector grad_q(dim);
2414  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2415  grad_e_c1.MultTranspose(shape, grad_q);
2416 
2417  const double aspectratio = shape * par_vals;
2418  dD_rho = 0.;
2419  dD_rho(0,0) = -0.5*pow(aspectratio,-1.5);
2420  dD_rho(1,1) = 0.5*pow(aspectratio,-0.5);
2421 
2422  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
2423  Mult(work1, Jtrcomp_q, work2); // z*R*Q
2424 
2425  for (int d = 0; d < dim; d++)
2426  {
2427  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2428  work1 = dD_rho;
2429  work1 *= grad_q(d); // work1 = dD/drho*drho/dx
2430  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
2431  }
2432  }
2433  else // 3D
2434  {
2435  par_vals.SetDataAndSize(tspec_vals.GetData()+
2436  aspectratioidx*ndofs, ndofs*3);
2437  par_vals_c1.SetData(par_vals.GetData());
2438  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2439  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2440 
2441  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2442  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2443  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2444  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2445  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2446  grad_e_c1.MultTranspose(shape, grad_q1);
2447  grad_e_c2.MultTranspose(shape, grad_q2);
2448  grad_e_c3.MultTranspose(shape, grad_q3);
2449 
2450  const double rho1 = shape * par_vals_c1;
2451  const double rho2 = shape * par_vals_c2;
2452  const double rho3 = shape * par_vals_c3;
2453  dD_rho = 0.;
2454  dD_rho(0,0) = (2./3.)*pow(rho1,-1./3.);
2455  dD_rho(1,1) = (2./3.)*pow(rho2,-1./3.);
2456  dD_rho(2,2) = (2./3.)*pow(rho3,-1./3.);
2457 
2458  Mult(Jtrcomp_s, Jtrcomp_r, work1); // z*R
2459  Mult(work1, Jtrcomp_q, work2); // z*R*Q
2460 
2461 
2462  for (int d = 0; d < dim; d++)
2463  {
2464  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2465  work1 = dD_rho;
2466  work1(0,0) *= grad_q1(d);
2467  work1(1,2) *= grad_q2(d);
2468  work1(2,2) *= grad_q3(d);
2469  // work1 = dD/dx = dD/drho1*drho1/dx + dD/drho2*drho2/dx
2470  AddMult(work2, work1, dJtr_i); // z*R*Q*dD/dx
2471  }
2472  }
2473  } // Done aspect ratio
2474 
2475  if (skewidx != -1) // Set skew
2476  {
2477  if (dim == 2)
2478  {
2479  par_vals.SetDataAndSize(tspec_vals.GetData()+
2480  skewidx*ndofs, ndofs);
2481 
2482  grad_phys.Mult(par_vals, grad_ptr_c1);
2483  Vector grad_q(dim);
2484  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2485  grad_e_c1.MultTranspose(shape, grad_q);
2486 
2487  const double skew = shape * par_vals;
2488 
2489  dQ_phi = 0.;
2490  dQ_phi(0,0) = 1.;
2491  dQ_phi(0,1) = -sin(skew);
2492  dQ_phi(1,1) = cos(skew);
2493 
2494  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
2495 
2496  for (int d = 0; d < dim; d++)
2497  {
2498  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2499  work1 = dQ_phi;
2500  work1 *= grad_q(d); // work1 = dQ/dphi*dphi/dx
2501  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
2502  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
2503  }
2504  }
2505  else
2506  {
2507  par_vals.SetDataAndSize(tspec_vals.GetData()+
2508  skewidx*ndofs, ndofs*3);
2509  par_vals_c1.SetData(par_vals.GetData());
2510  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2511  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2512 
2513  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2514  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2515  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2516  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2517  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2518  grad_e_c1.MultTranspose(shape, grad_q1);
2519  grad_e_c2.MultTranspose(shape, grad_q2);
2520  grad_e_c3.MultTranspose(shape, grad_q3);
2521 
2522  const double phi12 = shape * par_vals_c1;
2523  const double phi13 = shape * par_vals_c2;
2524  const double chi = shape * par_vals_c3;
2525 
2526  dQ_phi = 0.;
2527  dQ_phi(0,0) = 1.;
2528  dQ_phi(0,1) = -sin(phi12);
2529  dQ_phi(1,1) = cos(phi12);
2530 
2531  dQ_phi13 = 0.;
2532  dQ_phi13(0,2) = -sin(phi13);
2533  dQ_phi13(1,2) = cos(phi13)*cos(chi);
2534  dQ_phi13(2,2) = cos(phi13)*sin(chi);
2535 
2536  dQ_phichi = 0.;
2537  dQ_phichi(1,2) = -sin(phi13)*sin(chi);
2538  dQ_phichi(2,2) = sin(phi13)*cos(chi);
2539 
2540  Mult(Jtrcomp_s, Jtrcomp_r, work2); // z*R
2541 
2542  for (int d = 0; d < dim; d++)
2543  {
2544  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2545  work1 = dQ_phi;
2546  work1 *= grad_q1(d); // work1 = dQ/dphi12*dphi12/dx
2547  work1.Add(grad_q2(d), dQ_phi13); // + dQ/dphi13*dphi13/dx
2548  work1.Add(grad_q3(d), dQ_phichi); // + dQ/dchi*dchi/dx
2549  Mult(work1, Jtrcomp_d, work3); // dQ/dx*D
2550  AddMult(work2, work3, dJtr_i); // z*R*dQ/dx*D
2551  }
2552  }
2553  } // Done skew
2554 
2555  if (orientationidx != -1) // Set orientation
2556  {
2557  if (dim == 2)
2558  {
2559  par_vals.SetDataAndSize(tspec_vals.GetData()+
2560  orientationidx*ndofs, ndofs);
2561 
2562  grad_phys.Mult(par_vals, grad_ptr_c1);
2563  Vector grad_q(dim);
2564  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2565  grad_e_c1.MultTranspose(shape, grad_q);
2566 
2567  const double theta = shape * par_vals;
2568  dR_theta(0,0) = -sin(theta);
2569  dR_theta(0,1) = -cos(theta);
2570  dR_theta(1,0) = cos(theta);
2571  dR_theta(1,1) = -sin(theta);
2572 
2573  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2574  Mult(Jtrcomp_s, work1, work2); // z*Q*D
2575  for (int d = 0; d < dim; d++)
2576  {
2577  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2578  work1 = dR_theta;
2579  work1 *= grad_q(d); // work1 = dR/dtheta*dtheta/dx
2580  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
2581  }
2582  }
2583  else
2584  {
2585  par_vals.SetDataAndSize(tspec_vals.GetData()+
2586  orientationidx*ndofs, ndofs*3);
2587  par_vals_c1.SetData(par_vals.GetData());
2588  par_vals_c2.SetData(par_vals.GetData()+ndofs);
2589  par_vals_c3.SetData(par_vals.GetData()+2*ndofs);
2590 
2591  grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2592  grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2593  grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2594  Vector grad_q1(dim), grad_q2(dim), grad_q3(dim);
2595  tspec_fesv->GetFE(e_id)->CalcShape(ip, shape);
2596  grad_e_c1.MultTranspose(shape, grad_q1);
2597  grad_e_c2.MultTranspose(shape, grad_q2);
2598  grad_e_c3.MultTranspose(shape, grad_q3);
2599 
2600  const double theta = shape * par_vals_c1;
2601  const double psi = shape * par_vals_c2;
2602  const double beta = shape * par_vals_c3;
2603 
2604  const double ct = cos(theta), st = sin(theta),
2605  cp = cos(psi), sp = sin(psi),
2606  cb = cos(beta), sb = sin(beta);
2607 
2608  dR_theta = 0.;
2609  dR_theta(0,0) = -st*sp;
2610  dR_theta(1,0) = ct*sp;
2611  dR_theta(2,0) = 0;
2612 
2613  dR_theta(0,1) = -ct*cb - st*cp*sb;
2614  dR_theta(1,1) = -st*cb + ct*cp*sb;
2615  dR_theta(2,1) = 0.;
2616 
2617  dR_theta(0,0) = -ct*sb + st*cp*cb;
2618  dR_theta(1,0) = -st*sb - ct*cp*cb;
2619  dR_theta(2,0) = 0.;
2620 
2621  dR_beta = 0.;
2622  dR_beta(0,0) = 0.;
2623  dR_beta(1,0) = 0.;
2624  dR_beta(2,0) = 0.;
2625 
2626  dR_beta(0,1) = st*sb + ct*cp*cb;
2627  dR_beta(1,1) = -ct*sb + st*cp*cb;
2628  dR_beta(2,1) = -sp*cb;
2629 
2630  dR_beta(0,0) = -st*cb + ct*cp*sb;
2631  dR_beta(1,0) = ct*cb + st*cp*sb;
2632  dR_beta(2,0) = 0.;
2633 
2634  dR_psi = 0.;
2635  dR_psi(0,0) = ct*cp;
2636  dR_psi(1,0) = st*cp;
2637  dR_psi(2,0) = -sp;
2638 
2639  dR_psi(0,1) = 0. - ct*sp*sb;
2640  dR_psi(1,1) = 0. + st*sp*sb;
2641  dR_psi(2,1) = -cp*sb;
2642 
2643  dR_psi(0,0) = 0. + ct*sp*cb;
2644  dR_psi(1,0) = 0. + st*sp*cb;
2645  dR_psi(2,0) = cp*cb;
2646 
2647  Mult(Jtrcomp_q, Jtrcomp_d, work1); // Q*D
2648  Mult(Jtrcomp_s, work1, work2); // z*Q*D
2649  for (int d = 0; d < dim; d++)
2650  {
2651  DenseMatrix &dJtr_i = dJtr(i + d*ir.GetNPoints());
2652  work1 = dR_theta;
2653  work1 *= grad_q1(d); // work1 = dR/dtheta*dtheta/dx
2654  work1.Add(grad_q2(d), dR_psi); // +dR/dpsi*dpsi/dx
2655  work1.Add(grad_q3(d), dR_beta); // +dR/dbeta*dbeta/dx
2656  AddMult(work1, work2, dJtr_i); // z*dR/dx*Q*D
2657  }
2658  }
2659  } // Done orientation
2660  }
2661  break;
2662  }
2663  default:
2664  MFEM_ABORT("Incompatible target type for discrete adaptation!");
2665  }
2666  Jtrcomp.Clear();
2667 }
2668 
2669 void DiscreteAdaptTC:: UpdateGradientTargetSpecification(const Vector &x,
2670  const double dx,
2671  bool reuse_flag,
2672  int x_ordering)
2673 {
2674  if (reuse_flag && good_tspec_grad) { return; }
2675 
2676  const int dim = tspec_fesv->GetFE(0)->GetDim(),
2677  cnt = x.Size()/dim;
2678 
2679  tspec_pert1h.SetSize(x.Size()*ncomp);
2680 
2681  Vector TSpecTemp;
2682  Vector xtemp = x;
2683  for (int j = 0; j < dim; j++)
2684  {
2685  for (int i = 0; i < cnt; i++)
2686  {
2687  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2688  xtemp(idx) += dx;
2689  }
2690 
2691  TSpecTemp.NewDataAndSize(tspec_pert1h.GetData() + j*cnt*ncomp, cnt*ncomp);
2692  UpdateTargetSpecification(xtemp, TSpecTemp, x_ordering);
2693 
2694  for (int i = 0; i < cnt; i++)
2695  {
2696  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2697  xtemp(idx) -= dx;
2698  }
2699  }
2700 
2701  good_tspec_grad = reuse_flag;
2702 }
2703 
2704 void DiscreteAdaptTC::
2705 UpdateHessianTargetSpecification(const Vector &x,double dx,
2706  bool reuse_flag, int x_ordering)
2707 {
2708 
2709  if (reuse_flag && good_tspec_hess) { return; }
2710 
2711  const int dim = tspec_fesv->GetFE(0)->GetDim(),
2712  cnt = x.Size()/dim,
2713  totmix = 1+2*(dim-2);
2714 
2715  tspec_pert2h.SetSize(cnt*dim*ncomp);
2716  tspec_pertmix.SetSize(cnt*totmix*ncomp);
2717 
2718  Vector TSpecTemp;
2719  Vector xtemp = x;
2720 
2721  // T(x+2h)
2722  for (int j = 0; j < dim; j++)
2723  {
2724  for (int i = 0; i < cnt; i++)
2725  {
2726  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2727  xtemp(idx) += 2*dx;
2728  }
2729 
2730  TSpecTemp.NewDataAndSize(tspec_pert2h.GetData() + j*cnt*ncomp, cnt*ncomp);
2731  UpdateTargetSpecification(xtemp, TSpecTemp, x_ordering);
2732 
2733  for (int i = 0; i < cnt; i++)
2734  {
2735  int idx = x_ordering == Ordering::byNODES ? j*cnt + i : i*dim + j;
2736  xtemp(idx) -= 2*dx;
2737  }
2738  }
2739 
2740  // T(x+h,y+h)
2741  int j = 0;
2742  for (int k1 = 0; k1 < dim; k1++)
2743  {
2744  for (int k2 = 0; (k1 != k2) && (k2 < dim); k2++)
2745  {
2746  for (int i = 0; i < cnt; i++)
2747  {
2748  int idx1 = x_ordering == Ordering::byNODES ? k1*cnt+i : i*dim + k1;
2749  int idx2 = x_ordering == Ordering::byNODES ? k2*cnt+i : i*dim + k2;
2750  xtemp(idx1) += dx;
2751  xtemp(idx2) += dx;
2752  }
2753 
2754  TSpecTemp.NewDataAndSize(tspec_pertmix.GetData() + j*cnt*ncomp, cnt*ncomp);
2755  UpdateTargetSpecification(xtemp, TSpecTemp, x_ordering);
2756 
2757  for (int i = 0; i < cnt; i++)
2758  {
2759  int idx1 = x_ordering == Ordering::byNODES ? k1*cnt+i : i*dim + k1;
2760  int idx2 = x_ordering == Ordering::byNODES ? k2*cnt+i : i*dim + k2;
2761  xtemp(idx1) -= dx;
2762  xtemp(idx2) -= dx;
2763  }
2764  j++;
2765  }
2766  }
2767 
2768  good_tspec_hess = reuse_flag;
2769 }
2770 
2771 DiscreteAdaptTC::~DiscreteAdaptTC()
2772 {
2773  delete tspec_gf;
2774  delete adapt_eval;
2775  delete tspec_fesv;
2776 #ifdef MFEM_USE_MPI
2777  delete ptspec_fesv;
2778 #endif
2779 }
2780 
2781 void AdaptivityEvaluator::SetSerialMetaInfo(const Mesh &m,
2782  const FiniteElementSpace &f)
2783 {
2784  delete fes;
2785  delete mesh;
2786  mesh = new Mesh(m, true);
2787  fes = new FiniteElementSpace(mesh, f.FEColl(),
2788  f.GetVDim(), f.GetOrdering());
2789 }
2790 
2791 #ifdef MFEM_USE_MPI
2792 void AdaptivityEvaluator::SetParMetaInfo(const ParMesh &m,
2793  const ParFiniteElementSpace &f)
2794 {
2795  delete pfes;
2796  delete pmesh;
2797  pmesh = new ParMesh(m, true);
2798  pfes = new ParFiniteElementSpace(pmesh, f.FEColl(),
2799  f.GetVDim(), f.GetOrdering());
2800 }
2801 #endif
2802 
2803 void AdaptivityEvaluator::ClearGeometricFactors()
2804 {
2805 #ifdef MFEM_USE_MPI
2806  if (pmesh) { pmesh->DeleteGeometricFactors(); }
2807 #else
2808  if (mesh) { mesh->DeleteGeometricFactors(); }
2809 #endif
2810 }
2811 
2812 AdaptivityEvaluator::~AdaptivityEvaluator()
2813 {
2814  delete fes;
2815  delete mesh;
2816 #ifdef MFEM_USE_MPI
2817  delete pfes;
2818  delete pmesh;
2819 #endif
2820 }
2821 
2822 void TMOP_Integrator::ReleasePADeviceMemory(bool copy_to_host)
2823 {
2824  if (PA.enabled)
2825  {
2826  PA.H.GetMemory().DeleteDevice(copy_to_host);
2827  PA.H0.GetMemory().DeleteDevice(copy_to_host);
2828  if (!copy_to_host && !PA.Jtr.GetMemory().HostIsValid())
2829  {
2830  PA.Jtr_needs_update = true;
2831  }
2832  PA.Jtr.GetMemory().DeleteDevice(copy_to_host);
2833  }
2834 }
2835 
2836 TMOP_Integrator::~TMOP_Integrator()
2837 {
2838  delete lim_func;
2839  delete adapt_lim_gf;
2840  delete surf_fit_gf;
2841  delete surf_fit_limiter;
2842  delete surf_fit_grad;
2843  delete surf_fit_hess;
2844  for (int i = 0; i < ElemDer.Size(); i++)
2845  {
2846  delete ElemDer[i];
2847  delete ElemPertEnergy[i];
2848  }
2849 }
2850 
2851 void TMOP_Integrator::EnableLimiting(const GridFunction &n0,
2852  const GridFunction &dist, Coefficient &w0,
2853  TMOP_LimiterFunction *lfunc)
2854 {
2855  lim_nodes0 = &n0;
2856  lim_coeff = &w0;
2857  lim_dist = &dist;
2858  MFEM_VERIFY(lim_dist->FESpace()->GetVDim() == 1,
2859  "'dist' must be a scalar GridFunction!");
2860 
2861  delete lim_func;
2862  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2863 }
2864 
2865 void TMOP_Integrator::EnableLimiting(const GridFunction &n0, Coefficient &w0,
2866  TMOP_LimiterFunction *lfunc)
2867 {
2868  lim_nodes0 = &n0;
2869  lim_coeff = &w0;
2870  lim_dist = NULL;
2871 
2872  delete lim_func;
2873  lim_func = (lfunc) ? lfunc : new TMOP_QuadraticLimiter;
2874 }
2875 
2876 void TMOP_Integrator::EnableAdaptiveLimiting(const GridFunction &z0,
2877  Coefficient &coeff,
2878  AdaptivityEvaluator &ae)
2879 {
2880  adapt_lim_gf0 = &z0;
2881  delete adapt_lim_gf;
2882  adapt_lim_gf = new GridFunction(z0);
2883  adapt_lim_coeff = &coeff;
2884  adapt_lim_eval = &ae;
2885 
2886  adapt_lim_eval->SetSerialMetaInfo(*z0.FESpace()->GetMesh(),
2887  *z0.FESpace());
2888  adapt_lim_eval->SetInitialField
2889  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2890 }
2891 
2892 #ifdef MFEM_USE_MPI
2893 void TMOP_Integrator::EnableAdaptiveLimiting(const ParGridFunction &z0,
2894  Coefficient &coeff,
2895  AdaptivityEvaluator &ae)
2896 {
2897  adapt_lim_gf0 = &z0;
2898  adapt_lim_pgf0 = &z0;
2899  delete adapt_lim_gf;
2900  adapt_lim_gf = new GridFunction(z0);
2901  adapt_lim_coeff = &coeff;
2902  adapt_lim_eval = &ae;
2903 
2904  adapt_lim_eval->SetParMetaInfo(*z0.ParFESpace()->GetParMesh(),
2905  *z0.ParFESpace());
2906  adapt_lim_eval->SetInitialField
2907  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
2908 }
2909 #endif
2910 
2911 void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &s0,
2912  const Array<bool> &smarker,
2913  Coefficient &coeff,
2914  AdaptivityEvaluator &ae)
2915 {
2916  // To have both we must duplicate the markers.
2917  MFEM_VERIFY(surf_fit_pos == NULL,
2918  "Using both fitting approaches is not supported.");
2919 
2920  delete surf_fit_gf;
2921  surf_fit_gf = new GridFunction(s0);
2922  surf_fit_gf->CountElementsPerVDof(surf_fit_dof_count);
2923  surf_fit_marker = &smarker;
2924  surf_fit_coeff = &coeff;
2925  surf_fit_eval = &ae;
2926 
2927  surf_fit_eval->SetSerialMetaInfo(*s0.FESpace()->GetMesh(),
2928  *s0.FESpace());
2929  surf_fit_eval->SetInitialField
2930  (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
2931 }
2932 
2933 void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &pos,
2934  const Array<bool> &smarker,
2935  Coefficient &coeff)
2936 {
2937  // To have both we must duplicate the markers.
2938  MFEM_VERIFY(surf_fit_gf == NULL,
2939  "Using both fitting approaches is not supported.");
2940  MFEM_VERIFY(pos.FESpace()->GetMesh()->GetNodes(),
2941  "Positions on a mesh without Nodes is not supported.");
2942  MFEM_VERIFY(pos.FESpace()->GetOrdering() ==
2943  pos.FESpace()->GetMesh()->GetNodes()->FESpace()->GetOrdering(),
2944  "Incompatible ordering of spaces!");
2945 
2946  surf_fit_pos = &pos;
2947  pos.CountElementsPerVDof(surf_fit_dof_count);
2948  surf_fit_marker = &smarker;
2949  surf_fit_coeff = &coeff;
2950  delete surf_fit_limiter;
2951  surf_fit_limiter = new TMOP_QuadraticLimiter;
2952 }
2953 
2954 #ifdef MFEM_USE_MPI
2955 void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
2956  const Array<bool> &smarker,
2957  Coefficient &coeff,
2958  AdaptivityEvaluator &ae)
2959 {
2960  // To have both we must duplicate the markers.
2961  MFEM_VERIFY(surf_fit_pos == NULL,
2962  "Using both fitting approaches is not supported.");
2963 
2964  delete surf_fit_gf;
2965  surf_fit_gf = new GridFunction(s0);
2966  s0.CountElementsPerVDof(surf_fit_dof_count);
2967  surf_fit_marker = &smarker;
2968  surf_fit_coeff = &coeff;
2969  surf_fit_eval = &ae;
2970 
2971  surf_fit_eval->SetParMetaInfo(*s0.ParFESpace()->GetParMesh(),
2972  *s0.ParFESpace());
2973  surf_fit_eval->SetInitialField
2974  (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
2975  surf_fit_gf_bg = false;
2976 }
2977 
2978 void TMOP_Integrator::EnableSurfaceFittingFromSource(
2979  const ParGridFunction &s_bg, ParGridFunction &s0,
2980  const Array<bool> &smarker, Coefficient &coeff, AdaptivityEvaluator &ae,
2981  const ParGridFunction &s_bg_grad,
2982  ParGridFunction &s0_grad, AdaptivityEvaluator &age,
2983  const ParGridFunction &s_bg_hess,
2984  ParGridFunction &s0_hess, AdaptivityEvaluator &ahe)
2985 {
2986 #ifndef MFEM_USE_GSLIB
2987  MFEM_ABORT("Surface fitting from source requires GSLIB!");
2988 #endif
2989 
2990  // Setup for level set function
2991  delete surf_fit_gf;
2992  surf_fit_gf = new GridFunction(s0);
2993  *surf_fit_gf = 0.0;
2994  surf_fit_marker = &smarker;
2995  surf_fit_coeff = &coeff;
2996  surf_fit_eval = &ae;
2997 
2998  surf_fit_gf_bg = true;
2999  surf_fit_eval->SetParMetaInfo(*s_bg.ParFESpace()->GetParMesh(),
3000  *s_bg.ParFESpace());
3001  surf_fit_eval->SetInitialField
3002  (*s_bg.FESpace()->GetMesh()->GetNodes(), s_bg);
3003 
3004  // Setup for gradient on background mesh
3005  MFEM_VERIFY(s_bg_grad.ParFESpace()->GetOrdering() ==
3006  s0_grad.ParFESpace()->GetOrdering(),
3007  "Nodal ordering for gridfunction on source mesh and current mesh"
3008  "should be the same.");
3009  delete surf_fit_grad;
3010  surf_fit_grad = new GridFunction(s0_grad);
3011  *surf_fit_grad = 0.0;
3012  surf_fit_eval_bg_grad = &age;
3013  surf_fit_eval_bg_hess = &ahe;
3014  surf_fit_eval_bg_grad->SetParMetaInfo(*s_bg_grad.ParFESpace()->GetParMesh(),
3015  *s_bg_grad.ParFESpace());
3016  surf_fit_eval_bg_grad->SetInitialField
3017  (*s_bg_grad.FESpace()->GetMesh()->GetNodes(), s_bg_grad);
3018 
3019  // Setup for Hessian on background mesh
3020  MFEM_VERIFY(s_bg_hess.ParFESpace()->GetOrdering() ==
3021  s0_hess.ParFESpace()->GetOrdering(),
3022  "Nodal ordering for gridfunction on source mesh and current mesh"
3023  "should be the same.");
3024  delete surf_fit_hess;
3025  surf_fit_hess = new GridFunction(s0_hess);
3026  *surf_fit_hess = 0.0;
3027  surf_fit_eval_bg_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(),
3028  *s_bg_hess.ParFESpace());
3029  surf_fit_eval_bg_hess->SetInitialField
3030  (*s_bg_hess.FESpace()->GetMesh()->GetNodes(), s_bg_hess);
3031 
3032  // Count number of zones that share each of the DOFs
3033  s0.CountElementsPerVDof(surf_fit_dof_count);
3034  // Store DOF indices that are marked for fitting. Used to reduce work for
3035  // transferring information between source/background and current mesh.
3036  surf_fit_marker_dof_index.SetSize(0);
3037  for (int i = 0; i < surf_fit_marker->Size(); i++)
3038  {
3039  if ((*surf_fit_marker)[i] == true)
3040  {
3041  surf_fit_marker_dof_index.Append(i);
3042  }
3043  }
3044 }
3045 #endif
3046 
3047 void TMOP_Integrator::GetSurfaceFittingErrors(const Vector &pos,
3048  double &err_avg, double &err_max)
3049 {
3050  MFEM_VERIFY(surf_fit_marker, "Surface fitting has not been enabled.");
3051 
3052  const FiniteElementSpace *fes =
3053  (surf_fit_gf) ? surf_fit_gf->FESpace() : surf_fit_pos->FESpace();
3054 #ifdef MFEM_USE_MPI
3055  auto pfes = dynamic_cast<const ParFiniteElementSpace *>(fes);
3056  bool parallel = (pfes) ? true : false;
3057 #endif
3058 
3059  int dim = fes->GetMesh()->Dimension();
3060  const int node_cnt = surf_fit_marker->Size();
3061  err_max = 0.0;
3062  int dof_cnt = 0;
3063  double err_sum = 0.0;
3064  for (int i = 0; i < node_cnt; i++)
3065  {
3066  if ((*surf_fit_marker)[i] == false) { continue; }
3067 
3068 #ifdef MFEM_USE_MPI
3069  // Don't count the overlapping DOFs in parallel.
3070  // The pfes might be ordered byVDIM, while the loop goes consecutively.
3071  const int dof_i = pfes->DofToVDof(i, 0);
3072  if (parallel && pfes->GetLocalTDofNumber(dof_i) < 0) { continue; }
3073 #endif
3074 
3075  dof_cnt++;
3076  double sigma_s = 0.0;
3077  if (surf_fit_gf) { sigma_s = fabs((*surf_fit_gf)(i)); }
3078  if (surf_fit_pos)
3079  {
3080  Vector pos_s(dim), pos_s_target(dim);
3081  for (int d = 0; d < dim; d++)
3082  {
3083  pos_s(d) = (fes->GetOrdering() == Ordering::byNODES) ?
3084  pos(d*node_cnt + i) : pos(i*dim + d);
3085  pos_s_target(d) = (fes->GetOrdering() == Ordering::byNODES)
3086  ? (*surf_fit_pos)(d*node_cnt + i)
3087  : (*surf_fit_pos)(i*dim + d);
3088  }
3089  sigma_s = pos_s.DistanceTo(pos_s_target);
3090  }
3091 
3092  err_max = fmax(err_max, sigma_s);
3093  err_sum += sigma_s;
3094  }
3095 
3096 #ifdef MFEM_USE_MPI
3097  if (parallel)
3098  {
3099  MPI_Comm comm = pfes->GetComm();
3100  MPI_Allreduce(MPI_IN_PLACE, &err_max, 1, MPI_DOUBLE, MPI_MAX, comm);
3101  MPI_Allreduce(MPI_IN_PLACE, &dof_cnt, 1, MPI_INT, MPI_SUM, comm);
3102  MPI_Allreduce(MPI_IN_PLACE, &err_sum, 1, MPI_DOUBLE, MPI_SUM, comm);
3103  }
3104 #endif
3105 
3106  err_avg = (dof_cnt > 0) ? err_sum / dof_cnt : 0.0;
3107 }
3108 
3109 void TMOP_Integrator::UpdateAfterMeshTopologyChange()
3110 {
3111  if (adapt_lim_gf)
3112  {
3113  adapt_lim_gf->Update();
3114  adapt_lim_eval->SetSerialMetaInfo(*adapt_lim_gf->FESpace()->GetMesh(),
3115  *adapt_lim_gf->FESpace());
3116  adapt_lim_eval->SetInitialField
3117  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
3118  }
3119 }
3120 
3121 #ifdef MFEM_USE_MPI
3122 void TMOP_Integrator::ParUpdateAfterMeshTopologyChange()
3123 {
3124  if (adapt_lim_gf)
3125  {
3126  adapt_lim_gf->Update();
3127  adapt_lim_eval->SetParMetaInfo(*adapt_lim_pgf0->ParFESpace()->GetParMesh(),
3128  *adapt_lim_pgf0->ParFESpace());
3129  adapt_lim_eval->SetInitialField
3130  (*adapt_lim_gf->FESpace()->GetMesh()->GetNodes(), *adapt_lim_gf);
3131  }
3132 }
3133 #endif
3134 
3135 double TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
3136  ElementTransformation &T,
3137  const Vector &elfun)
3138 {
3139  const int dof = el.GetDof(), dim = el.GetDim();
3140  const int el_id = T.ElementNo;
3141  double energy;
3142 
3143  // No adaptive limiting / surface fitting terms if the function is called
3144  // as part of a FD derivative computation (because we include the exact
3145  // derivatives of these terms in FD computations).
3146  const bool adaptive_limiting = (adapt_lim_gf && fd_call_flag == false);
3147  const bool surface_fit = (surf_fit_marker && fd_call_flag == false);
3148 
3149  DSh.SetSize(dof, dim);
3150  Jrt.SetSize(dim);
3151  Jpr.SetSize(dim);
3152  Jpt.SetSize(dim);
3153  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3154 
3155  const IntegrationRule &ir = EnergyIntegrationRule(el);
3156 
3157  energy = 0.0;
3158  DenseTensor Jtr(dim, dim, ir.GetNPoints());
3159  targetC->ComputeElementTargets(el_id, el, ir, elfun, Jtr);
3160 
3161  // Limited case.
3162  Vector shape, p, p0, d_vals;
3163  DenseMatrix pos0;
3164  if (lim_coeff)
3165  {
3166  shape.SetSize(dof);
3167  p.SetSize(dim);
3168  p0.SetSize(dim);
3169  pos0.SetSize(dof, dim);
3170  Vector pos0V(pos0.Data(), dof * dim);
3171  Array<int> pos_dofs;
3172  lim_nodes0->FESpace()->GetElementVDofs(el_id, pos_dofs);
3173  lim_nodes0->GetSubVector(pos_dofs, pos0V);
3174  if (lim_dist)
3175  {
3176  lim_dist->GetValues(el_id, ir, d_vals);
3177  }
3178  else
3179  {
3180  d_vals.SetSize(ir.GetNPoints()); d_vals = 1.0;
3181  }
3182  }
3183 
3184  // Define ref->physical transformation, when a Coefficient is specified.
3185  IsoparametricTransformation *Tpr = NULL;
3186  if (metric_coeff || lim_coeff || adaptive_limiting || surface_fit)
3187  {
3188  Tpr = new IsoparametricTransformation;
3189  Tpr->SetFE(&el);
3190  Tpr->ElementNo = el_id;
3191  Tpr->ElementType = ElementTransformation::ELEMENT;
3192  Tpr->Attribute = T.Attribute;
3193  Tpr->mesh = T.mesh;
3194  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3195  }
3196  // TODO: computing the coefficients 'metric_coeff' and 'lim_coeff' in physical
3197  // coordinates means that, generally, the gradient and Hessian of the
3198  // TMOP_Integrator will depend on the derivatives of the coefficients.
3199  //
3200  // In some cases the coefficients are independent of any movement of
3201  // the physical coordinates (i.e. changes in 'elfun'), e.g. when the
3202  // coefficient is a ConstantCoefficient or a GridFunctionCoefficient.
3203 
3204  Vector adapt_lim_gf_q, adapt_lim_gf0_q;
3205  if (adaptive_limiting)
3206  {
3207  adapt_lim_gf->GetValues(el_id, ir, adapt_lim_gf_q);
3208  adapt_lim_gf0->GetValues(el_id, ir, adapt_lim_gf0_q);
3209  }
3210 
3211  for (int i = 0; i < ir.GetNPoints(); i++)
3212  {
3213  const IntegrationPoint &ip = ir.IntPoint(i);
3214 
3215  metric->SetTargetJacobian(Jtr(i));
3216  CalcInverse(Jtr(i), Jrt);
3217  const double weight =
3218  (integ_over_target) ? ip.weight * Jtr(i).Det() : ip.weight;
3219 
3220  el.CalcDShape(ip, DSh);
3221  MultAtB(PMatI, DSh, Jpr);
3222  Mult(Jpr, Jrt, Jpt);
3223 
3224  double val = metric_normal * metric->EvalW(Jpt);
3225  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
3226 
3227  if (lim_coeff)
3228  {
3229  el.CalcShape(ip, shape);
3230  PMatI.MultTranspose(shape, p);
3231  pos0.MultTranspose(shape, p0);
3232  val += lim_normal *
3233  lim_func->Eval(p, p0, d_vals(i)) *
3234  lim_coeff->Eval(*Tpr, ip);
3235  }
3236 
3237  // Contribution from the adaptive limiting term.
3238  if (adaptive_limiting)
3239  {
3240  const double diff = adapt_lim_gf_q(i) - adapt_lim_gf0_q(i);
3241  val += adapt_lim_coeff->Eval(*Tpr, ip) * lim_normal * diff * diff;
3242  }
3243 
3244  energy += weight * val;
3245  }
3246 
3247  // Contribution from the surface fitting term.
3248  if (surface_fit)
3249  {
3250  // Scalar for surf_fit_gf, vector for surf_fit_pos, but that's ok.
3251  const FiniteElementSpace *fes_fit =
3252  (surf_fit_gf) ? surf_fit_gf->FESpace() : surf_fit_pos->FESpace();
3253  const IntegrationRule *ir_s = &fes_fit->GetFE(el_id)->GetNodes();
3254  Array<int> vdofs;
3255  fes_fit->GetElementVDofs(el_id, vdofs);
3256 
3257  Vector sigma_e(dof);
3258  if (surf_fit_gf) { surf_fit_gf->GetSubVector(vdofs, sigma_e); }
3259 
3260  for (int s = 0; s < dof; s++)
3261  {
3262  // Because surf_fit_pos.fes might be ordered byVDIM.
3263  const int scalar_dof_id = fes_fit->VDofToDof(vdofs[s]);
3264  if ((*surf_fit_marker)[scalar_dof_id] == false) { continue; }
3265 
3266  const IntegrationPoint &ip_s = ir_s->IntPoint(s);
3267  Tpr->SetIntPoint(&ip_s);
3268 
3269  if (surf_fit_gf)
3270  {
3271  energy += surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
3272  sigma_e(s) * sigma_e(s);
3273  }
3274  if (surf_fit_pos)
3275  {
3276  // Fitting to exact positions.
3277  Vector pos(dim), pos_target(dim);
3278  for (int d = 0; d < dim; d++)
3279  {
3280  pos(d) = PMatI(s, d);
3281  pos_target(d) = (*surf_fit_pos)(vdofs[d*dof + s]);
3282  }
3283  energy += surf_fit_coeff->Eval(*Tpr, ip_s) * surf_fit_normal *
3284  surf_fit_limiter->Eval(pos, pos_target, 1.0);
3285  }
3286  }
3287  }
3288 
3289  delete Tpr;
3290  return energy;
3291 }
3292 
3293 double TMOP_Integrator::GetRefinementElementEnergy(const FiniteElement &el,
3294  ElementTransformation &T,
3295  const Vector &elfun,
3296  const IntegrationRule &irule)
3297 {
3298  int dof = el.GetDof(), dim = el.GetDim(),
3299  NEsplit = elfun.Size() / (dof*dim), el_id = T.ElementNo;
3300  double energy = 0.;
3301 
3302  TargetConstructor *tc = const_cast<TargetConstructor *>(targetC);
3303  DiscreteAdaptTC *dtc = dynamic_cast<DiscreteAdaptTC *>(tc);
3304  // For DiscreteAdaptTC the GridFunctions used to set the targets must be
3305  // mapped onto the fine elements.
3306  if (dtc) { dtc->SetTspecFromIntRule(el_id, irule); }
3307 
3308  for (int e = 0; e < NEsplit; e++)
3309  {
3310  DSh.SetSize(dof, dim);
3311  Jrt.SetSize(dim);
3312  Jpr.SetSize(dim);
3313  Jpt.SetSize(dim);
3314  Vector elfun_child(dof*dim);
3315  for (int i = 0; i < dof; i++)
3316  {
3317  for (int d = 0; d < dim; d++)
3318  {
3319  // elfun is (xe1,xe2,...xen,ye1,ye2...yen) and has nodal coordinates
3320  // for all the children element of the parent element being considered.
3321  // So we must index and get (xek, yek) i.e. nodal coordinates for
3322  // the fine element being considered.
3323  elfun_child(i + d*dof) = elfun(i + e*dof + d*dof*NEsplit);
3324  }
3325  }
3326  PMatI.UseExternalData(elfun_child.GetData(), dof, dim);
3327 
3328  const IntegrationRule &ir = EnergyIntegrationRule(el);
3329 
3330  double el_energy = 0;
3331  DenseTensor Jtr(dim, dim, ir.GetNPoints());
3332  if (dtc)
3333  {
3334  // This is used to index into the tspec vector inside DiscreteAdaptTC.
3335  dtc->SetRefinementSubElement(e);
3336  }
3337  targetC->ComputeElementTargets(el_id, el, ir, elfun_child, Jtr);
3338 
3339  // Define ref->physical transformation, wn a Coefficient is specified.
3340  IsoparametricTransformation *Tpr = NULL;
3341  if (metric_coeff || lim_coeff)
3342  {
3343  Tpr = new IsoparametricTransformation;
3344  Tpr->SetFE(&el);
3345  Tpr->ElementNo = T.ElementNo;
3346  Tpr->ElementType = ElementTransformation::ELEMENT;
3347  Tpr->Attribute = T.Attribute;
3348  Tpr->mesh = T.mesh;
3349  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3350  }
3351 
3352  for (int i = 0; i < ir.GetNPoints(); i++)
3353  {
3354  const IntegrationPoint &ip = ir.IntPoint(i);
3355  h_metric->SetTargetJacobian(Jtr(i));
3356  CalcInverse(Jtr(i), Jrt);
3357  const double weight =
3358  (integ_over_target) ? ip.weight * Jtr(i).Det() : ip.weight;
3359 
3360  el.CalcDShape(ip, DSh);
3361  MultAtB(PMatI, DSh, Jpr);
3362  Mult(Jpr, Jrt, Jpt);
3363 
3364  double val = metric_normal * h_metric->EvalW(Jpt);
3365  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
3366 
3367  el_energy += weight * val;
3368  delete Tpr;
3369  }
3370  energy += el_energy;
3371  }
3372  energy /= NEsplit;
3373 
3374  if (dtc) { dtc->ResetRefinementTspecData(); }
3375 
3376  return energy;
3377 }
3378 
3379 double TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
3380  ElementTransformation &T,
3381  const Vector &elfun)
3382 {
3383  int dof = el.GetDof(), dim = el.GetDim();
3384  double energy = 0.;
3385 
3386  DSh.SetSize(dof, dim);
3387  Jrt.SetSize(dim);
3388  Jpr.SetSize(dim);
3389  Jpt.SetSize(dim);
3390  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3391 
3392  const IntegrationRule &ir = EnergyIntegrationRule(el);
3393 
3394  energy = 0.0;
3395  DenseTensor Jtr(dim, dim, ir.GetNPoints());
3396  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3397 
3398  // Define ref->physical transformation, wn a Coefficient is specified.
3399  IsoparametricTransformation *Tpr = NULL;
3400  if (metric_coeff)
3401  {
3402  Tpr = new IsoparametricTransformation;
3403  Tpr->SetFE(&el);
3404  Tpr->ElementNo = T.ElementNo;
3405  Tpr->ElementType = ElementTransformation::ELEMENT;
3406  Tpr->Attribute = T.Attribute;
3407  Tpr->mesh = T.mesh;
3408  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3409  }
3410 
3411  for (int i = 0; i < ir.GetNPoints(); i++)
3412  {
3413  const IntegrationPoint &ip = ir.IntPoint(i);
3414  h_metric->SetTargetJacobian(Jtr(i));
3415  CalcInverse(Jtr(i), Jrt);
3416  const double weight =
3417  (integ_over_target) ? ip.weight * Jtr(i).Det() : ip.weight;
3418 
3419  el.CalcDShape(ip, DSh);
3420  MultAtB(PMatI, DSh, Jpr);
3421  Mult(Jpr, Jrt, Jpt);
3422 
3423  double val = metric_normal * h_metric->EvalW(Jpt);
3424  if (metric_coeff) { val *= metric_coeff->Eval(*Tpr, ip); }
3425 
3426  energy += weight * val;
3427  }
3428 
3429  delete Tpr;
3430  return energy;
3431 }
3432 
3433 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
3434  ElementTransformation &T,
3435  const Vector &elfun, Vector &elvect)
3436 {
3437  if (!fdflag)
3438  {
3439  AssembleElementVectorExact(el, T, elfun, elvect);
3440  }
3441  else
3442  {
3443  AssembleElementVectorFD(el, T, elfun, elvect);
3444  }
3445 }
3446 
3447 void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,
3448  ElementTransformation &T,
3449  const Vector &elfun,
3450  DenseMatrix &elmat)
3451 {
3452  if (!fdflag)
3453  {
3454  AssembleElementGradExact(el, T, elfun, elmat);
3455  }
3456  else
3457  {
3458  AssembleElementGradFD(el, T, elfun, elmat);
3459  }
3460 }
3461 
3462 void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
3463  ElementTransformation &T,
3464  const Vector &elfun,
3465  Vector &elvect)
3466 {
3467  const int dof = el.GetDof(), dim = el.GetDim();
3468 
3469  DenseMatrix Amat(dim), work1(dim), work2(dim);
3470  DSh.SetSize(dof, dim);
3471  DS.SetSize(dof, dim);
3472  Jrt.SetSize(dim);
3473  Jpt.SetSize(dim);
3474  P.SetSize(dim);
3475  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3476  elvect.SetSize(dof*dim);
3477  PMatO.UseExternalData(elvect.GetData(), dof, dim);
3478 
3479  const IntegrationRule &ir = ActionIntegrationRule(el);
3480  const int nqp = ir.GetNPoints();
3481 
3482  elvect = 0.0;
3483  Vector weights(nqp);
3484  DenseTensor Jtr(dim, dim, nqp);
3485  DenseTensor dJtr(dim, dim, dim*nqp);
3486  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3487 
3488  // Limited case.
3489  DenseMatrix pos0;
3490  Vector shape, p, p0, d_vals, grad;
3491  shape.SetSize(dof);
3492  if (lim_coeff)
3493  {
3494  p.SetSize(dim);
3495  p0.SetSize(dim);
3496  pos0.SetSize(dof, dim);
3497  Vector pos0V(pos0.Data(), dof * dim);
3498  Array<int> pos_dofs;
3499  lim_nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
3500  lim_nodes0->GetSubVector(pos_dofs, pos0V);
3501  if (lim_dist)
3502  {
3503  lim_dist->GetValues(T.ElementNo, ir, d_vals);
3504  }
3505  else
3506  {
3507  d_vals.SetSize(nqp); d_vals = 1.0;
3508  }
3509  }
3510 
3511  // Define ref->physical transformation, when a Coefficient is specified.
3512  IsoparametricTransformation *Tpr = NULL;
3513  if (metric_coeff || lim_coeff || adapt_lim_gf ||
3514  surf_fit_gf || surf_fit_pos || exact_action)
3515  {
3516  Tpr = new IsoparametricTransformation;
3517  Tpr->SetFE(&el);
3518  Tpr->ElementNo = T.ElementNo;
3519  Tpr->ElementType = ElementTransformation::ELEMENT;
3520  Tpr->Attribute = T.Attribute;
3521  Tpr->mesh = T.mesh;
3522  Tpr->GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
3523  if (exact_action)
3524  {
3525  targetC->ComputeElementTargetsGradient(ir, elfun, *Tpr, dJtr);
3526  }
3527  }
3528 
3529  Vector d_detW_dx(dim);
3530  Vector d_Winv_dx(dim);
3531 
3532  for (int q = 0; q < nqp; q++)
3533  {
3534  const IntegrationPoint &ip = ir.IntPoint(q);
3535  metric->SetTargetJacobian(Jtr(q));
3536  CalcInverse(Jtr(q), Jrt);
3537  weights(q) = (integ_over_target) ? ip.weight * Jtr(q).Det() : ip.weight;
3538  double weight_m = weights(q) * metric_normal;
3539 
3540  el.CalcDShape(ip, DSh);
3541  Mult(DSh, Jrt, DS);
3542  MultAtB(PMatI, DS, Jpt);
3543 
3544  metric->EvalP(Jpt, P);
3545 
3546  if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
3547 
3548  P *= weight_m;
3549  AddMultABt(DS, P, PMatO); // w_q det(W) dmu/dx : dA/dx Winv
3550 
3551  if (exact_action)
3552  {
3553  el.CalcShape(ip, shape);
3554  // Derivatives of adaptivity-based targets.
3555  // First term: w_q d*(Det W)/dx * mu(T)
3556  // d(Det W)/dx = det(W)*Tr[Winv*dW/dx]
3557  DenseMatrix dwdx(dim);
3558  for (int d = 0; d < dim; d++)
3559  {
3560  const DenseMatrix &dJtr_q = dJtr(q + d * nqp);
3561  Mult(Jrt, dJtr_q, dwdx );
3562  d_detW_dx(d) = dwdx.Trace();
3563  }
3564  d_detW_dx *= weight_m*metric->EvalW(Jpt); // *[w_q*det(W)]*mu(T)
3565 
3566  // Second term: w_q det(W) dmu/dx : AdWinv/dx
3567  // dWinv/dx = -Winv*dW/dx*Winv
3568  MultAtB(PMatI, DSh, Amat);
3569  for (int d = 0; d < dim; d++)
3570  {
3571  const DenseMatrix &dJtr_q = dJtr(q + d*nqp);
3572  Mult(Jrt, dJtr_q, work1); // Winv*dw/dx
3573  Mult(work1, Jrt, work2); // Winv*dw/dx*Winv
3574  Mult(Amat, work2, work1); // A*Winv*dw/dx*Winv
3575  MultAtB(P, work1, work2); // dmu/dT^T*A*Winv*dw/dx*Winv
3576  d_Winv_dx(d) = work2.Trace(); // Tr[dmu/dT : AWinv*dw/dx*Winv]
3577  }
3578  d_Winv_dx *= -weight_m; // Include (-) factor as well
3579 
3580  d_detW_dx += d_Winv_dx;
3581  AddMultVWt(shape, d_detW_dx, PMatO);
3582  }
3583 
3584  if (lim_coeff)
3585  {
3586  if (!exact_action) { el.CalcShape(ip, shape); }
3587  PMatI.MultTranspose(shape, p);
3588  pos0.MultTranspose(shape, p0);
3589  lim_func->Eval_d1(p, p0, d_vals(q), grad);
3590  grad *= weights(q) * lim_normal * lim_coeff->Eval(*Tpr, ip);
3591  AddMultVWt(shape, grad, PMatO);
3592  }
3593  }
3594 
3595  if (adapt_lim_gf) { AssembleElemVecAdaptLim(el, *Tpr, ir, weights, PMatO); }
3596  if (surf_fit_gf || surf_fit_pos) { AssembleElemVecSurfFit(el, *Tpr, PMatO); }
3597 
3598  delete Tpr;
3599 }
3600 
3601 void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
3602  ElementTransformation &T,
3603  const Vector &elfun,
3604  DenseMatrix &elmat)
3605 {
3606  const int dof = el.GetDof(), dim = el.GetDim();
3607 
3608  DSh.SetSize(dof, dim);
3609  DS.SetSize(dof, dim);
3610  Jrt.SetSize(dim);
3611  Jpt.SetSize(dim);
3612  PMatI.UseExternalData(elfun.GetData(), dof, dim);
3613  elmat.SetSize(dof*dim);
3614 
3615  const IntegrationRule &ir = GradientIntegrationRule(el);
3616  const int nqp = ir.GetNPoints();
3617 
3618  elmat = 0.0;
3619  Vector weights(nqp);
3620  DenseTensor Jtr(dim, dim, nqp);
3621  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
3622 
3623  // Limited case.
3624  DenseMatrix pos0, hess;
3625  Vector shape, p, p0, d_vals;
3626  if (lim_coeff)
3627  {
3628  shape.SetSize(dof);
3629  p.SetSize(dim);
3630  p0.SetSize(dim);
3631  pos0.SetSize(dof, dim);
3632  Vector pos0V(pos0.Data(), dof * dim);
3633  Array<int> pos_dofs;
3634  lim_nodes0->FESpace()->GetElementVDofs(T.ElementNo, pos_dofs);
3635  lim_nodes0->GetSubVector(pos_dofs, pos0V);
3636  if (lim_dist)
3637  {
3638  lim_dist->GetValues(T.ElementNo, ir, d_vals);
3639  }
3640  else
3641  {
3642  d_vals.SetSize(nqp); d_vals = 1.0;
3643  }
3644  }
3645 
3646  // Define ref->physical transformation, when a Coefficient is specified.
3647  IsoparametricTransformation *Tpr = NULL;
3648  if (metric_coeff || lim_coeff || adapt_lim_gf || surf_fit_gf || surf_fit_pos)
3649  {
3650  Tpr = new IsoparametricTransformation;
3651  Tpr->SetFE(&el);
3652  Tpr->ElementNo = T.ElementNo;
3653  Tpr->ElementType = ElementTransformation::ELEMENT;
3654  Tpr->Attribute = T.Attribute;
3655  Tpr->mesh = T.mesh;
3656  Tpr->GetPointMat().Transpose(PMatI);
3657  }
3658 
3659  for (int q = 0; q < nqp; q++)
3660  {
3661  const IntegrationPoint &ip = ir.IntPoint(q);
3662  const DenseMatrix &Jtr_q = Jtr(q);
3663  metric->SetTargetJacobian(Jtr_q);
3664  CalcInverse(Jtr_q, Jrt);
3665  weights(q) = (integ_over_target) ? ip.weight * Jtr_q.Det() : ip.weight;
3666  double weight_m = weights(q) * metric_normal;
3667 
3668  el.CalcDShape(ip, DSh);
3669  Mult(DSh, Jrt, DS);
3670  MultAtB(PMatI, DS, Jpt);
3671 
3672  if (metric_coeff) { weight_m *= metric_coeff->Eval(*Tpr, ip); }
3673 
3674  metric->AssembleH(Jpt, DS, weight_m, elmat);
3675 
3676  // TODO: derivatives of adaptivity-based targets.
3677 
3678  if (lim_coeff)
3679  {
3680  el.CalcShape(ip, shape);
3681  PMatI.MultTranspose(shape, p);
3682  pos0.MultTranspose(shape, p0);
3683  weight_m = weights(q) * lim_normal * lim_coeff->Eval(*Tpr, ip);
3684  lim_func->Eval_d2(p, p0, d_vals(q), hess);
3685  for (int i = 0; i < dof; i++)
3686  {
3687  const double w_shape_i = weight_m * shape(i);
3688  for (int j = 0; j < dof; j++)
3689  {
3690  const double w = w_shape_i * shape(j);
3691  for (int d1 = 0; d1 < dim; d1++)
3692  {
3693  for (int d2 = 0; d2 < dim; d2++)
3694  {
3695  elmat(d1*dof + i, d2*dof + j) += w * hess(d1, d2);
3696  }
3697  }
3698  }
3699  }
3700  }
3701  }
3702 
3703  if (adapt_lim_gf) { AssembleElemGradAdaptLim(el, *Tpr, ir, weights, elmat); }
3704  if (surf_fit_gf || surf_fit_pos) { AssembleElemGradSurfFit(el, *Tpr, elmat);}
3705 
3706  delete Tpr;
3707 }
3708 
3709 void TMOP_Integrator::AssembleElemVecAdaptLim(const FiniteElement &el,
3710  IsoparametricTransformation &Tpr,
3711  const IntegrationRule &ir,
3712  const Vector &weights,
3713  DenseMatrix &mat)
3714 {
3715  const int dof = el.GetDof(), dim = el.GetDim(), nqp = weights.Size();
3716  Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3717 
3718  Array<int> dofs;
3719  adapt_lim_gf->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
3720  adapt_lim_gf->GetSubVector(dofs, adapt_lim_gf_e);
3721  adapt_lim_gf->GetValues(Tpr.ElementNo, ir, adapt_lim_gf_q);
3722  adapt_lim_gf0->GetValues(Tpr.ElementNo, ir, adapt_lim_gf0_q);
3723 
3724  // Project the gradient of adapt_lim_gf in the same space.
3725  // The FE coefficients of the gradient go in adapt_lim_gf_grad_e.
3726  DenseMatrix adapt_lim_gf_grad_e(dof, dim);
3727  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3728  el.ProjectGrad(el, Tpr, grad_phys);
3729  Vector grad_ptr(adapt_lim_gf_grad_e.GetData(), dof*dim);
3730  grad_phys.Mult(adapt_lim_gf_e, grad_ptr);
3731 
3732  Vector adapt_lim_gf_grad_q(dim);
3733 
3734  for (int q = 0; q < nqp; q++)
3735  {
3736  const IntegrationPoint &ip = ir.IntPoint(q);
3737  el.CalcShape(ip, shape);
3738  adapt_lim_gf_grad_e.MultTranspose(shape, adapt_lim_gf_grad_q);
3739  adapt_lim_gf_grad_q *= 2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q));
3740  adapt_lim_gf_grad_q *= weights(q) * lim_normal * adapt_lim_coeff->Eval(Tpr, ip);
3741  AddMultVWt(shape, adapt_lim_gf_grad_q, mat);
3742  }
3743 }
3744 
3745 void TMOP_Integrator::AssembleElemGradAdaptLim(const FiniteElement &el,
3746  IsoparametricTransformation &Tpr,
3747  const IntegrationRule &ir,
3748  const Vector &weights,
3749  DenseMatrix &mat)
3750 {
3751  const int dof = el.GetDof(), dim = el.GetDim(), nqp = weights.Size();
3752  Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3753 
3754  Array<int> dofs;
3755  adapt_lim_gf->FESpace()->GetElementDofs(Tpr.ElementNo, dofs);
3756  adapt_lim_gf->GetSubVector(dofs, adapt_lim_gf_e);
3757  adapt_lim_gf->GetValues(Tpr.ElementNo, ir, adapt_lim_gf_q);
3758  adapt_lim_gf0->GetValues(Tpr.ElementNo, ir, adapt_lim_gf0_q);
3759 
3760  // Project the gradient of adapt_lim_gf in the same space.
3761  // The FE coefficients of the gradient go in adapt_lim_gf_grad_e.
3762  DenseMatrix adapt_lim_gf_grad_e(dof, dim);
3763  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3764  el.ProjectGrad(el, Tpr, grad_phys);
3765  Vector grad_ptr(adapt_lim_gf_grad_e.GetData(), dof*dim);
3766  grad_phys.Mult(adapt_lim_gf_e, grad_ptr);
3767 
3768  // Project the gradient of each gradient of adapt_lim_gf in the same space.
3769  // The FE coefficients of the second derivatives go in adapt_lim_gf_hess_e.
3770  DenseMatrix adapt_lim_gf_hess_e(dof*dim, dim);
3771  Mult(grad_phys, adapt_lim_gf_grad_e, adapt_lim_gf_hess_e);
3772  // Reshape to be more convenient later (no change in the data).
3773  adapt_lim_gf_hess_e.SetSize(dof, dim*dim);
3774 
3775  Vector adapt_lim_gf_grad_q(dim);
3776  DenseMatrix adapt_lim_gf_hess_q(dim, dim);
3777 
3778  for (int q = 0; q < nqp; q++)
3779  {
3780  const IntegrationPoint &ip = ir.IntPoint(q);
3781  el.CalcShape(ip, shape);
3782 
3783  adapt_lim_gf_grad_e.MultTranspose(shape, adapt_lim_gf_grad_q);
3784  Vector gg_ptr(adapt_lim_gf_hess_q.GetData(), dim*dim);
3785  adapt_lim_gf_hess_e.MultTranspose(shape, gg_ptr);
3786 
3787  const double w = weights(q) * lim_normal * adapt_lim_coeff->Eval(Tpr, ip);
3788  for (int i = 0; i < dof * dim; i++)
3789  {
3790  const int idof = i % dof, idim = i / dof;
3791  for (int j = 0; j <= i; j++)
3792  {
3793  const int jdof = j % dof, jdim = j / dof;
3794  const double entry =
3795  w * ( 2.0 * adapt_lim_gf_grad_q(idim) * shape(idof) *
3796  /* */ adapt_lim_gf_grad_q(jdim) * shape(jdof) +
3797  2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q)) *
3798  adapt_lim_gf_hess_q(idim, jdim) * shape(idof) * shape(jdof));
3799  mat(i, j) += entry;
3800  if (i != j) { mat(j, i) += entry; }
3801  }
3802  }
3803  }
3804 }
3805 
3806 void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x,
3807  IsoparametricTransformation &Tpr,
3808  DenseMatrix &mat)
3809 {
3810  const int el_id = Tpr.ElementNo;
3811 
3812  // Scalar for surf_fit_gf, vector for surf_fit_pos, but that's ok.
3813  const FiniteElementSpace *fes_fit =
3814  (surf_fit_gf) ? surf_fit_gf->FESpace() : surf_fit_pos->FESpace();
3815  const FiniteElement &el_s = *fes_fit->GetFE(el_id);
3816  const int dof_s = el_s.GetDof(), dim = el_x.GetDim();
3817 
3818  // Check if the element has any DOFs marked for surface fitting.
3819  Array<int> dofs, vdofs;
3820  fes_fit->GetElementVDofs(el_id, vdofs);
3821  int count = 0;
3822  for (int s = 0; s < dof_s; s++)
3823  {
3824  // Because surf_fit_pos.fes might be ordered byVDIM.
3825  const int scalar_dof_id = fes_fit->VDofToDof(vdofs[s]);
3826  count += ((*surf_fit_marker)[scalar_dof_id]) ? 1 : 0;
3827  }
3828  if (count == 0) { return; }
3829 
3830  Vector sigma_e(dof_s);
3831  DenseMatrix surf_fit_grad_e(dof_s, dim);
3832  if (surf_fit_gf || surf_fit_gf_bg)
3833  {
3834  surf_fit_gf->GetSubVector(vdofs, sigma_e);
3835 
3836  // Project the gradient of sigma in the same space.
3837  // The FE coefficients of the gradient go in surf_fit_grad_e.
3838  Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
3839  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3840  if (surf_fit_gf_bg)
3841  {
3842  surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
3843  surf_fit_grad->GetSubVector(dofs, grad_ptr);
3844  }
3845  else
3846  {
3847  el_s.ProjectGrad(el_s, Tpr, grad_phys);
3848  grad_phys.Mult(sigma_e, grad_ptr);
3849  }
3850  }
3851  else { Tpr.GetPointMat().Transpose(PMatI); }
3852 
3853  const IntegrationRule &ir = el_s.GetNodes();
3854 
3855  for (int s = 0; s < dof_s; s++)
3856  {
3857  // Because surf_fit_pos.fes might be ordered byVDIM.
3858  const int scalar_dof_id = fes_fit->VDofToDof(vdofs[s]);
3859  if ((*surf_fit_marker)[scalar_dof_id] == false) { continue; }
3860 
3861  const IntegrationPoint &ip = ir.IntPoint(s);
3862  Tpr.SetIntPoint(&ip);
3863  double w = surf_fit_normal * surf_fit_coeff->Eval(Tpr, ip) *
3864  1.0 / surf_fit_dof_count[vdofs[s]];
3865 
3866  if (surf_fit_gf) { w *= 2.0 * sigma_e(s); }
3867  if (surf_fit_pos)
3868  {
3869  Vector pos(dim), pos_target(dim);
3870  for (int d = 0; d < dim; d++)
3871  {
3872  pos(d) = PMatI(s, d);
3873  pos_target(d) = (*surf_fit_pos)(vdofs[d*dof_s + s]);
3874  }
3875  Vector grad_s(dim);
3876  surf_fit_limiter->Eval_d1(pos, pos_target, 1.0, grad_s);
3877  for (int d = 0; d < dim; d++) { surf_fit_grad_e(s, d) = grad_s(d); }
3878  }
3879 
3880  for (int d = 0; d < dim; d++)
3881  {
3882  mat(s, d) += w * surf_fit_grad_e(s, d);
3883  }
3884  }
3885 }
3886 
3887 void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
3888  IsoparametricTransformation &Tpr,
3889  DenseMatrix &mat)
3890 {
3891  const int el_id = Tpr.ElementNo;
3892 
3893  // Scalar for surf_fit_gf, vector for surf_fit_pos, but that's ok.
3894  const FiniteElementSpace *fes_fit =
3895  (surf_fit_gf) ? surf_fit_gf->FESpace() : surf_fit_pos->FESpace();
3896  const FiniteElement &el_s = *fes_fit->GetFE(el_id);
3897  const int dof_s = el_s.GetDof(), dim = el_x.GetDim();
3898 
3899  // Check if the element has any DOFs marked for surface fitting.
3900  Array<int> dofs, vdofs;
3901  fes_fit->GetElementVDofs(el_id, vdofs);
3902  int count = 0;
3903  for (int s = 0; s < dof_s; s++)
3904  {
3905  // Because surf_fit_pos.fes might be ordered byVDIM.
3906  const int scalar_dof_id = fes_fit->VDofToDof(vdofs[s]);
3907  count += ((*surf_fit_marker)[scalar_dof_id]) ? 1 : 0;
3908  }
3909  if (count == 0) { return; }
3910 
3911  Vector sigma_e(dof_s);
3912  DenseMatrix surf_fit_grad_e(dof_s, dim);
3913  DenseMatrix surf_fit_hess_e(dof_s, dim*dim);
3914  if (surf_fit_gf || surf_fit_gf_bg)
3915  {
3916  surf_fit_gf->GetSubVector(vdofs, sigma_e);
3917 
3918  // Project the gradient of sigma in the same space.
3919  // The FE coefficients of the gradient go in surf_fit_grad_e.
3920  Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
3921  DenseMatrix grad_phys; // This will be (dof x dim, dof).
3922  if (surf_fit_gf_bg)
3923  {
3924  surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
3925  surf_fit_grad->GetSubVector(dofs, grad_ptr);
3926  }
3927  else
3928  {
3929  el_s.ProjectGrad(el_s, Tpr, grad_phys);
3930  grad_phys.Mult(sigma_e, grad_ptr);
3931  }
3932 
3933  // Project the Hessian of sigma in the same space.
3934  // The FE coefficients of the Hessian go in surf_fit_hess_e.
3935  Vector hess_ptr(surf_fit_hess_e.GetData(), dof_s*dim*dim);
3936  if (surf_fit_gf_bg)
3937  {
3938  surf_fit_hess->FESpace()->GetElementVDofs(el_id, dofs);
3939  surf_fit_hess->GetSubVector(dofs, hess_ptr);
3940  }
3941  else
3942  {
3943  surf_fit_hess_e.SetSize(dof_s*dim, dim);
3944  Mult(grad_phys, surf_fit_grad_e, surf_fit_hess_e);
3945  surf_fit_hess_e.SetSize(dof_s, dim * dim);
3946  }
3947  }
3948  else { Tpr.GetPointMat().Transpose(PMatI); }
3949 
3950  const IntegrationRule &ir = el_s.GetNodes();
3951 
3952  DenseMatrix surf_fit_hess_s(dim, dim);
3953  for (int s = 0; s < dof_s; s++)
3954  {
3955  // Because surf_fit_pos.fes might be ordered byVDIM.
3956  const int scalar_dof_id = fes_fit->VDofToDof(vdofs[s]);
3957  if ((*surf_fit_marker)[scalar_dof_id] == false) { continue; }
3958 
3959  const IntegrationPoint &ip = ir.IntPoint(s);
3960  Tpr.SetIntPoint(&ip);
3961  double w = surf_fit_normal * surf_fit_coeff->Eval(Tpr, ip);
3962 
3963  if (surf_fit_gf || surf_fit_gf_bg)
3964  {
3965  Vector gg_ptr(surf_fit_hess_s.GetData(), dim * dim);
3966  surf_fit_hess_e.GetRow(s, gg_ptr);
3967  w *= 2.0;
3968  }
3969  if (surf_fit_pos)
3970  {
3971  Vector pos(dim), pos_target(dim);
3972  for (int d = 0; d < dim; d++)
3973  {
3974  pos(d) = PMatI(s, d);
3975  pos_target(d) = (*surf_fit_pos)(vdofs[d*dof_s + s]);
3976  }
3977  // Eval_d2 returns the full Hessian, but we still use the general
3978  // computation that's in the dim x dim loop below.
3979  sigma_e(s) = 1.0;
3980  for (int d = 0; d < dim; d++) { surf_fit_grad_e(s, d) = 0.0; }
3981  surf_fit_limiter->Eval_d2(pos, pos_target, 1.0, surf_fit_hess_s);
3982  }
3983 
3984  // Loops over the local matrix.
3985  for (int idim = 0; idim < dim; idim++)
3986  {
3987  for (int jdim = 0; jdim <= idim; jdim++)
3988  {
3989  double entry = w * ( surf_fit_grad_e(s, idim) *
3990  surf_fit_grad_e(s, jdim) +
3991  sigma_e(s) * surf_fit_hess_s(idim, jdim));
3992  entry *= 1.0 / surf_fit_dof_count[vdofs[s]];
3993  int idx = s + idim*dof_s;
3994  int jdx = s + jdim*dof_s;
3995  mat(idx, jdx) += entry;
3996  if (idx != jdx) { mat(jdx, idx) += entry; }
3997  }
3998  }
3999  }
4000 }
4001 
4002 double TMOP_Integrator::GetFDDerivative(const FiniteElement &el,
4003  ElementTransformation &T,
4004  Vector &elfun, const int dofidx,
4005  const int dir, const double e_fx,
4006  bool update_stored)
4007 {
4008  int dof = el.GetDof();
4009  int idx = dir*dof+dofidx;
4010  elfun[idx] += dx;
4011  double e_fxph = GetElementEnergy(el, T, elfun);
4012  elfun[idx] -= dx;
4013  double dfdx = (e_fxph-e_fx)/dx;
4014 
4015  if (update_stored)
4016  {
4017  (*(ElemPertEnergy[T.ElementNo]))(idx) = e_fxph;
4018  (*(ElemDer[T.ElementNo]))(idx) = dfdx;
4019  }
4020 
4021  return dfdx;
4022 }
4023 
4024 void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
4025  ElementTransformation &T,
4026  const Vector &elfun,
4027  Vector &elvect)
4028 {
4029  const int dof = el.GetDof(), dim = el.GetDim(), elnum = T.ElementNo;
4030  if (elnum>=ElemDer.Size())
4031  {
4032  ElemDer.Append(new Vector);
4033  ElemPertEnergy.Append(new Vector);
4034  ElemDer[elnum]->SetSize(dof*dim);
4035  ElemPertEnergy[elnum]->SetSize(dof*dim);
4036  }
4037 
4038  elvect.SetSize(dof*dim);
4039  Vector elfunmod(elfun);
4040 
4041  // In GetElementEnergy(), skip terms that have exact derivative calculations.
4042  fd_call_flag = true;
4043 
4044  // Energy for unperturbed configuration.
4045  const double e_fx = GetElementEnergy(el, T, elfun);
4046 
4047  for (int j = 0; j < dim; j++)
4048  {
4049  for (int i = 0; i < dof; i++)
4050  {
4051  if (discr_tc)
4052  {
4053  discr_tc->UpdateTargetSpecificationAtNode(
4054  el, T, i, j, discr_tc->GetTspecPert1H());
4055  }
4056  elvect(j*dof+i) = GetFDDerivative(el, T, elfunmod, i, j, e_fx, true);
4057  if (discr_tc) { discr_tc->RestoreTargetSpecificationAtNode(T, i); }
4058  }
4059  }
4060  fd_call_flag = false;
4061 
4062  // Contributions from adaptive limiting, surface fitting (exact derivatives).
4063  if (adapt_lim_gf || surf_fit_gf || surf_fit_pos)
4064  {
4065  const IntegrationRule &ir = ActionIntegrationRule(el);
4066  const int nqp = ir.GetNPoints();
4067  DenseTensor Jtr(dim, dim, nqp);
4068  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
4069 
4070  IsoparametricTransformation Tpr;
4071  Tpr.SetFE(&el);
4072  Tpr.ElementNo = T.ElementNo;
4073  Tpr.Attribute = T.Attribute;
4074  Tpr.mesh = T.mesh;
4075  PMatI.UseExternalData(elfun.GetData(), dof, dim);
4076  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
4077 
4078  Vector weights(nqp);
4079  for (int q = 0; q < nqp; q++)
4080  {
4081  weights(q) = (integ_over_target) ?
4082  ir.IntPoint(q).weight * Jtr(q).Det() :
4083  ir.IntPoint(q).weight;
4084  }
4085 
4086  PMatO.UseExternalData(elvect.GetData(), dof, dim);
4087  if (adapt_lim_gf) { AssembleElemVecAdaptLim(el, Tpr, ir, weights, PMatO); }
4088  if (surf_fit_gf || surf_fit_pos) { AssembleElemVecSurfFit(el, Tpr, PMatO); }
4089  }
4090 }
4091 
4092 void TMOP_Integrator::AssembleElementGradFD(const FiniteElement &el,
4093  ElementTransformation &T,
4094  const Vector &elfun,
4095  DenseMatrix &elmat)
4096 {
4097  const int dof = el.GetDof(), dim = el.GetDim();
4098 
4099  elmat.SetSize(dof*dim);
4100  Vector elfunmod(elfun);
4101 
4102  const Vector &ElemDerLoc = *(ElemDer[T.ElementNo]);
4103  const Vector &ElemPertLoc = *(ElemPertEnergy[T.ElementNo]);
4104 
4105  // In GetElementEnergy(), skip terms that have exact derivative calculations.
4106  fd_call_flag = true;
4107  for (int i = 0; i < dof; i++)
4108  {
4109  for (int j = 0; j < i+1; j++)
4110  {
4111  for (int k1 = 0; k1 < dim; k1++)
4112  {
4113  for (int k2 = 0; k2 < dim; k2++)
4114  {
4115  elfunmod(k2*dof+j) += dx;
4116 
4117  if (discr_tc)
4118  {
4119  discr_tc->UpdateTargetSpecificationAtNode(
4120  el, T, j, k2, discr_tc->GetTspecPert1H());
4121  if (j != i)
4122  {
4123  discr_tc->UpdateTargetSpecificationAtNode(
4124  el, T, i, k1, discr_tc->GetTspecPert1H());
4125  }
4126  else // j==i
4127  {
4128  if (k1 != k2)
4129  {
4130  int idx = k1+k2-1;
4131  discr_tc->UpdateTargetSpecificationAtNode(
4132  el, T, i, idx, discr_tc->GetTspecPertMixH());
4133  }
4134  else // j==i && k1==k2
4135  {
4136  discr_tc->UpdateTargetSpecificationAtNode(
4137  el, T, i, k1, discr_tc->GetTspecPert2H());
4138  }
4139  }
4140  }
4141 
4142  double e_fx = ElemPertLoc(k2*dof+j);
4143  double e_fpxph = GetFDDerivative(el, T, elfunmod, i, k1, e_fx,
4144  false);
4145  elfunmod(k2*dof+j) -= dx;
4146  double e_fpx = ElemDerLoc(k1*dof+i);
4147 
4148  elmat(k1*dof+i, k2*dof+j) = (e_fpxph - e_fpx) / dx;
4149  elmat(k2*dof+j, k1*dof+i) = (e_fpxph - e_fpx) / dx;
4150 
4151  if (discr_tc)
4152  {
4153  discr_tc->RestoreTargetSpecificationAtNode(T, i);
4154  discr_tc->RestoreTargetSpecificationAtNode(T, j);
4155  }
4156  }
4157  }
4158  }
4159  }
4160  fd_call_flag = false;
4161 
4162  // Contributions from adaptive limiting.
4163  if (adapt_lim_gf || surf_fit_gf || surf_fit_pos)
4164  {
4165  const IntegrationRule &ir = GradientIntegrationRule(el);
4166  const int nqp = ir.GetNPoints();
4167  DenseTensor Jtr(dim, dim, nqp);
4168  targetC->ComputeElementTargets(T.ElementNo, el, ir, elfun, Jtr);
4169 
4170  IsoparametricTransformation Tpr;
4171  Tpr.SetFE(&el);
4172  Tpr.ElementNo = T.ElementNo;
4173  Tpr.Attribute = T.Attribute;
4174  Tpr.mesh = T.mesh;
4175  PMatI.UseExternalData(elfun.GetData(), dof, dim);
4176  Tpr.GetPointMat().Transpose(PMatI); // PointMat = PMatI^T
4177 
4178  Vector weights(nqp);
4179  for (int q = 0; q < nqp; q++)
4180  {
4181  weights(q) = (integ_over_target) ?
4182  ir.IntPoint(q).weight * Jtr(q).Det() :
4183  ir.IntPoint(q).weight;
4184  }
4185 
4186  if (adapt_lim_gf) { AssembleElemGradAdaptLim(el, Tpr, ir, weights, elmat); }
4187  if (surf_fit_gf || surf_fit_pos) { AssembleElemGradSurfFit(el, Tpr, elmat); }
4188  }
4189 }
4190 
4191 void TMOP_Integrator::UpdateSurfaceFittingWeight(double factor)
4192 {
4193  if (!surf_fit_coeff) { return; }
4194 
4195  if (surf_fit_coeff)
4196  {
4197  auto cf = dynamic_cast<ConstantCoefficient *>(surf_fit_coeff);
4198  MFEM_VERIFY(cf, "Dynamic weight works only with a ConstantCoefficient.");
4199  cf->constant *= factor;
4200  }
4201 }
4202 
4203 double TMOP_Integrator::GetSurfaceFittingWeight()
4204 {
4205  if (surf_fit_coeff)
4206  {
4207  auto cf = dynamic_cast<ConstantCoefficient *>(surf_fit_coeff);
4208  MFEM_VERIFY(cf, "Dynamic weight works only with a ConstantCoefficient.");
4209  return cf->constant;
4210  }
4211  return 0.0;
4212 }
4213 
4214 void TMOP_Integrator::EnableNormalization(const GridFunction &x)
4215 {
4216  ComputeNormalizationEnergies(x, metric_normal, lim_normal, surf_fit_normal);
4217  metric_normal = 1.0 / metric_normal;
4218  lim_normal = 1.0 / lim_normal;
4219  //if (surf_fit_gf) { surf_fit_normal = 1.0 / surf_fit_normal; }
4220  if (surf_fit_gf || surf_fit_pos) { surf_fit_normal = lim_normal; }
4221 }
4222 
4223 #ifdef MFEM_USE_MPI
4224 void TMOP_Integrator::ParEnableNormalization(const ParGridFunction &x)
4225 {
4226  double loc[3];
4227  ComputeNormalizationEnergies(x, loc[0], loc[1], loc[2]);
4228  double rdc[3];
4229  MPI_Allreduce(loc, rdc, 3, MPI_DOUBLE, MPI_SUM, x.ParFESpace()->GetComm());
4230  metric_normal = 1.0 / rdc[0];
4231  lim_normal = 1.0 / rdc[1];
4232  // if (surf_fit_gf) { surf_fit_normal = 1.0 / rdc[2]; }
4233  if (surf_fit_gf || surf_fit_pos) { surf_fit_normal = lim_normal; }
4234 }
4235 #endif
4236 
4237 void TMOP_Integrator::ComputeNormalizationEnergies(const GridFunction &x,
4238  double &metric_energy,
4239  double &lim_energy,
4240  double &surf_fit_gf_energy)
4241 {
4242  Array<int> vdofs;
4243  Vector x_vals;
4244  const FiniteElementSpace* const fes = x.FESpace();
4245 
4246  const int dim = fes->GetMesh()->Dimension();
4247  Jrt.SetSize(dim);
4248  Jpr.SetSize(dim);
4249  Jpt.SetSize(dim);
4250 
4251  metric_energy = 0.0;
4252  lim_energy = 0.0;
4253  surf_fit_gf_energy = 0.0;
4254  for (int i = 0; i < fes->GetNE(); i++)
4255  {
4256  const FiniteElement *fe = fes->GetFE(i);
4257  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
4258  const int nqp = ir.GetNPoints();
4259  DenseTensor Jtr(dim, dim, nqp);
4260  const int dof = fe->GetDof();
4261  DSh.SetSize(dof, dim);
4262 
4263  fes->GetElementVDofs(i, vdofs);
4264  x.GetSubVector(vdofs, x_vals);
4265  PMatI.UseExternalData(x_vals.GetData(), dof, dim);
4266 
4267  targetC->ComputeElementTargets(i, *fe, ir, x_vals, Jtr);
4268 
4269  for (int q = 0; q < nqp; q++)
4270  {
4271  const IntegrationPoint &ip = ir.IntPoint(q);
4272  metric->SetTargetJacobian(Jtr(q));
4273  CalcInverse(Jtr(q), Jrt);
4274  const double weight =
4275  (integ_over_target) ? ip.weight * Jtr(q).Det() : ip.weight;
4276 
4277  fe->CalcDShape(ip, DSh);
4278  MultAtB(PMatI, DSh, Jpr);
4279  Mult(Jpr, Jrt, Jpt);
4280 
4281  metric_energy += weight * metric->EvalW(Jpt);
4282  lim_energy += weight;
4283  }
4284 
4285  // Normalization of the surface fitting term.
4286  if (surf_fit_gf)
4287  {
4288  Array<int> dofs;
4289  Vector sigma_e;
4290  surf_fit_gf->FESpace()->GetElementDofs(i, dofs);
4291  surf_fit_gf->GetSubVector(dofs, sigma_e);
4292  for (int s = 0; s < dofs.Size(); s++)
4293  {
4294  if ((*surf_fit_marker)[dofs[s]] == true)
4295  {
4296  surf_fit_gf_energy += sigma_e(s) * sigma_e(s);
4297  }
4298  }
4299  }
4300  }
4301 
4302  // Cases when integration is not over the target element, or when the
4303  // targets don't contain volumetric information.
4304  if (integ_over_target == false || targetC->ContainsVolumeInfo() == false)
4305  {
4306  lim_energy = fes->GetNE();
4307  }
4308 }
4309 
4310 void TMOP_Integrator::ComputeMinJac(const Vector &x,
4311  const FiniteElementSpace &fes)
4312 {
4313  const FiniteElement *fe = fes.GetFE(0);
4314  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
4315  const int NE = fes.GetMesh()->GetNE(), dim = fe->GetDim(),
4316  dof = fe->GetDof(), nsp = ir.GetNPoints();
4317 
4318  Array<int> xdofs(dof * dim);
4319  DenseMatrix dshape(dof, dim), pos(dof, dim);
4320  Vector posV(pos.Data(), dof * dim);
4321  Jpr.SetSize(dim);
4322 
4323  dx = std::numeric_limits<float>::max();
4324 
4325  double detv_sum;
4326  double detv_avg_min = std::numeric_limits<float>::max();
4327  for (int i = 0; i < NE; i++)
4328  {
4329  fes.GetElementVDofs(i, xdofs);
4330  x.GetSubVector(xdofs, posV);
4331  detv_sum = 0.;
4332  for (int j = 0; j < nsp; j++)
4333  {
4334  fes.GetFE(i)->CalcDShape(ir.IntPoint(j), dshape);
4335  MultAtB(pos, dshape, Jpr);
4336  detv_sum += std::fabs(Jpr.Det());
4337  }
4338  double detv_avg = pow(detv_sum/nsp, 1./dim);
4339  detv_avg_min = std::min(detv_avg, detv_avg_min);
4340  }
4341  dx = detv_avg_min / dxscale;
4342 }
4343 
4344 void TMOP_Integrator::
4345 UpdateAfterMeshPositionChange(const Vector &x_new,
4346  const FiniteElementSpace &x_fes)
4347 {
4348  if (discr_tc) { PA.Jtr_needs_update = true; }
4349 
4350  Ordering::Type ordering = x_fes.GetOrdering();
4351 
4352  // Update the finite difference delta if FD are used.
4353  if (fdflag) { ComputeFDh(x_new, x_fes); }
4354 
4355  // Update the target constructor if it's a discrete one.
4356  if (discr_tc)
4357  {
4358  discr_tc->UpdateTargetSpecification(x_new, true, ordering);
4359  if (fdflag)
4360  {
4361  discr_tc->UpdateGradientTargetSpecification(x_new, dx, true, ordering);
4362  discr_tc->UpdateHessianTargetSpecification(x_new, dx, true, ordering);
4363  }
4364  }
4365 
4366  // Update adapt_lim_gf if adaptive limiting is enabled.
4367  if (adapt_lim_gf)
4368  {
4369  adapt_lim_eval->ComputeAtNewPosition(x_new, *adapt_lim_gf, ordering);
4370  }
4371 
4372  // Update surf_fit_gf if surface fitting is enabled.
4373  if (surf_fit_gf)
4374  {
4375  if (surf_fit_gf_bg)
4376  {
4377  // Interpolate information for only DOFs marked for fitting.
4378  const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension();
4379  const int cnt = surf_fit_marker_dof_index.Size();
4380  const int total_cnt = x_new.Size()/dim;
4381  Vector new_x_sorted(cnt*dim);
4382  if (ordering == 0)
4383  {
4384  for (int d = 0; d < dim; d++)
4385  {
4386  for (int i = 0; i < cnt; i++)
4387  {
4388  int dof_index = surf_fit_marker_dof_index[i];
4389  new_x_sorted(i + d*cnt) = x_new(dof_index + d*total_cnt);
4390  }
4391  }
4392  }
4393  else
4394  {
4395  for (int i = 0; i < cnt; i++)
4396  {
4397  int dof_index = surf_fit_marker_dof_index[i];
4398  for (int d = 0; d < dim; d++)
4399  {
4400  new_x_sorted(d + i*dim) = x_new(d + dof_index*dim);
4401  }
4402  }
4403  }
4404 
4405  Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int;
4406  surf_fit_eval->ComputeAtNewPosition(
4407  new_x_sorted, surf_fit_gf_int, ordering);
4408  for (int i = 0; i < cnt; i++)
4409  {
4410  int dof_index = surf_fit_marker_dof_index[i];
4411  (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i);
4412  }
4413 
4414  surf_fit_eval_bg_grad->ComputeAtNewPosition(
4415  new_x_sorted, surf_fit_grad_int, ordering);
4416  // Assumes surf_fit_grad and surf_fit_gf share the same space
4417  const int grad_dim = surf_fit_grad->VectorDim();
4418  const int grad_cnt = surf_fit_grad->Size()/grad_dim;
4419  if (surf_fit_grad->FESpace()->GetOrdering() == Ordering::byNODES)
4420  {
4421  for (int d = 0; d < grad_dim; d++)
4422  {
4423  for (int i = 0; i < cnt; i++)
4424  {
4425  int dof_index = surf_fit_marker_dof_index[i];
4426  (*surf_fit_grad)[dof_index + d*grad_cnt] =
4427  surf_fit_grad_int(i + d*cnt);
4428  }
4429  }
4430  }
4431  else
4432  {
4433  for (int i = 0; i < cnt; i++)
4434  {
4435  int dof_index = surf_fit_marker_dof_index[i];
4436  for (int d = 0; d < grad_dim; d++)
4437  {
4438  (*surf_fit_grad)[dof_index*dim + d] =
4439  surf_fit_grad_int(i*dim + d);
4440  }
4441  }
4442  }
4443 
4444  surf_fit_eval_bg_hess->ComputeAtNewPosition(
4445  new_x_sorted, surf_fit_hess_int, ordering);
4446  // Assumes surf_fit_hess and surf_fit_gf share the same space
4447  const int hess_dim = surf_fit_hess->VectorDim();
4448  const int hess_cnt = surf_fit_hess->Size()/hess_dim;
4449  if (surf_fit_hess->FESpace()->GetOrdering() == Ordering::byNODES)
4450  {
4451  for (int d = 0; d < hess_dim; d++)
4452  {
4453  for (int i = 0; i < cnt; i++)
4454  {
4455  int dof_index = surf_fit_marker_dof_index[i];
4456  (*surf_fit_hess)[dof_index + d*hess_cnt] =
4457  surf_fit_hess_int(i + d*cnt);
4458  }
4459  }
4460  }
4461  else
4462  {
4463  for (int i = 0; i < cnt; i++)
4464  {
4465  int dof_index = surf_fit_marker_dof_index[i];
4466  for (int d = 0; d < hess_dim; d++)
4467  {
4468  (*surf_fit_hess)[dof_index*dim + d] =
4469  surf_fit_hess_int(i*dim + d);
4470  }
4471  }
4472  }
4473  }
4474  else
4475  {
4476  surf_fit_eval->ComputeAtNewPosition(x_new, *surf_fit_gf, ordering);
4477  }
4478  }
4479 }
4480 
4481 void TMOP_Integrator::ComputeFDh(const Vector &x, const FiniteElementSpace &fes)
4482 {
4483  if (!fdflag) { return; }
4484  ComputeMinJac(x, fes);
4485 #ifdef MFEM_USE_MPI
4486  const ParFiniteElementSpace *pfes =
4487  dynamic_cast<const ParFiniteElementSpace *>(&fes);
4488  if (pfes)
4489  {
4490  double min_jac_all;
4491  MPI_Allreduce(&dx, &min_jac_all, 1, MPI_DOUBLE, MPI_MIN, pfes->GetComm());
4492  dx = min_jac_all;
4493  }
4494 #endif
4495 }
4496 
4497 void TMOP_Integrator::EnableFiniteDifferences(const GridFunction &x)
4498 {
4499  fdflag = true;
4500  const FiniteElementSpace *fes = x.FESpace();
4501  ComputeFDh(x,*fes);
4502  if (discr_tc)
4503  {
4504 #ifdef MFEM_USE_GSLIB
4505  const AdaptivityEvaluator *ae = discr_tc->GetAdaptivityEvaluator();
4506  if (dynamic_cast<const InterpolatorFP *>(ae))
4507  {
4508  MFEM_ABORT("Using GSLIB-based interpolation with finite differences"
4509  "requires careful consideration. Contact TMOP team.");
4510  }
4511 #endif
4512  discr_tc->UpdateTargetSpecification(x, false, fes->GetOrdering());
4513  discr_tc->UpdateGradientTargetSpecification(x, dx, false, fes->GetOrdering());
4514  discr_tc->UpdateHessianTargetSpecification(x, dx, false, fes->GetOrdering());
4515  }
4516 }
4517 
4518 #ifdef MFEM_USE_MPI
4519 void TMOP_Integrator::EnableFiniteDifferences(const ParGridFunction &x)
4520 {
4521  fdflag = true;
4522  const ParFiniteElementSpace *pfes = x.ParFESpace();
4523  ComputeFDh(x,*pfes);
4524  if (discr_tc)
4525  {
4526 #ifdef MFEM_USE_GSLIB
4527  const AdaptivityEvaluator *ae = discr_tc->GetAdaptivityEvaluator();
4528  if (dynamic_cast<const InterpolatorFP *>(ae))
4529  {
4530  MFEM_ABORT("Using GSLIB-based interpolation with finite differences"
4531  "requires careful consideration. Contact TMOP team.");
4532  }
4533 #endif
4534  discr_tc->UpdateTargetSpecification(x, false, pfes->GetOrdering());
4535  discr_tc->UpdateGradientTargetSpecification(x, dx, false, pfes->GetOrdering());
4536  discr_tc->UpdateHessianTargetSpecification(x, dx, false, pfes->GetOrdering());
4537  }
4538 }
4539 #endif
4540 
4541 double TMOP_Integrator::ComputeMinDetT(const Vector &x,
4542  const FiniteElementSpace &fes)
4543 {
4544  double min_detT = std::numeric_limits<double>::infinity();
4545  const int NE = fes.GetMesh()->GetNE();
4546  const int dim = fes.GetMesh()->Dimension();
4547  Array<int> xdofs;
4548  Jpr.SetSize(dim);
4549  Jpt.SetSize(dim);
4550  Jrt.SetSize(dim);
4551 
4552  for (int i = 0; i < NE; i++)
4553  {
4554  const FiniteElement *fe = fes.GetFE(i);
4555  const IntegrationRule &ir = EnergyIntegrationRule(*fe);
4556  const int dof = fe->GetDof(), nsp = ir.GetNPoints();
4557 
4558  DSh.SetSize(dof, dim);
4559  Vector posV(dof * dim);
4560  PMatI.UseExternalData(posV.GetData(), dof, dim);
4561 
4562  fes.GetElementVDofs(i, xdofs);
4563  x.GetSubVector(xdofs, posV);
4564 
4565  DenseTensor Jtr(dim, dim, ir.GetNPoints());
4566  targetC->ComputeElementTargets(i, *fe, ir, posV, Jtr);
4567 
4568  for (int q = 0; q < nsp; q++)
4569  {
4570  const IntegrationPoint &ip = ir.IntPoint(q);
4571  const DenseMatrix &Jtr_q = Jtr(q);
4572  CalcInverse(Jtr_q, Jrt);
4573  fe->CalcDShape(ip, DSh);
4574  MultAtB(PMatI, DSh, Jpr);
4575  Mult(Jpr,