MFEM  v4.6.0
Finite element discretization library
distance.cpp
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 // ------------------------------------------------
13 // Distance Miniapp: Finite element distance solver
14 // ------------------------------------------------
15 //
16 // This miniapp computes the "distance" to a given point source or to the zero
17 // level set of a given function. Here "distance" refers to the length of the
18 // shortest path through the mesh. The input can be a DeltaCoefficient (for a
19 // point source), or any Coefficient (for the case of a level set). The output
20 // is a GridFunction that can be scalar (representing the scalar distance), or a
21 // vector (its magnitude is the distance, and its direction is the starting
22 // direction of the shortest path). The miniapp supports 3 solvers:
23 //
24 // 1. Heat solver:
25 // K. Crane, C. Weischedel, M. Weischedel
26 // Geodesics in Heat: A New Approach to Computing Distance Based on Heat Flow
27 // ACM Transactions on Graphics, Vol. 32, No. 5, October, 2013
28 //
29 // 2. p-Laplacian solver:
30 // A. Belyaev, P. Fayolle
31 // On Variational and PDE-based Distance Function Approximations,
32 // Computer Graphics Forum, 34: 104-118, 2015, Section 7.
33 //
34 // 3. Rvachev normalization solver: same paper as p-Laplacian, Section 6.
35 // This solver is computationally cheap, but is accurate for distance
36 // approximations only near the zero level set.
37 //
38 //
39 // The solution of the p-Laplacian solver approaches the signed distance when
40 // p->\infinity. Therefore, increasing p will improve the computed distance and,
41 // of course, will increase the computational cost. The discretization of the
42 // p-Laplacian equation utilizes ideas from:
43 //
44 // L. V. Kantorovich, V. I. Krylov
45 // Approximate Methods of Higher Analysis, Interscience Publishers, Inc., 1958
46 //
47 // J. Melenk, I. Babuska
48 // The partition of unity finite element method: Basic theory and applications,
49 // Computer Methods in Applied Mechanics and Engineering, 1996, 139, 289-314
50 //
51 // Resolving highly oscillatory input fields requires refining the mesh or
52 // increasing the order of the approximation. The above requirement for mesh
53 // resolution is independent of the conditions imposed on the mesh by the
54 // discretization of the actual distance solver. On the other hand, it is often
55 // enough to compute the distance field to a mean zero level of a smoothed
56 // version of the input field. In such cases, one can use a low-pass filter that
57 // removes the high-frequency content of the input field, such as the one in the
58 // class PDEFilter, based on the Screened Poisson equation. The radius specifies
59 // the minimal feature size in the filter output and, in the current example, is
60 // linked to the average mesh size. Theoretical description of the filter and
61 // discussion of the parameters can be found in:
62 //
63 // B. S. Lazarov, O. Sigmund
64 // Filters in topology optimization based on Helmholtz-type differential equations
65 // International Journal for Numerical Methods in Engineering, 2011, 86, 765-781
66 //
67 // Compile with: make distance
68 //
69 // Sample runs:
70 //
71 // Problem 0: point source.
72 // mpirun -np 4 distance -m ./corners.mesh -p 0 -rs 3 -t 200.0
73 //
74 // Problem 1: zero level set: ball at the center of the domain - the exact
75 // distance is known, the code computes global and local errors.
76 // mpirun -np 4 distance -m ../../data/inline-segment.mesh -rs 3 -o 2 -t 1.0 -p 1
77 // mpirun -np 4 distance -m ../../data/inline-quad.mesh -rs 3 -o 2 -t 1.0 -p 1
78 // mpirun -np 4 distance -m ../../data/inline-hex.mesh -rs 1 -o 2 -p 1 -s 1
79 //
80 // Problem 2: zero level set: perturbed sine
81 // mpirun -np 4 distance -m ../../data/inline-quad.mesh -rs 3 -o 2 -t 1.0 -p 2
82 // mpirun -np 4 distance -m ../../data/amr-quad.mesh -rs 3 -o 2 -t 1.0 -p 2
83 //
84 // Problem 3: level set: Gyroid
85 // mpirun -np 4 distance -m ../../data/periodic-square.mesh -rs 5 -o 2 -t 1.0 -p 3
86 // mpirun -np 4 distance -m ../../data/periodic-cube.mesh -rs 3 -o 2 -t 1.0 -p 3 -s 2
87 //
88 // Problem 4: level set: Union of doughnut and swiss cheese shapes
89 // mpirun -np 4 distance -m ../../data/inline-hex.mesh -rs 3 -o 2 -t 1.0 -p 4
90 
91 #include <fstream>
92 #include <iostream>
93 #include "../common/mfem-common.hpp"
94 #include "sbm_aux.hpp"
95 
96 using namespace std;
97 using namespace mfem;
98 using namespace common;
99 
100 double sine_ls(const Vector &x)
101 {
102  const double sine = 0.25 * std::sin(4 * M_PI * x(0)) +
103  0.05 * std::sin(16 * M_PI * x(0));
104  return (x(1) >= sine + 0.5) ? -1.0 : 1.0;
105 }
106 
107 const double radius = 0.4;
108 
109 double sphere_ls(const Vector &x)
110 {
111  const int dim = x.Size();
112  const double xc = x(0) - 0.5;
113  const double yc = (dim > 1) ? x(1) - 0.5 : 0.0;
114  const double zc = (dim > 2) ? x(2) - 0.5 : 0.0;
115  const double r = sqrt(xc*xc + yc*yc + zc*zc);
116 
117  return (r >= radius) ? -1.0 : 1.0;
118 }
119 
120 double exact_dist_sphere(const Vector &x)
121 {
122  const int dim = x.Size();
123  const double xc = x(0) - 0.5;
124  const double yc = (dim > 1) ? x(1) - 0.5 : 0.0;
125  const double zc = (dim > 2) ? x(2) - 0.5 : 0.0;
126  const double r = sqrt(xc*xc + yc*yc + zc*zc);
127 
128  return fabs(r - radius);
129 }
130 
131 class ExactDistSphereLoc : public Coefficient
132 {
133 private:
134  ParGridFunction &dist;
135  const double dx;
136 
137 public:
138  ExactDistSphereLoc(ParGridFunction &d)
139  : dist(d), dx(dist.ParFESpace()->GetParMesh()->GetElementSize(0)) { }
140 
141  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
142  {
143  Vector pos(T.GetDimension());
144  T.Transform(ip, pos);
145  pos -= 0.5;
146  const double r = sqrt(pos * pos);
147 
148  // One zone length in every direction.
149  if (fabs(r - radius) < dx) { return fabs(r - radius); }
150  else { return dist.GetValue(T, ip); }
151  }
152 };
153 
154 
155 double Gyroid(const Vector &xx)
156 {
157  const double period = 2.0 * M_PI;
158  double x = xx[0]*period;
159  double y = xx[1]*period;
160  double z = (xx.Size()==3) ? xx[2]*period : 0.0;
161 
162  return std::sin(x)*std::cos(y) +
163  std::sin(y)*std::cos(z) +
164  std::sin(z)*std::cos(x);
165 }
166 
167 double Sph(const mfem::Vector &xx)
168 {
169  double R=0.4;
170  mfem::Vector lvec(3);
171  lvec=0.0;
172  for (int i=0; i<xx.Size(); i++)
173  {
174  lvec[i]=xx[i];
175  }
176 
177  return lvec[0]*lvec[0]+lvec[1]*lvec[1]+lvec[2]*lvec[2]-R*R;
178 }
179 
180 void DGyroid(const mfem::Vector &xx, mfem::Vector &vals)
181 {
182  vals.SetSize(xx.Size());
183  vals=0.0;
184 
185  double pp=4*M_PI;
186 
187  mfem::Vector lvec(3);
188  lvec=0.0;
189  for (int i=0; i<xx.Size(); i++)
190  {
191  lvec[i]=xx[i]*pp;
192  }
193 
194  vals[0]=cos(lvec[0])*cos(lvec[1])-sin(lvec[2])*sin(lvec[0]);
195  vals[1]=-sin(lvec[0])*sin(lvec[1])+cos(lvec[1])*cos(lvec[2]);
196  if (xx.Size()>2)
197  {
198  vals[2]=-sin(lvec[1])*sin(lvec[2])+cos(lvec[2])*cos(lvec[0]);
199  }
200 
201  vals*=pp;
202 }
203 
204 int main(int argc, char *argv[])
205 {
206  // Initialize MPI and HYPRE.
207  Mpi::Init(argc, argv);
208  int myid = Mpi::WorldRank();
209  Hypre::Init();
210 
211  // Parse command-line options.
212  const char *mesh_file = "../../data/inline-quad.mesh";
213  int solver_type = 0;
214  int problem = 1;
215  int rs_levels = 2;
216  int order = 2;
217  double t_param = 1.0;
218  const char *device_config = "cpu";
219  bool visualization = true;
220 
221  OptionsParser args(argc, argv);
222  args.AddOption(&mesh_file, "-m", "--mesh",
223  "Mesh file to use.");
224  args.AddOption(&solver_type, "-s", "--solver",
225  "Solver type:\n\t"
226  "0: Heat\n\t"
227  "1: P-Laplacian\n\t"
228  "2: Rvachev scaling");
229  args.AddOption(&problem, "-p", "--problem",
230  "Problem type:\n\t"
231  "0: Point source\n\t"
232  "1: Circle / sphere level set in 2D / 3D\n\t"
233  "2: 2D sine-looking level set\n\t"
234  "3: Gyroid level set in 2D or 3D\n\t"
235  "4: Combo of a doughnut and swiss cheese shapes in 3D.");
236  args.AddOption(&rs_levels, "-rs", "--refine-serial",
237  "Number of times to refine the mesh uniformly in serial.");
238  args.AddOption(&order, "-o", "--order",
239  "Finite element order (polynomial degree) or -1 for"
240  " isoparametric space.");
241  args.AddOption(&t_param, "-t", "--t-param",
242  "Diffusion time step (scaled internally scaled by dx*dx).");
243  args.AddOption(&device_config, "-d", "--device",
244  "Device configuration string, see Device::Configure().");
245  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
246  "--no-visualization",
247  "Enable or disable GLVis visualization.");
248  args.Parse();
249  if (!args.Good())
250  {
251  if (myid == 0)
252  {
253  args.PrintUsage(cout);
254  }
255  return 1;
256  }
257  if (myid == 0) { args.PrintOptions(cout); }
258 
259  // Enable hardware devices such as GPUs, and programming models such as CUDA,
260  // OCCA, RAJA and OpenMP based on command line options.
261  Device device(device_config);
262  if (myid == 0) { device.Print(); }
263 
264  // Refine the mesh.
265  Mesh mesh(mesh_file, 1, 1);
266  const int dim = mesh.Dimension();
267  for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
268 
269  // MPI distribution.
270  ParMesh pmesh(MPI_COMM_WORLD, mesh);
271  mesh.Clear();
272 
273  Coefficient *ls_coeff = nullptr;
274  int smooth_steps;
275  if (problem == 0)
276  {
277  ls_coeff = new DeltaCoefficient(0.5, -0.5, 1000.0);
278  smooth_steps = 0;
279  }
280  else if (problem == 1)
281  {
282  ls_coeff = new FunctionCoefficient(sphere_ls);
283  smooth_steps = 0;
284  }
285  else if (problem == 2)
286  {
287  ls_coeff = new FunctionCoefficient(sine_ls);
288  smooth_steps = 0;
289  }
290  else if (problem == 3)
291  {
292  ls_coeff = new FunctionCoefficient(Gyroid);
293  smooth_steps = 0;
294  }
295  else if (problem == 4)
296  {
297  ls_coeff = new FunctionCoefficient(doughnut_cheese);
298  smooth_steps = 0;
299  }
300  else { MFEM_ABORT("Unrecognized -problem option."); }
301 
302  const double dx = AvgElementSize(pmesh);
303  DistanceSolver *dist_solver = NULL;
304  if (solver_type == 0)
305  {
306  auto ds = new HeatDistanceSolver(t_param * dx * dx);
307  if (problem == 0)
308  {
309  ds->transform = false;
310  }
311  ds->smooth_steps = smooth_steps;
312  ds->vis_glvis = false;
313  dist_solver = ds;
314  }
315  else if (solver_type == 1)
316  {
317  const int p = 10;
318  const int newton_iter = 50;
319  dist_solver = new PLapDistanceSolver(p, newton_iter);
320  }
321  else if (solver_type == 2)
322  {
323  dist_solver = new NormalizationDistanceSolver;
324  }
325  else { MFEM_ABORT("Wrong solver option."); }
326  dist_solver->print_level.FirstAndLast().Summary();
327 
328  H1_FECollection fec(order, dim);
329  ParFiniteElementSpace pfes_s(&pmesh, &fec), pfes_v(&pmesh, &fec, dim);
330  ParGridFunction distance_s(&pfes_s), distance_v(&pfes_v);
331 
332  // Smooth-out Gibbs oscillations from the input level set. The smoothing
333  // parameter here is specified to be mesh dependent with length scale dx.
334  ParGridFunction filt_gf(&pfes_s);
335  if (problem != 0)
336  {
337  double filter_weight = dx;
338  // The normalization-based solver needs a more diffused input.
339  if (solver_type == 2) { filter_weight *= 4.0; }
340  PDEFilter filter(pmesh, filter_weight);
341  filter.Filter(*ls_coeff, filt_gf);
342  }
343  else { filt_gf.ProjectCoefficient(*ls_coeff); }
344  delete ls_coeff;
345  GridFunctionCoefficient ls_filt_coeff(&filt_gf);
346 
347  dist_solver->ComputeScalarDistance(ls_filt_coeff, distance_s);
348  dist_solver->ComputeVectorDistance(ls_filt_coeff, distance_v);
349 
350  // Send the solution by socket to a GLVis server.
351  if (visualization)
352  {
353  int size = 500;
354  char vishost[] = "localhost";
355  int visport = 19916;
356 
357  socketstream sol_sock_w;
358  common::VisualizeField(sol_sock_w, vishost, visport, filt_gf,
359  "Input Level Set", 0, 0, size, size);
360 
361  MPI_Barrier(pmesh.GetComm());
362 
363  socketstream sol_sock_ds;
364  common::VisualizeField(sol_sock_ds, vishost, visport, distance_s,
365  "Distance", size, 0, size, size,
366  "rRjmm********A");
367 
368  MPI_Barrier(pmesh.GetComm());
369 
370  socketstream sol_sock_dv;
371  common::VisualizeField(sol_sock_dv, vishost, visport, distance_v,
372  "Directions", 2*size, 0, size, size,
373  "rRjmm********vveA");
374  }
375 
376  // ParaView output.
377  ParaViewDataCollection dacol("ParaViewDistance", &pmesh);
378  dacol.SetLevelsOfDetail(order);
379  dacol.RegisterField("filtered_level_set", &filt_gf);
380  dacol.RegisterField("distance", &distance_s);
381  dacol.SetTime(1.0);
382  dacol.SetCycle(1);
383  dacol.Save();
384 
385  ConstantCoefficient zero(0.0);
386  const double s_norm = distance_s.ComputeL2Error(zero),
387  v_norm = distance_v.ComputeL2Error(zero);
388  if (myid == 0)
389  {
390  cout << fixed << setprecision(10) << "Norms: "
391  << s_norm << " " << v_norm << endl;
392  }
393 
394  if (problem == 1)
395  {
396  FunctionCoefficient exact_dist_coeff(exact_dist_sphere);
397  const double error_l1 = distance_s.ComputeL1Error(exact_dist_coeff),
398  error_li = distance_s.ComputeMaxError(exact_dist_coeff);
399  if (myid == 0)
400  {
401  cout << "Global L1 error: " << error_l1 << endl
402  << "Global Linf error: " << error_li << endl;
403  }
404 
405  ExactDistSphereLoc exact_dist_coeff_loc(distance_s);
406  const double error_l1_loc = distance_s.ComputeL1Error(exact_dist_coeff_loc),
407  error_li_loc = distance_s.ComputeMaxError(exact_dist_coeff_loc);
408  if (myid == 0)
409  {
410  cout << "Local L1 error: " << error_l1_loc << endl
411  << "Local Linf error: " << error_li_loc << endl;
412  }
413  }
414 
415  delete dist_solver;
416 
417  return 0;
418 }
int visport
int main(int argc, char *argv[])
Definition: distance.cpp:204
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
Helper class for ParaView visualization data.
double exact_dist_sphere(const Vector &x)
Definition: distance.cpp:120
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void Filter(ParGridFunction &func, ParGridFunction &ffield)
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
virtual void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)=0
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
double sine_ls(const Vector &x)
Definition: distance.cpp:100
double Sph(const mfem::Vector &xx)
Definition: distance.cpp:167
double sphere_ls(const Vector &x)
Definition: distance.cpp:109
void SetTime(double t)
Set physical time (for time-dependent simulations)
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:285
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
IterativeSolver::PrintLevel print_level
Definition: dist_solver.hpp:34
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
int problem
Definition: ex15.cpp:62
void DGyroid(const mfem::Vector &xx, mfem::Vector &vals)
Definition: distance.cpp:180
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: pgridfunc.cpp:273
Class for integration point with weight.
Definition: intrules.hpp:31
const double radius
Definition: distance.cpp:107
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
int dim
Definition: ex24.cpp:53
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:678
A general function coefficient.
Vector data type.
Definition: vector.hpp:58
double doughnut_cheese(const Vector &coord)
Definition: sbm_aux.hpp:32
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:259
virtual void Save() override
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
double Gyroid(const Vector &xx)
Definition: distance.cpp:155
Class for parallel meshes.
Definition: pmesh.hpp:32
double AvgElementSize(ParMesh &pmesh)
Definition: dist_solver.cpp:47