MFEM  v4.5.2 Finite element discretization library
mesh-optimizer.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 void discrete_aspr_3d(const Vector &x, Vector &v)
88 {
89  int dim = x.Size();
90  v.SetSize(dim);
91  double l1, l2, l3;
92  l1 = 1.;
93  l2 = 1. + 5*x(1);
94  l3 = 1. + 10*x(2);
95  v[0] = l1/pow(l2*l3,0.5);
96  v[1] = l2/pow(l1*l3,0.5);
97  v[2] = l3/pow(l2*l1,0.5);
98 }
99
101 {
102 private:
103  int metric;
104
105 public:
106  HessianCoefficient(int dim, int metric_id)
107  : TMOPMatrixCoefficient(dim), metric(metric_id) { }
108
110  const IntegrationPoint &ip)
111  {
112  Vector pos(3);
113  T.Transform(ip, pos);
114  if (metric != 14 && metric != 36 && metric != 85)
115  {
116  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
117  const double r = sqrt(xc*xc + yc*yc);
118  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
119  const double eps = 0.5;
120
121  const double tan1 = std::tanh(sf*(r-r1)),
122  tan2 = std::tanh(sf*(r-r2));
123
124  K(0, 0) = eps + 1.0 * (tan1 - tan2);
125  K(0, 1) = 0.0;
126  K(1, 0) = 0.0;
127  K(1, 1) = 1.0;
128  }
129  else if (metric == 14 || metric == 36) // Size + Alignment
130  {
131  const double xc = pos(0), yc = pos(1);
132  double theta = M_PI * yc * (1.0 - yc) * cos(2 * M_PI * xc);
133  double alpha_bar = 0.1;
134
135  K(0, 0) = cos(theta);
136  K(1, 0) = sin(theta);
137  K(0, 1) = -sin(theta);
138  K(1, 1) = cos(theta);
139
140  K *= alpha_bar;
141  }
142  else if (metric == 85) // Shape + Alignment
143  {
144  Vector x = pos;
145  double xc = x(0)-0.5, yc = x(1)-0.5;
146  double th = 22.5*M_PI/180.;
147  double xn = cos(th)*xc + sin(th)*yc;
148  double yn = -sin(th)*xc + cos(th)*yc;
149  xc = xn; yc=yn;
150
151  double tfac = 20;
152  double s1 = 3;
153  double s2 = 2;
154  double wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1)
155  - std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) - 1);
156  if (wgt > 1) { wgt = 1; }
157  if (wgt < 0) { wgt = 0; }
158
159  xc = pos(0), yc = pos(1);
160  double theta = M_PI * (yc) * (1.0 - yc) * cos(2 * M_PI * xc);
161
162  K(0, 0) = cos(theta);
163  K(1, 0) = sin(theta);
164  K(0, 1) = -sin(theta);
165  K(1, 1) = cos(theta);
166
167  double asp_ratio_tar = 0.1 + 1*(1-wgt)*(1-wgt);
168
169  K(0, 0) *= 1/pow(asp_ratio_tar,0.5);
170  K(1, 0) *= 1/pow(asp_ratio_tar,0.5);
171  K(0, 1) *= pow(asp_ratio_tar,0.5);
172  K(1, 1) *= pow(asp_ratio_tar,0.5);
173  }
174  }
175
177  const IntegrationPoint &ip, int comp)
178  {
179  Vector pos(3);
180  T.Transform(ip, pos);
181  K = 0.;
182  if (metric != 14 && metric != 85)
183  {
184  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
185  const double r = sqrt(xc*xc + yc*yc);
186  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
187
188  const double tan1 = std::tanh(sf*(r-r1)),
189  tan2 = std::tanh(sf*(r-r2));
190  double tan1d = 0., tan2d = 0.;
191  if (r > 0.001)
192  {
193  tan1d = (1.-tan1*tan1)*(sf)/r,
194  tan2d = (1.-tan2*tan2)*(sf)/r;
195  }
196
197  K(0, 1) = 0.0;
198  K(1, 0) = 0.0;
199  K(1, 1) = 1.0;
200  if (comp == 0) { K(0, 0) = tan1d*xc - tan2d*xc; }
201  else if (comp == 1) { K(0, 0) = tan1d*yc - tan2d*yc; }
202  }
203  }
204 };
205
207 {
208 private:
209  int dim;
210  // 0 - size target in an annular region,
211  // 1 - size+aspect-ratio in an annular region,
212  // 2 - size+aspect-ratio target for a rotate sine wave.
213  int hr_target_type;
214
215 public:
216  HRHessianCoefficient(int dim_, int hr_target_type_ = 0)
217  : TMOPMatrixCoefficient(dim_), dim(dim_),
218  hr_target_type(hr_target_type_) { }
219
221  const IntegrationPoint &ip)
222  {
223  Vector pos(3);
224  T.Transform(ip, pos);
225  if (hr_target_type == 0) // size only circle
226  {
227  double small = 0.001, big = 0.01;
228  if (dim == 3) { small = 0.005, big = 0.1; }
229  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
230  double zc;
231  if (dim == 3) { zc = pos(2) - 0.5; }
232  double r = sqrt(xc*xc + yc*yc);
233  if (dim == 3) { r = sqrt(xc*xc + yc*yc + zc*zc); }
234  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
235
236  const double tan1 = std::tanh(sf*(r-r1)),
237  tan2 = std::tanh(sf*(r-r2));
238
239  double ind = (tan1 - tan2);
240  if (ind > 1.0) {ind = 1.;}
241  if (ind < 0.0) {ind = 0.;}
242  double val = ind * small + (1.0 - ind) * big;
243  K = 0.0;
244  K(0, 0) = 1.0;
245  K(0, 1) = 0.0;
246  K(1, 0) = 0.0;
247  K(1, 1) = 1.0;
248  K(0, 0) *= pow(val,0.5);
249  K(1, 1) *= pow(val,0.5);
250  if (dim == 3) { K(2, 2) = pow(val,0.5); }
251  }
252  else if (hr_target_type == 1) // circle with size and AR
253  {
254  const double small = 0.001, big = 0.01;
255  const double xc = pos(0)-0.5, yc = pos(1)-0.5;
256  const double rv = xc*xc + yc*yc;
257  double r = 0;
258  if (rv>0.) {r = sqrt(rv);}
259
260  double r1 = 0.2; double r2 = 0.3; double sf=30.0;
261  const double szfac = 1;
262  const double asfac = 4;
263  const double eps2 = szfac/asfac;
264  const double eps1 = szfac;
265
266  double tan1 = std::tanh(sf*(r-r1)+1),
267  tan2 = std::tanh(sf*(r-r2)-1);
268  double wgt = 0.5*(tan1-tan2);
269
270  tan1 = std::tanh(sf*(r-r1)),
271  tan2 = std::tanh(sf*(r-r2));
272
273  double ind = (tan1 - tan2);
274  if (ind > 1.0) {ind = 1.;}
275  if (ind < 0.0) {ind = 0.;}
276  double szval = ind * small + (1.0 - ind) * big;
277
278  double th = std::atan2(yc,xc)*180./M_PI;
279  if (wgt > 1) { wgt = 1; }
280  if (wgt < 0) { wgt = 0; }
281
282  double maxval = eps2 + eps1*(1-wgt)*(1-wgt);
283  double minval = eps1;
284  double avgval = 0.5*(maxval+minval);
285  double ampval = 0.5*(maxval-minval);
286  double val1 = avgval + ampval*sin(2.*th*M_PI/180.+90*M_PI/180.);
287  double val2 = avgval + ampval*sin(2.*th*M_PI/180.-90*M_PI/180.);
288
289  K(0,1) = 0.0;
290  K(1,0) = 0.0;
291  K(0,0) = val1;
292  K(1,1) = val2;
293
294  K(0,0) *= pow(szval,0.5);
295  K(1,1) *= pow(szval,0.5);
296  }
297  else if (hr_target_type == 2) // sharp rotated sine wave
298  {
299  double xc = pos(0)-0.5, yc = pos(1)-0.5;
300  double th = 15.5*M_PI/180.;
301  double xn = cos(th)*xc + sin(th)*yc;
302  double yn = -sin(th)*xc + cos(th)*yc;
303  double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
304  double stretch = 1/cos(th2);
305  xc = xn/stretch;
306  yc = yn;
307  double tfac = 20;
308  double s1 = 3;
309  double s2 = 2;
310  double yl1 = -0.025;
311  double yl2 = 0.025;
312  double wgt = std::tanh((tfac*(yc-yl1) + s2*std::sin(s1*M_PI*xc)) + 1) -
313  std::tanh((tfac*(yc-yl2) + s2*std::sin(s1*M_PI*xc)) - 1);
314  if (wgt > 1) { wgt = 1; }
315  if (wgt < 0) { wgt = 0; }
316
317  const double eps2 = 25;
318  const double eps1 = 1;
319  K(1,1) = eps1/eps2 + eps1*(1-wgt)*(1-wgt);
320  K(0,0) = eps1;
321  K(0,1) = 0.0;
322  K(1,0) = 0.0;
323  }
324  else { MFEM_ABORT("Unsupported option / wrong input."); }
325  }
326
328  const IntegrationPoint &ip, int comp)
329  {
330  K = 0.;
331  }
332 };
333
334 // Additional IntegrationRules that can be used with the --quad-type option.
337
338 // Defined with respect to the icf mesh.
339 double weight_fun(const Vector &x)
340 {
341  const double r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
342  const double den = 0.002;
343  double l2 = 0.2 + 0.5*std::tanh((r-0.16)/den) - 0.5*std::tanh((r-0.17)/den)
344  + 0.5*std::tanh((r-0.23)/den) - 0.5*std::tanh((r-0.24)/den);
345  return l2;
346 }
347
348 // Used for the adaptive limiting examples.
350 {
351  const double xc = x(0) - 0.1, yc = x(1) - 0.2;
352  const double r = sqrt(xc*xc + yc*yc);
353  double r1 = 0.45; double r2 = 0.55; double sf=30.0;
354  double val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
355
356  val = std::max(0.,val);
357  val = std::min(1.,val);
358  return val;
359 }
360
361 // Used for exact surface alignment
362 double surface_level_set(const Vector &x)
363 {
364  const int type = 1;
365
366  const int dim = x.Size();
367  if (type == 0)
368  {
369  const double sine = 0.25 * std::sin(4 * M_PI * x(0));
370  return (x(1) >= sine + 0.5) ? 1.0 : -1.0;
371  }
372  else
373  {
374  if (dim == 2)
375  {
376  const double xc = x(0) - 0.5, yc = x(1) - 0.5;
377  const double r = sqrt(xc*xc + yc*yc);
378  return std::tanh(2.0*(r-0.3));
379  }
380  else
381  {
382  const double xc = x(0) - 0.5, yc = x(1) - 0.5, zc = x(2) - 0.5;
383  const double r = sqrt(xc*xc + yc*yc + zc*zc);
384  return std::tanh(2.0*(r-0.3));
385  }
386  }
387 }
388
389 int material_id(int el_id, const GridFunction &g)
390 {
391  const FiniteElementSpace *fes = g.FESpace();
392  const FiniteElement *fe = fes->GetFE(el_id);
393  Vector g_vals;
394  const IntegrationRule &ir =
395  IntRules.Get(fe->GetGeomType(), fes->GetOrder(el_id) + 2);
396
397  double integral = 0.0;
398  g.GetValues(el_id, ir, g_vals);
400  for (int q = 0; q < ir.GetNPoints(); q++)
401  {
402  const IntegrationPoint &ip = ir.IntPoint(q);
403  Tr->SetIntPoint(&ip);
404  integral += ip.weight * g_vals(q) * Tr->Weight();
405  }
406  return (integral > 0.0) ? 1.0 : 0.0;
407 }
408
409 void DiffuseField(GridFunction &field, int smooth_steps)
410 {
411  // Setup the Laplacian operator
412  BilinearForm *Lap = new BilinearForm(field.FESpace());
414  Lap->Assemble();
415  Lap->Finalize();
416
417  // Setup the smoothing operator
418  DSmoother *S = new DSmoother(0,1.0,smooth_steps);
419  S->iterative_mode = true;
420  S->SetOperator(Lap->SpMat());
421
422  Vector tmp(field.Size());
423  tmp = 0.0;
424  S->Mult(tmp, field);
425
426  delete S;
427  delete Lap;
428 }
429
430 #ifdef MFEM_USE_MPI
431 void DiffuseField(ParGridFunction &field, int smooth_steps)
432 {
433  // Setup the Laplacian operator
434  ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
436  Lap->Assemble();
437  Lap->Finalize();
438  HypreParMatrix *A = Lap->ParallelAssemble();
439
440  HypreSmoother *S = new HypreSmoother(*A,0,smooth_steps);
441  S->iterative_mode = true;
442
443  Vector tmp(A->Width());
444  field.SetTrueVector();
445  Vector fieldtrue = field.GetTrueVector();
446  tmp = 0.0;
447  S->Mult(tmp, fieldtrue);
448
449  field.SetFromTrueDofs(fieldtrue);
450
451  delete S;
452  delete Lap;
453 }
454 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:232
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 ...
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
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)
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.
const SparseMatrix & SpMat() const
Returns a const reference to the sparse matrix: .
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:512
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)
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Container class for integration rules.
Definition: intrules.hpp:311
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:686
STL namespace.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:3637
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
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
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
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:524
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:169
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:971
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
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 ...
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
HessianCoefficient(int dim, int metric_id)
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
double weight_fun(const Vector &x)
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
int GetOrder(int i) const
Returns the polynomial degree of the i&#39;th finite element.
Definition: fespace.hpp:571
const Vector & GetTrueVector() const
Definition: gridfunc.hpp:130
double material_indicator_2d(const Vector &x)
aka closed Newton-Cotes
Definition: intrules.hpp:298
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
double discrete_size_3d(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.
int dim
Definition: ex24.cpp:53