MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
dist_solver.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 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).
29  int print_level = 0;
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 
74 // A. Belyaev et al: "On Variational and PDE-based Distance Function
75 // Approximations", Section 7, DOI:10.1111/cgf.12611.
77 {
78 public:
79  PLapDistanceSolver(int maxp_ = 30, int newton_iter_ = 10,
80  double rtol = 1e-7, double atol = 1e-12, int print_lv = 0)
81  : maxp(maxp_), newton_iter(newton_iter_),
82  newton_rel_tol(rtol), newton_abs_tol(atol)
83  {
84  print_level = print_lv;
85  }
86 
87  void SetMaxPower(int new_pp) { maxp = new_pp; }
88 
89  // The computed distance is "signed".
91 
92 private:
93  int maxp; // maximum value of the power p
94  const int newton_iter;
95  const double newton_rel_tol, newton_abs_tol;
96 };
97 
98 
100 {
101 private:
102  const GridFunction &u;
103 
104 public:
106  : VectorCoefficient(dim), u(u_gf) { }
107 
109 
111  {
112  T.SetIntPoint(&ip);
113 
114  u.GetGradient(T, V);
115  const double norm = V.Norml2() + 1e-12;
116  V /= -norm;
117  }
118 };
119 
120 
121 // Product of the modulus of the first coefficient and the second coefficient
123 {
124 private:
125  Coefficient &basef, &corrf;
126 
127 public:
129  : basef(basec_), corrf(corrc_) { }
130 
131  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
132  {
133  T.SetIntPoint(&ip);
134  double u = basef.Eval(T,ip);
135  double c = corrf.Eval(T,ip);
136  if (u<0.0) { u*=-1.0; }
137  return u*c;
138  }
139 };
140 
141 
142 // Formulation for the ScreenedPoisson equation. The positive part of the input
143 // coefficient supply unit volumetric loading, the negative part - negative unit
144 // volumetric loading. The parameter rh is the radius of a linear cone filter
145 // which will deliver similar smoothing effect as the Screened Poisson
146 // equation. It determines the length scale of the smoothing.
148 {
149 protected:
150  double diffcoef;
152 
153 public:
154  ScreenedPoisson(Coefficient &nfunc, double rh):func(&nfunc)
155  {
156  double rd=rh/(2*std::sqrt(3.0));
157  diffcoef= rd*rd;
158  }
159 
161 
162  void SetInput(Coefficient &nfunc) { func = &nfunc; }
163 
164  virtual double GetElementEnergy(const FiniteElement &el,
166  const Vector &elfun) override;
167 
168  virtual void AssembleElementVector(const FiniteElement &el,
170  const Vector &elfun,
171  Vector &elvect) override;
172 
173  virtual void AssembleElementGrad(const FiniteElement &el,
175  const Vector &elfun,
176  DenseMatrix &elmat) override;
177 };
178 
179 
181 {
182 
183 protected:
186  bool ownership;
187  double pp, ee;
188 
189 public:
190  // The VectorCoefficent should contain a vector with entries:
191  // [0] - derivative with respect to x
192  // [1] - derivative with respect to y
193  // [2] - derivative with respect to z
195  bool ownership_=true)
196  : func(nfunc), fgrad(nfgrad), ownership(ownership_), pp(2.0), ee(1e-7) { }
197 
198  void SetPower(double pp_) { pp = pp_; }
199  void SetReg(double ee_) { ee = ee_; }
200 
201  virtual ~PUMPLaplacian()
202  {
203  if (ownership)
204  {
205  delete func;
206  delete fgrad;
207  }
208  }
209 
210  virtual double GetElementEnergy(const FiniteElement &el,
212  const Vector &elfun) override;
213 
214  virtual void AssembleElementVector(const FiniteElement &el,
216  const Vector &elfun,
217  Vector &elvect) override;
218 
219  virtual void AssembleElementGrad(const FiniteElement &el,
221  const Vector &elfun,
222  DenseMatrix &elmat) override;
223 };
224 
225 // Low-pass filter based on the Screened Poisson equation.
226 // B. S. Lazarov, O. Sigmund: "Filters in topology optimization based on
227 // Helmholtz-type differential equations", DOI:10.1002/nme.3072.
229 {
230 public:
231  PDEFilter(ParMesh &mesh, double rh, int order = 2,
232  int maxiter = 100, double rtol = 1e-12,
233  double atol = 1e-15, int print_lv = 0)
234  : rr(rh),
235  fecp(order, mesh.Dimension()),
236  fesp(&mesh, &fecp, 1),
237  gf(&fesp)
238  {
239  sv = fesp.NewTrueDofVector();
240 
241  nf = new ParNonlinearForm(&fesp);
242  prec = new HypreBoomerAMG();
243  prec->SetPrintLevel(print_lv);
244 
245  gmres = new GMRESSolver(mesh.GetComm());
246 
247  gmres->SetAbsTol(atol);
248  gmres->SetRelTol(rtol);
249  gmres->SetMaxIter(maxiter);
250  gmres->SetPrintLevel(print_lv);
251  gmres->SetPreconditioner(*prec);
252 
253  sint=nullptr;
254  }
255 
257  {
258  delete gmres;
259  delete prec;
260  delete nf;
261  delete sv;
262  }
263 
265  {
266  GridFunctionCoefficient gfc(&func);
267  Filter(gfc, ffield);
268  }
269 
270  void Filter(Coefficient &func, ParGridFunction &ffield);
271 
272 private:
273  const double rr;
274  H1_FECollection fecp;
276  ParGridFunction gf;
277 
278  ParNonlinearForm* nf;
279  HypreBoomerAMG* prec;
280  GMRESSolver *gmres;
281  HypreParVector *sv;
282 
283  ScreenedPoisson* sint;
284 };
285 
286 } // namespace mfem
287 
288 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:235
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 ...
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:792
virtual ~DistanceSolver()
Definition: dist_solver.hpp:32
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
virtual void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
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
PLapDistanceSolver(int maxp_=30, int newton_iter_=10, double rtol=1e-7, double atol=1e-12, int print_lv=0)
Definition: dist_solver.hpp:79
void ScalarDistToVector(ParGridFunction &dist_s, ParGridFunction &dist_v)
Definition: dist_solver.cpp:72
VectorCoefficient * fgrad
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:1443
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:1523
void SetInput(Coefficient &nfunc)
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun) override
Compute the local energy.
ScreenedPoisson(Coefficient &nfunc, double rh)
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:148
PProductCoefficient(Coefficient &basec_, Coefficient &corrc_)
void SetAbsTol(double atol)
Definition: solvers.hpp:199
void SetRelTol(double rtol)
Definition: solvers.hpp:198
MPI_Comm GetComm() const
Definition: pmesh.hpp:288
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
GMRES method.
Definition: solvers.hpp:497
NormalizedGradCoefficient(const GridFunction &u_gf, int dim)
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
Assemble the local gradient matrix.
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
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
int dim
Definition: ex24.cpp:53
void SetReg(double ee_)
virtual double GetElementEnergy(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun) override
Compute the local energy.
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:216
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)
Definition: dist_solver.hpp:87
Class for parallel meshes.
Definition: pmesh.hpp:32
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)
void GetGradient(ElementTransformation &tr, Vector &grad) const
Definition: gridfunc.cpp:1641
double AvgElementSize(ParMesh &pmesh)
Definition: dist_solver.cpp:42
virtual void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)=0