MFEM  v4.3.0 Finite element discretization library
mesh-optimizer.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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  // double th2 = (th > 45.*M_PI/180) ? M_PI/2 - th : th;
94  // double stretch = 1/cos(th2);
95  xc = xn; yc = yn;
96
97  double tfac = 20;
98  double s1 = 3;
99  double s2 = 2;
100  double wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1)
101  - std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) - 1);
102  if (wgt > 1) { wgt = 1; }
103  if (wgt < 0) { wgt = 0; }
104  return 0.1 + 1*(1-wgt)*(1-wgt);
105 }
106
107 void discrete_aspr_3d(const Vector &x, Vector &v)
108 {
109  int dim = x.Size();
110  v.SetSize(dim);
111  double l1, l2, l3;
112  l1 = 1.;
113  l2 = 1. + 5*x(1);
114  l3 = 1. + 10*x(2);
115  v[0] = l1/pow(l2*l3,0.5);
116  v[1] = l2/pow(l1*l3,0.5);
117  v[2] = l3/pow(l2*l1,0.5);
118 }
119
121 {
122 private:
123  int metric;
124
125 public:
126  HessianCoefficient(int dim, int metric_id)
127  : TMOPMatrixCoefficient(dim), metric(metric_id) { }
128
130  const IntegrationPoint &ip)
131  {
132  Vector pos(3);
133  T.Transform(ip, pos);
134  if (metric != 14 && metric != 36 && metric != 85)
135  {
136  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
137  const double r = sqrt(xc*xc + yc*yc);
138  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
139  const double eps = 0.5;
140
141  const double tan1 = std::tanh(sf*(r-r1)),
142  tan2 = std::tanh(sf*(r-r2));
143
144  K(0, 0) = eps + 1.0 * (tan1 - tan2);
145  K(0, 1) = 0.0;
146  K(1, 0) = 0.0;
147  K(1, 1) = 1.0;
148  }
149  else if (metric == 14 || metric == 36) // Size + Alignment
150  {
151  const double xc = pos(0), yc = pos(1);
152  double theta = M_PI * yc * (1.0 - yc) * cos(2 * M_PI * xc);
153  double alpha_bar = 0.1;
154
155  K(0, 0) = cos(theta);
156  K(1, 0) = sin(theta);
157  K(0, 1) = -sin(theta);
158  K(1, 1) = cos(theta);
159
160  K *= alpha_bar;
161  }
162  else if (metric == 85) // Shape + Alignment
163  {
164  Vector x = pos;
165  double xc = x(0)-0.5, yc = x(1)-0.5;
166  double th = 22.5*M_PI/180.;
167  double xn = cos(th)*xc + sin(th)*yc;
168  double yn = -sin(th)*xc + cos(th)*yc;
169  xc = xn; yc=yn;
170
171  double tfac = 20;
172  double s1 = 3;
173  double s2 = 2;
174  double wgt = std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) + 1)
175  - std::tanh((tfac*(yc) + s2*std::sin(s1*M_PI*xc)) - 1);
176  if (wgt > 1) { wgt = 1; }
177  if (wgt < 0) { wgt = 0; }
178
179  xc = pos(0), yc = pos(1);
180  double theta = M_PI * (yc) * (1.0 - yc) * cos(2 * M_PI * xc);
181
182  K(0, 0) = cos(theta);
183  K(1, 0) = sin(theta);
184  K(0, 1) = -sin(theta);
185  K(1, 1) = cos(theta);
186
187  double asp_ratio_tar = 0.1 + 1*(1-wgt)*(1-wgt);
188
189  K(0, 0) *= 1/pow(asp_ratio_tar,0.5);
190  K(1, 0) *= 1/pow(asp_ratio_tar,0.5);
191  K(0, 1) *= pow(asp_ratio_tar,0.5);
192  K(1, 1) *= pow(asp_ratio_tar,0.5);
193  }
194  }
195
197  const IntegrationPoint &ip, int comp)
198  {
199  Vector pos(3);
200  T.Transform(ip, pos);
201  K = 0.;
202  if (metric != 14 && metric != 85)
203  {
204  const double xc = pos(0) - 0.5, yc = pos(1) - 0.5;
205  const double r = sqrt(xc*xc + yc*yc);
206  double r1 = 0.15; double r2 = 0.35; double sf=30.0;
207
208  const double tan1 = std::tanh(sf*(r-r1)),
209  tan2 = std::tanh(sf*(r-r2));
210  double tan1d = 0., tan2d = 0.;
211  if (r > 0.001)
212  {
213  tan1d = (1.-tan1*tan1)*(sf)/r,
214  tan2d = (1.-tan2*tan2)*(sf)/r;
215  }
216
217  K(0, 1) = 0.0;
218  K(1, 0) = 0.0;
219  K(1, 1) = 1.0;
220  if (comp == 0) { K(0, 0) = tan1d*xc - tan2d*xc; }
221  else if (comp == 1) { K(0, 0) = tan1d*yc - tan2d*yc; }
222  }
223  }
224 };
225
226 // Additional IntegrationRules that can be used with the --quad-type option.
229
230 // Defined with respect to the icf mesh.
231 double weight_fun(const Vector &x)
232 {
233  const double r = sqrt(x(0)*x(0) + x(1)*x(1) + 1e-12);
234  const double den = 0.002;
235  double l2 = 0.2 + 0.5*std::tanh((r-0.16)/den) - 0.5*std::tanh((r-0.17)/den)
236  + 0.5*std::tanh((r-0.23)/den) - 0.5*std::tanh((r-0.24)/den);
237  return l2;
238 }
239
240 // Used for the adaptive limiting examples.
242 {
243  const double xc = x(0) - 0.1, yc = x(1) - 0.2;
244  const double r = sqrt(xc*xc + yc*yc);
245  double r1 = 0.45; double r2 = 0.55; double sf=30.0;
246  double val = 0.5*(1+std::tanh(sf*(r-r1))) - 0.5*(1+std::tanh(sf*(r-r2)));
247
248  val = std::max(0.,val);
249  val = std::min(1.,val);
250  return val;
251 }
252
253 void DiffuseField(GridFunction &field, int smooth_steps)
254 {
255  //Setup the Laplacian operator
256  BilinearForm *Lap = new BilinearForm(field.FESpace());
258  Lap->Assemble();
259  Lap->Finalize();
260
261  //Setup the smoothing operator
262  DSmoother *S = new DSmoother(0,1.0,smooth_steps);
263  S->iterative_mode = true;
264  S->SetOperator(Lap->SpMat());
265
266  Vector tmp(field.Size());
267  tmp = 0.0;
268  S->Mult(tmp, field);
269
270  delete S;
271  delete Lap;
272 }
273
274 #ifdef MFEM_USE_MPI
275 void DiffuseField(ParGridFunction &field, int smooth_steps)
276 {
277  //Setup the Laplacian operator
278  ParBilinearForm *Lap = new ParBilinearForm(field.ParFESpace());
280  Lap->Assemble();
281  Lap->Finalize();
282  HypreParMatrix *A = Lap->ParallelAssemble();
283
284  HypreSmoother *S = new HypreSmoother(*A,0,smooth_steps);
285  S->iterative_mode = true;
286
287  Vector tmp(A->Width());
288  field.SetTrueVector();
289  Vector fieldtrue = field.GetTrueVector();
290  tmp = 0.0;
291  S->Mult(tmp, fieldtrue);
292
293  field.SetFromTrueDofs(fieldtrue);
294
295  delete S;
296  delete Lap;
297 }
298 #endif
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 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.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
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.
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:190
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:652
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Relax the linear system Ax=b.
Definition: hypre.cpp:3222
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition: gridfunc.hpp:137
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: pgridfunc.hpp:164
void Assemble(int skip_zeros=1)
Assemble the local matrix.
double discrete_ori_2d(const Vector &x)
Parallel smoothers in hypre.
Definition: hypre.hpp:840
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629
void DiffuseField(ParGridFunction &field, int smooth_steps)
Definition: dist_solver.cpp:17
const Vector & GetTrueVector() const
Definition: gridfunc.hpp:127
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 ...
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication with Jacobi smoother.
HessianCoefficient(int dim, int metric_id)
double weight_fun(const Vector &x)
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
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
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
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:108