MFEM  v4.5.2 Finite element discretization library
dist_solver.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 #ifndef MFEM_DIST_SOLVER_HPP
13 #define MFEM_DIST_SOLVER_HPP
14
15 #include "mfem.hpp"
16
17 namespace mfem
18 {
19
20 double AvgElementSize(ParMesh &pmesh);
21
23 {
24 protected:
25  void ScalarDistToVector(ParGridFunction &dist_s, ParGridFunction &dist_v);
26
27 public:
28  // 0 = nothing, 1 = main solver only, 2 = full (solver + preconditioner).
30
32  virtual ~DistanceSolver() { }
33
34  // Computes a scalar ParGridFunction which is the length of the shortest path
35  // to the zero level set of the given Coefficient. It is expected that the
36  // given [distance] has a valid (scalar) ParFiniteElementSpace, and that the
37  // result is computed in the same space. Some implementations may output a
38  // "signed" distance, i.e., the distance has different signs on both sides of
39  // the zero level set.
40  virtual void ComputeScalarDistance(Coefficient &zero_level_set,
41  ParGridFunction &distance) = 0;
42
43  // Computes a vector ParGridFunction where the magnitude is the length of the
44  // shortest path to the zero level set of the given Coefficient, and the
45  // direction is the starting direction of the shortest path. It is expected
46  // that the given [distance] has a valid (vector) ParFiniteElementSpace, and
47  // that the result is computed in the same space.
48  virtual void ComputeVectorDistance(Coefficient &zero_level_set,
49  ParGridFunction &distance);
50 };
51
52
53 // K. Crane et al: "Geodesics in Heat: A New Approach to Computing Distance
54 // Based on Heat Flow", DOI:10.1145/2516971.2516977.
56 {
57 public:
58  HeatDistanceSolver(double diff_coeff)
59  : DistanceSolver(), parameter_t(diff_coeff), smooth_steps(0),
60  diffuse_iter(1), transform(true), vis_glvis(false) { }
61
62  // The computed distance is not "signed". In addition to the standard usage
63  // (with zero level sets), this function can be applied to point sources when
64  // transform = false.
65  void ComputeScalarDistance(Coefficient &zero_level_set,
66  ParGridFunction &distance);
67
68  double parameter_t;
71 };
72
73 // A. Belyaev et al: "On Variational and PDE-based Distance Function
74 // Approximations", Section 6, DOI:10.1111/cgf.12611.
75 // This solver is computationally cheap, but is accurate for distance
76 // approximations only near the zero level set.
78 {
79 private:
80
81  class NormalizationCoeff : public Coefficient
82  {
83  private:
85
86  public:
87  NormalizationCoeff(ParGridFunction &u_gf) : u(u_gf) { }
88  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
89  };
90
91 public:
93
95 };
96
97
98 // A. Belyaev et al: "On Variational and PDE-based Distance Function
99 // Approximations", Section 7, DOI:10.1111/cgf.12611.
101 {
102 public:
103  PLapDistanceSolver(int maxp_ = 30, int newton_iter_ = 10,
104  double rtol = 1e-7, double atol = 1e-12)
105  : maxp(maxp_), newton_iter(newton_iter_),
106  newton_rel_tol(rtol), newton_abs_tol(atol) { }
107
108  void SetMaxPower(int new_pp) { maxp = new_pp; }
109
110  // The computed distance is "signed".
112
113 private:
114  int maxp; // maximum value of the power p
115  const int newton_iter;
116  const double newton_rel_tol, newton_abs_tol;
117 };
118
120 {
121 private:
122  const GridFunction &u;
123
124 public:
126  : VectorCoefficient(dim), u(u_gf) { }
127
129
131  {
132  T.SetIntPoint(&ip);
133
135  const double norm = V.Norml2() + 1e-12;
136  V /= -norm;
137  }
138 };
139
140
141 // Product of the modulus of the first coefficient and the second coefficient
143 {
144 private:
145  Coefficient &basef, &corrf;
146
147 public:
149  : basef(basec_), corrf(corrc_) { }
150
151  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
152  {
153  T.SetIntPoint(&ip);
154  double u = basef.Eval(T,ip);
155  double c = corrf.Eval(T,ip);
156  if (u<0.0) { u*=-1.0; }
157  return u*c;
158  }
159 };
160
161
162 // Formulation for the ScreenedPoisson equation. The positive part of the input
163 // coefficient supply unit volumetric loading, the negative part - negative unit
165 // which will deliver similar smoothing effect as the Screened Poisson
166 // equation. It determines the length scale of the smoothing.
168 {
169 protected:
170  double diffcoef;
172
173 public:
174  ScreenedPoisson(Coefficient &nfunc, double rh):func(&nfunc)
175  {
176  double rd=rh/(2*std::sqrt(3.0));
177  diffcoef= rd*rd;
178  }
179
181
182  void SetInput(Coefficient &nfunc) { func = &nfunc; }
183
184  virtual double GetElementEnergy(const FiniteElement &el,
186  const Vector &elfun) override;
187
188  virtual void AssembleElementVector(const FiniteElement &el,
190  const Vector &elfun,
191  Vector &elvect) override;
192
193  virtual void AssembleElementGrad(const FiniteElement &el,
195  const Vector &elfun,
196  DenseMatrix &elmat) override;
197 };
198
199
201 {
202
203 protected:
206  bool ownership;
207  double pp, ee;
208
209 public:
210  // The VectorCoefficent should contain a vector with entries:
211  // [0] - derivative with respect to x
212  // [1] - derivative with respect to y
213  // [2] - derivative with respect to z
215  bool ownership_=true)
217
218  void SetPower(double pp_) { pp = pp_; }
219  void SetReg(double ee_) { ee = ee_; }
220
221  virtual ~PUMPLaplacian()
222  {
223  if (ownership)
224  {
225  delete func;
227  }
228  }
229
230  virtual double GetElementEnergy(const FiniteElement &el,
232  const Vector &elfun) override;
233
234  virtual void AssembleElementVector(const FiniteElement &el,
236  const Vector &elfun,
237  Vector &elvect) override;
238
239  virtual void AssembleElementGrad(const FiniteElement &el,
241  const Vector &elfun,
242  DenseMatrix &elmat) override;
243 };
244
245 // Low-pass filter based on the Screened Poisson equation.
246 // B. S. Lazarov, O. Sigmund: "Filters in topology optimization based on
247 // Helmholtz-type differential equations", DOI:10.1002/nme.3072.
249 {
250 public:
251  PDEFilter(ParMesh &mesh, double rh, int order = 2,
252  int maxiter = 100, double rtol = 1e-12,
253  double atol = 1e-15, int print_lv = 0)
254  : rr(rh),
255  fecp(order, mesh.Dimension()),
256  fesp(&mesh, &fecp, 1),
257  gf(&fesp)
258  {
259  sv = fesp.NewTrueDofVector();
260
261  nf = new ParNonlinearForm(&fesp);
262  prec = new HypreBoomerAMG();
263  prec->SetPrintLevel(print_lv);
264
265  gmres = new GMRESSolver(mesh.GetComm());
266
267  gmres->SetAbsTol(atol);
268  gmres->SetRelTol(rtol);
269  gmres->SetMaxIter(maxiter);
270  gmres->SetPrintLevel(print_lv);
271  gmres->SetPreconditioner(*prec);
272
273  sint=nullptr;
274  }
275
277  {
278  delete gmres;
279  delete prec;
280  delete nf;
281  delete sv;
282  }
283
285  {
286  GridFunctionCoefficient gfc(&func);
287  Filter(gfc, ffield);
288  }
289
290  void Filter(Coefficient &func, ParGridFunction &ffield);
291
292 private:
293  const double rr;
294  H1_FECollection fecp;
296  ParGridFunction gf;
297
298  ParNonlinearForm* nf;
299  HypreBoomerAMG* prec;
300  GMRESSolver *gmres;
301  HypreParVector *sv;
302
303  ScreenedPoisson* sint;
304 };
305
306 } // namespace mfem
307
308 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:232
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
HypreParVector * NewTrueDofVector()
Definition: pfespace.hpp:331
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void ComputeScalarDistance(Coefficient &func, ParGridFunction &fdist)
Base class for vector Coefficients that optionally depend on time and space.
Coefficient * func
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
virtual ~DistanceSolver()
Definition: dist_solver.hpp:32
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Parallel non-linear operator on the true dofs.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
HeatDistanceSolver(double diff_coeff)
Definition: dist_solver.hpp:58
void ScalarDistToVector(ParGridFunction &dist_s, ParGridFunction &dist_v)
Definition: dist_solver.cpp:72
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
PUMPLaplacian(Coefficient *nfunc, VectorCoefficient *nfgrad, bool ownership_=true)
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
void SetInput(Coefficient &nfunc)
void SetMaxIter(int max_it)
Definition: solvers.hpp:201
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun) override
Compute the local energy.
ScreenedPoisson(Coefficient &nfunc, double rh)
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
PProductCoefficient(Coefficient &basec_, Coefficient &corrc_)
IterativeSolver::PrintLevel print_level
Definition: dist_solver.hpp:29
void SetAbsTol(double atol)
Definition: solvers.hpp:200
void SetRelTol(double rtol)
Definition: solvers.hpp:199
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
GMRES method.
Definition: solvers.hpp:525
virtual void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)
Class for integration point with weight.
Definition: intrules.hpp:25
This class is used to express the local action of a general nonlinear finite element operator...
Definition: nonlininteg.hpp:27
int dim
Definition: ex24.cpp:53
PLapDistanceSolver(int maxp_=30, int newton_iter_=10, double rtol=1e-7, double atol=1e-12)
void SetReg(double ee_)
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun) override
Compute the local energy.
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:804
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
Vector data type.
Definition: vector.hpp:60
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
virtual ~PUMPLaplacian()
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
void SetPower(double pp_)
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetMaxPower(int new_pp)
Class for parallel meshes.
Definition: pmesh.hpp:32
void ComputeScalarDistance(Coefficient &u_coeff, ParGridFunction &dist)
Settings for the output behavior of the IterativeSolver.
Definition: solvers.hpp:78
void Filter(ParGridFunction &func, ParGridFunction &ffield)
PDEFilter(ParMesh &mesh, double rh, int order=2, int maxiter=100, double rtol=1e-12, double atol=1e-15, int print_lv=0)
double AvgElementSize(ParMesh &pmesh)
Definition: dist_solver.cpp:42
virtual void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)=0