MFEM  v4.4.0 Finite element discretization library
mesh-optimizer.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
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 // MFEM Mesh Optimizer Miniapp - Serial/Parallel Shared Code
13
14 #include "mfem.hpp"
15 #include <fstream>
16 #include <iostream>
17
18 using namespace mfem;
19 using namespace std;
20
21 double discrete_size_2d(const Vector &x)
22 {
23  int opt = 2;
24  const double small = 0.001, big = 0.01;
25  double val = 0.;
26
27  if (opt == 1) // sine wave.
28  {
29  const double X = x(0), Y = x(1);
30  val = std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) + 1) -
31  std::tanh((10*(Y-0.5) + std::sin(4.0*M_PI*X)) - 1);
32  }
33  else if (opt == 2) // semi-circle
34  {
35  const double xc = x(0) - 0.0, yc = x(1) - 0.5;
36  const double r = sqrt(xc*xc + yc*yc);
37  double r1 = 0.45; double r2 = 0.55; double sf=30.0;
38  val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
39  }
40
41  val = std::max(0.,val);
42  val = std::min(1.,val);
43
44  return val * small + (1.0 - val) * big;
45 }
46
47 double discrete_size_3d(const Vector &x)
48 {
49  const double small = 0.0001, big = 0.01;
50  double val = 0.;
51
52  // semi-circle
53  const double xc = x(0) - 0.0, yc = x(1) - 0.5, zc = x(2) - 0.5;
54  const double r = sqrt(xc*xc + yc*yc + zc*zc);
55  double r1 = 0.45; double r2 = 0.55; double sf=30.0;
56  val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
57
58  val = std::max(0.,val);
59  val = std::min(1.,val);
60
61  return val * small + (1.0 - val) * big;
62 }
63
64 double material_indicator_2d(const Vector &x)
65 {
66  double xc = x(0)-0.5, yc = x(1)-0.5;
67  double th = 22.5*M_PI/180.;
68  double xn = cos(th)*xc + sin(th)*yc;
69  double yn = -sin(th)*xc + cos(th)*yc;
70  double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
71  double stretch = 1/cos(th2);
72  xc = xn/stretch; yc = yn/stretch;
73  double tfac = 20;
74  double s1 = 3;
75  double s2 = 3;
76  double wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1);
77  if (wgt > 1) { wgt = 1; }
78  if (wgt < 0) { wgt = 0; }
79  return wgt;
80 }
81
82 double discrete_ori_2d(const Vector &x)
83 {
84  return M_PI * x(1) * (1.0 - x(1)) * cos(2 * M_PI * x(0));
85 }
86
87 double discrete_aspr_2d(const Vector &x)
88 {
89  double xc = x(0)-0.5, yc = x(1)-0.5;
90  double th = 22.5*M_PI/180.;
91  double xn = cos(th)*xc + sin(th)*yc;
92  double yn = -sin(th)*xc + cos(th)*yc;
93  xc = xn; yc = yn;
94
95  double tfac = 20;
96  double s1 = 3;
97  double s2 = 2;
98  double wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1)
99  - std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) - 1);
100  if (wgt > 1) { wgt = 1; }
101  if (wgt < 0) { wgt = 0; }
102  return 0.1 + 1*(1-wgt)*(1-wgt);
103 }
104
105 void discrete_aspr_3d(const Vector &x, Vector &v)
106 {
107  int dim = x.Size();
108  v.SetSize(dim);
109  double l1, l2, l3;
110  l1 = 1.;
111  l2 = 1. + 5*x(1);
112  l3 = 1. + 10*x(2);
113  v[0] = l1/pow(l2*l3,0.5);
114  v[1] = l2/pow(l1*l3,0.5);
115  v[2] = l3/pow(l2*l1,0.5);
116 }
117
119 {
120 private:
121  int metric;
122
123 public:
124  HessianCoefficient(int dim, int metric_id)
125  : TMOPMatrixCoefficient(dim), metric(metric_id) { }
126
128  const IntegrationPoint &ip)
129  {
130  Vector pos(3);
131  T.Transform(ip, pos);
132  if (metric != 14 && metric != 36 && metric != 85)
133  {
134  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
135  const double r = sqrt(xc*xc + yc*yc);
136  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
137  const double eps = 0.5;
138
139  const double tan1 = std::tanh(sf*(r-r1)),
140  tan2 = std::tanh(sf*(r-r2));
141
142  K(0, 0) = eps + 1.0 * (tan1 - tan2);
143  K(0, 1) = 0.0;
144  K(1, 0) = 0.0;
145  K(1, 1) = 1.0;
146  }
147  else if (metric == 14 || metric == 36) // Size + Alignment
148  {
149  const double xc = pos(0), yc = pos(1);
150  double theta = M_PI * yc * (1.0 - yc) * cos(2 * M_PI * xc);
151  double alpha_bar = 0.1;
152
153  K(0, 0) = cos(theta);
154  K(1, 0) = sin(theta);
155  K(0, 1) = -sin(theta);
156  K(1, 1) = cos(theta);
157
158  K *= alpha_bar;
159  }
160  else if (metric == 85) // Shape + Alignment
161  {
162  Vector x = pos;
163  double xc = x(0)-0.5, yc = x(1)-0.5;
164  double th = 22.5*M_PI/180.;
165  double xn = cos(th)*xc + sin(th)*yc;
166  double yn = -sin(th)*xc + cos(th)*yc;
167  xc = xn; yc=yn;
168
169  double tfac = 20;
170  double s1 = 3;
171  double s2 = 2;
172  double wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1)
173  - std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) - 1);
174  if (wgt > 1) { wgt = 1; }
175  if (wgt < 0) { wgt = 0; }
176
177  xc = pos(0), yc = pos(1);
178  double theta = M_PI * (yc) * (1.0 - yc) * cos(2 * M_PI * xc);
179
180  K(0, 0) = cos(theta);
181  K(1, 0) = sin(theta);
182  K(0, 1) = -sin(theta);
183  K(1, 1) = cos(theta);
184
185  double asp_ratio_tar = 0.1 + 1*(1-wgt)*(1-wgt);
186
187  K(0, 0) *= 1/pow(asp_ratio_tar,0.5);
188  K(1, 0) *= 1/pow(asp_ratio_tar,0.5);
189  K(0, 1) *= pow(asp_ratio_tar,0.5);
190  K(1, 1) *= pow(asp_ratio_tar,0.5);
191  }
192  }
193
195  const IntegrationPoint &ip, int comp)
196  {
197  Vector pos(3);
198  T.Transform(ip, pos);
199  K = 0.;
200  if (metric != 14 && metric != 85)
201  {
202  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
203  const double r = sqrt(xc*xc + yc*yc);
204  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
205
206  const double tan1 = std::tanh(sf*(r-r1)),
207  tan2 = std::tanh(sf*(r-r2));
208  double tan1d = 0., tan2d = 0.;
209  if (r > 0.001)
210  {
211  tan1d = (1.-tan1*tan1)*(sf)/r,
212  tan2d = (1.-tan2*tan2)*(sf)/r;
213  }
214
215  K(0, 1) = 0.0;
216  K(1, 0) = 0.0;
217  K(1, 1) = 1.0;
218  if (comp == 0) { K(0, 0) = tan1d*xc - tan2d*xc; }
219  else if (comp == 1) { K(0, 0) = tan1d*yc - tan2d*yc; }
220  }
221  }
222 };
223
225 {
226 private:
227  int dim;
228  // 0 - size target in an annular region,
229  // 1 - size+aspect-ratio in an annular region,
230  // 2 - size+aspect-ratio target for a rotate sine wave.
231  int hr_target_type;
232
233 public:
234  HRHessianCoefficient(int dim_, int hr_target_type_ = 0)
235  : TMOPMatrixCoefficient(dim_), dim(dim_),
236  hr_target_type(hr_target_type_) { }
237
239  const IntegrationPoint &ip)
240  {
241  Vector pos(3);
242  T.Transform(ip, pos);
243  if (hr_target_type == 0) // size only circle
244  {
245  double small = 0.001, big = 0.01;
246  if (dim == 3) { small = 0.005, big = 0.1; }
247  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
248  double zc;
249  if (dim == 3) { zc = pos(2) - 0.5; }
250  double r = sqrt(xc*xc + yc*yc);
251  if (dim == 3) { r = sqrt(xc*xc + yc*yc + zc*zc); }
252  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
253
254  const double tan1 = std::tanh(sf*(r-r1)),
255  tan2 = std::tanh(sf*(r-r2));
256
257  double ind = (tan1 - tan2);
258  if (ind > 1.0) {ind = 1.;}
259  if (ind < 0.0) {ind = 0.;}
260  double val = ind * small + (1.0 - ind) * big;
261  K = 0.0;
262  K(0, 0) = 1.0;
263  K(0, 1) = 0.0;
264  K(1, 0) = 0.0;
265  K(1, 1) = 1.0;
266  K(0, 0) *= pow(val,0.5);
267  K(1, 1) *= pow(val,0.5);
268  if (dim == 3) { K(2, 2) = pow(val,0.5); }
269  }
270  else if (hr_target_type == 1) // circle with size and AR
271  {
272  const double small = 0.001, big = 0.01;
273  const double xc = pos(0)-0.5, yc = pos(1)-0.5;
274  const double rv = xc*xc + yc*yc;
275  double r = 0;
276  if (rv>0.) {r = sqrt(rv);}
277
278  double r1 = 0.2; double r2 = 0.3; double sf=30.0;
279  const double szfac = 1;
280  const double asfac = 4;
281  const double eps2 = szfac/asfac;
282  const double eps1 = szfac;
283
284  double tan1 = std::tanh(sf*(r-r1)+1),
285  tan2 = std::tanh(sf*(r-r2)-1);
286  double wgt = 0.5*(tan1-tan2);
287
288  tan1 = std::tanh(sf*(r-r1)),
289  tan2 = std::tanh(sf*(r-r2));
290
291  double ind = (tan1 - tan2);
292  if (ind > 1.0) {ind = 1.;}
293  if (ind < 0.0) {ind = 0.;}
294  double szval = ind * small + (1.0 - ind) * big;
295
296  double th = std::atan2(yc,xc)*180./M_PI;
297  if (wgt > 1) { wgt = 1; }
298  if (wgt < 0) { wgt = 0; }
299
300  double maxval = eps2 + eps1*(1-wgt)*(1-wgt);
301  double minval = eps1;
302  double avgval = 0.5*(maxval+minval);
303  double ampval = 0.5*(maxval-minval);
304  double val1 = avgval + ampval*sin(2.*th*M_PI/180.+90*M_PI/180.);
305  double val2 = avgval + ampval*sin(2.*th*M_PI/180.-90*M_PI/180.);
306
307  K(0,1) = 0.0;
308  K(1,0) = 0.0;
309  K(0,0) = val1;
310  K(1,1) = val2;
311
312  K(0,0) *= pow(szval,0.5);
313  K(1,1) *= pow(szval,0.5);
314  }
315  else if (hr_target_type == 2) // sharp rotated sine wave
316  {
317  double xc = pos(0)-0.5, yc = pos(1)-0.5;
318  double th = 15.5*M_PI/180.;
319  double xn = cos(th)*xc + sin(th)*yc;
320  double yn = -sin(th)*xc + cos(th)*yc;
321  double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
322  double stretch = 1/cos(th2);
323  xc = xn/stretch;
324  yc = yn;
325  double tfac = 20;
326  double s1 = 3;
327  double s2 = 2;
328  double yl1 = -0.025;
329  double yl2 = 0.025;
330  double wgt = std::tanh((tfac*(yc-yl1) + s2*std::sin(s1*M_PI*xc)) + 1) -
331  std::tanh((tfac*(yc-yl2) + s2*std::sin(s1*M_PI*xc)) - 1);
332  if (wgt > 1) { wgt = 1; }
333  if (wgt < 0) { wgt = 0; }
334
335  const double eps2 = 25;
336  const double eps1 = 1;
337  K(1,1) = eps1/eps2 + eps1*(1-wgt)*(1-wgt);
338  K(0,0) = eps1;
339  K(0,1) = 0.0;
340  K(1,0) = 0.0;
341  }
342  else { MFEM_ABORT("Unsupported option / wrong input."); }
343  }
344
346  const IntegrationPoint &ip, int comp)
347  {
348  K = 0.;
349  }
350 };
351
352 // Additional IntegrationRules that can be used with the --quad-type option.
355
356 // Defined with respect to the icf mesh.
357 double weight_fun(const Vector &x)
358 {
359  const double r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
360  const double den = 0.002;
361  double l2 = 0.2 + 0.5*std::tanh((r-0.16)/den) - 0.5*std::tanh((r-0.17)/den)
362  + 0.5*std::tanh((r-0.23)/den) - 0.5*std::tanh((r-0.24)/den);
363  return l2;
364 }
365
366 // Used for the adaptive limiting examples.
368 {
369  const double xc = x(0) - 0.1, yc = x(1) - 0.2;
370  const double r = sqrt(xc*xc + yc*yc);
371  double r1 = 0.45; double r2 = 0.55; double sf=30.0;
372  double val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
373
374  val = std::max(0.,val);
375  val = std::min(1.,val);
376  return val;
377 }
378
379 // Used for exact surface alignment
380 double surface_level_set(const Vector &x)
381 {
382  const int type = 1;
383
384  const int dim = x.Size();
385  if (type == 0)
386  {
387  const double sine = 0.25 * std::sin(4 * M_PI * x(0));
388  return (x(1) >= sine + 0.5) ? 1.0 : -1.0;
389  }
390  else
391  {
392  if (dim == 2)
393  {
394  const double xc = x(0) - 0.5, yc = x(1) - 0.5;
395  const double r = sqrt(xc*xc + yc*yc);
396  return std::tanh(2.0*(r-0.3));
397  }
398  else
399  {
400  const double xc = x(0) - 0.5, yc = x(1) - 0.5, zc = x(2) - 0.5;
401  const double r = sqrt(xc*xc + yc*yc + zc*zc);
402  return std::tanh(2.0*(r-0.3));
403  }
404  }
405 }
406
407 int material_id(int el_id, const GridFunction &g)
408 {
409  const FiniteElementSpace *fes = g.FESpace();
410  const FiniteElement *fe = fes->GetFE(el_id);
411  Vector g_vals;
412  const IntegrationRule &ir =
413  IntRules.Get(fe->GetGeomType(), fes->GetOrder(el_id) + 2);
414
415  double integral = 0.0;
416  g.GetValues(el_id, ir, g_vals);
418  for (int q = 0; q < ir.GetNPoints(); q++)
419  {
420  const IntegrationPoint &ip = ir.IntPoint(q);
421  Tr->SetIntPoint(&ip);
422  integral += ip.weight * g_vals(q) * Tr->Weight();
423  }
424  return (integral > 0.0) ? 1.0 : 0.0;
425 }
426
427 void DiffuseField(GridFunction &field, int smooth_steps)
428 {
429  // Setup the Laplacian operator
430  BilinearForm *Lap = new BilinearForm(field.FESpace());
432  Lap->Assemble();
433  Lap->Finalize();
434
435  // Setup the smoothing operator
436  DSmoother *S = new DSmoother(0,1.0,smooth_steps);
437  S->iterative_mode = true;
438  S->SetOperator(Lap->SpMat());
439
440  Vector tmp(field.Size());
441  tmp = 0.0;
442  S->Mult(tmp, field);
443
444  delete S;
445  delete Lap;
446 }
447
448 #ifdef MFEM_USE_MPI
449 void DiffuseField(ParGridFunction &field, int smooth_steps)
450 {
451  // Setup the Laplacian operator
452  ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
454  Lap->Assemble();
455  Lap->Finalize();
456  HypreParMatrix *A = Lap->ParallelAssemble();
457
458  HypreSmoother *S = new HypreSmoother(*A,0,smooth_steps);
459  S->iterative_mode = true;
460
461  Vector tmp(A->Width());
462  field.SetTrueVector();
463  Vector fieldtrue = field.GetTrueVector();
464  tmp = 0.0;
465  S->Mult(tmp, fieldtrue);
466
467  field.SetFromTrueDofs(fieldtrue);
468
469  delete S;
470  delete Lap;
471 }
472 #endif
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
void discrete_aspr_3d(const Vector &x, Vector &v)
virtual void EvalGrad(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip, int comp)
Evaluate the derivative of the matrix coefficient with respect to comp in the element described by T ...
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
double discrete_size_2d(const Vector &x)
FDualNumber< tbase > tanh(const FDualNumber< tbase > &f)
tanh([dual number])
Definition: fdual.hpp:632
Data type for scaled Jacobi-type smoother of sparse matrix.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void EvalGrad(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip, int comp)
Evaluate the derivative of the matrix coefficient with respect to comp in the element described by T ...
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual void SetOperator(const Operator &a)
Set/update the solver for the given operator.
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
HRHessianCoefficient(int dim_, int hr_target_type_=0)
double surface_level_set(const Vector &x)
aka closed Newton-Cotes
Definition: intrules.hpp:298
Container class for integration rules.
Definition: intrules.hpp:311
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:655
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:3426
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:144
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
void Assemble(int skip_zeros=1)
Assemble the local matrix.
double discrete_ori_2d(const Vector &x)
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
Parallel smoothers in hypre.
Definition: hypre.hpp:896
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:510
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:652
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
const Vector & GetTrueVector() const
Definition: gridfunc.hpp:134
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with Jacobi smoother.
HessianCoefficient(int dim, int metric_id)
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
double weight_fun(const Vector &x)
int GetOrder(int i) const
Returns the polynomial degree of the i&#39;th finite element.
Definition: fespace.hpp:547
double material_indicator_2d(const Vector &x)
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
double discrete_size_3d(const Vector &x)
double discrete_aspr_2d(const Vector &x)
Class for integration point with weight.
Definition: intrules.hpp:25
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2783
Class for parallel bilinear form.
Vector data type.
Definition: vector.hpp:60
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix.
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:337
int material_id(int el_id, const GridFunction &g)