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