MFEM  v4.4.0 Finite element discretization library
distance.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 // ------------------------------------------------
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 2 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
33 //
34 // The solution of the p-Laplacian solver approaches the signed distance when
35 // p->\infinity. Therefore, increasing p will improve the computed distance and,
36 // of course, will increase the computational cost. The discretization of the
37 // p-Laplacian equation utilizes ideas from:
38 //
39 // L. V. Kantorovich, V. I. Krylov
40 // Approximate Methods of Higher Analysis, Interscience Publishers, Inc., 1958
41 //
42 // J. Melenk, I. Babuska
43 // The partition of unity finite element method: Basic theory and applications,
44 // Computer Methods in Applied Mechanics and Engineering, 1996, 139, 289-314
45 //
46 // Resolving highly oscillatory input fields requires refining the mesh or
47 // increasing the order of the approximation. The above requirement for mesh
48 // resolution is independent of the conditions imposed on the mesh by the
49 // discretization of the actual distance solver. On the other hand, it is often
50 // enough to compute the distance field to a mean zero level of a smoothed
51 // version of the input field. In such cases, one can use a low-pass filter that
52 // removes the high-frequency content of the input field, such as the one in the
53 // class PDEFilter, based on the Screened Poisson equation. The radius specifies
54 // the minimal feature size in the filter output and, in the current example, is
55 // linked to the average mesh size. Theoretical description of the filter and
56 // discussion of the parameters can be found in:
57 //
58 // B. S. Lazarov, O. Sigmund
59 // Filters in topology optimization based on Helmholtz-type differential equations
60 // International Journal for Numerical Methods in Engineering, 2011, 86, 765-781
61 //
62 // Compile with: make distance
63 //
64 // Sample runs:
65 //
66 // Problem 0: point source.
67 // mpirun -np 4 distance -m ./corners.mesh -p 0 -rs 3 -t 200.0
68 //
69 // Problem 1: zero level set: circle / sphere at the center of the mesh
70 // mpirun -np 4 distance -m ../../data/inline-quad.mesh -rs 3 -o 2 -t 1.0 -p 1
71 // mpirun -np 4 distance -m ../../data/periodic-cube.mesh -rs 2 -o 2 -p 1 -s 1
72 //
73 // Problem 2: zero level set: perturbed sine
74 // mpirun -np 4 distance -m ../../data/inline-quad.mesh -rs 3 -o 2 -t 1.0 -p 2
75 // mpirun -np 4 distance -m ../../data/amr-quad.mesh -rs 3 -o 2 -t 1.0 -p 2
76 //
77 // Problem 3: level set: Gyroid
78 // mpirun -np 4 distance -m ../../data/periodic-square.mesh -rs 5 -o 2 -t 1.0 -p 3
79 // mpirun -np 4 distance -m ../../data/periodic-cube.mesh -rs 3 -o 2 -t 1.0 -p 3
80 //
81 // Problem 4: level set: Union of doughnut and swiss cheese shapes
82 // mpirun -np 4 distance -m ../../data/inline-hex.mesh -rs 3 -o 2 -t 1.0 -p 4
83
84 #include <fstream>
85 #include <iostream>
86 #include "../common/mfem-common.hpp"
87 #include "dist_solver.hpp"
88 #include "sbm_aux.hpp"
89
90 using namespace std;
91 using namespace mfem;
92
93 double sine_ls(const Vector &x)
94 {
95  const double sine = 0.25 * std::sin(4 * M_PI * x(0)) +
96  0.05 * std::sin(16 * M_PI * x(0));
97  return (x(1) >= sine + 0.5) ? -1.0 : 1.0;
98 }
99
100 double sphere_ls(const Vector &x)
101 {
102  const int dim = x.Size();
103  if (dim == 2)
104  {
105  const double xc = x(0) - 0.5, yc = x(1) - 0.5;
106  const double r = sqrt(xc*xc + yc*yc);
107  return (r >= 0.4) ? -1.0 : 1.0;
108  }
109  else if (dim == 3)
110  {
111  const double xc = x(0) - 0.0, yc = x(1) - 0.0, zc = x(2) - 0.0;
112  const double r = sqrt(xc*xc + yc*yc + zc*zc);
113  return (r >= 0.8) ? -1.0 : 1.0;
114  }
115  else
116  {
117  return (x(0) >= 0.5) ? -1.0 : 1.0;
118  }
119 }
120
121 double Gyroid(const Vector &xx)
122 {
123  const double period = 2.0 * M_PI;
124  double x=xx[0]*period;
125  double y=xx[1]*period;
126  double z=0.0;
127  if (xx.Size()==3)
128  {
129  z=xx[2]*period;
130  }
131  return std::sin(x)*std::cos(y) +
132  std::sin(y)*std::cos(z) +
133  std::sin(z)*std::cos(x);
134 }
135
136 double Sph(const mfem::Vector &xx)
137 {
138  double R=0.4;
139  mfem::Vector lvec(3);
140  lvec=0.0;
141  for (int i=0; i<xx.Size(); i++)
142  {
143  lvec[i]=xx[i];
144  }
145
146  return lvec[0]*lvec[0]+lvec[1]*lvec[1]+lvec[2]*lvec[2]-R*R;
147 }
148
149 void DGyroid(const mfem::Vector &xx, mfem::Vector &vals)
150 {
151  vals.SetSize(xx.Size());
152  vals=0.0;
153
154  double pp=4*M_PI;
155
156  mfem::Vector lvec(3);
157  lvec=0.0;
158  for (int i=0; i<xx.Size(); i++)
159  {
160  lvec[i]=xx[i]*pp;
161  }
162
163  vals[0]=cos(lvec[0])*cos(lvec[1])-sin(lvec[2])*sin(lvec[0]);
164  vals[1]=-sin(lvec[0])*sin(lvec[1])+cos(lvec[1])*cos(lvec[2]);
165  if (xx.Size()>2)
166  {
167  vals[2]=-sin(lvec[1])*sin(lvec[2])+cos(lvec[2])*cos(lvec[0]);
168  }
169
170  vals*=pp;
171 }
172
173 int main(int argc, char *argv[])
174 {
175  // Initialize MPI and HYPRE.
176  Mpi::Init(argc, argv);
177  int myid = Mpi::WorldRank();
178  Hypre::Init();
179
180  // Parse command-line options.
181  const char *mesh_file = "../../data/inline-quad.mesh";
182  int solver_type = 0;
183  int problem = 1;
184  int rs_levels = 2;
185  int order = 2;
186  double t_param = 1.0;
187  const char *device_config = "cpu";
188  bool visualization = true;
189
190  OptionsParser args(argc, argv);
192  "Mesh file to use.");
194  "Solver type:\n\t"
195  "0: Heat\n\t"
196  "1: P-Laplacian");
198  "Problem type:\n\t"
199  "0: Point source\n\t"
200  "1: Circle / sphere level set in 2D / 3D\n\t"
201  "2: 2D sine-looking level set\n\t"
202  "3: Gyroid level set in 2D or 3D\n\t"
203  "4: Combo of a doughnut and swiss cheese shapes in 3D.");
205  "Number of times to refine the mesh uniformly in serial.");
207  "Finite element order (polynomial degree) or -1 for"
208  " isoparametric space.");
210  "Diffusion time step (scaled internally scaled by dx*dx).");
212  "Device configuration string, see Device::Configure().");
214  "--no-visualization",
215  "Enable or disable GLVis visualization.");
216  args.Parse();
217  if (!args.Good())
218  {
219  if (myid == 0)
220  {
221  args.PrintUsage(cout);
222  }
223  return 1;
224  }
225  if (myid == 0) { args.PrintOptions(cout); }
226
227  // Enable hardware devices such as GPUs, and programming models such as CUDA,
228  // OCCA, RAJA and OpenMP based on command line options.
229  Device device(device_config);
230  if (myid == 0) { device.Print(); }
231
232  // Refine the mesh.
233  Mesh mesh(mesh_file, 1, 1);
234  const int dim = mesh.Dimension();
235  for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
236
237  // MPI distribution.
238  ParMesh pmesh(MPI_COMM_WORLD, mesh);
239  mesh.Clear();
240
241  Coefficient *ls_coeff = nullptr;
242  int smooth_steps;
243  if (problem == 0)
244  {
245  ls_coeff = new DeltaCoefficient(0.5, -0.5, 1000.0);
246  smooth_steps = 0;
247  }
248  else if (problem == 1)
249  {
250  ls_coeff = new FunctionCoefficient(sphere_ls);
251  smooth_steps = 0;
252  }
253  else if (problem == 2)
254  {
255  ls_coeff = new FunctionCoefficient(sine_ls);
256  smooth_steps = 0;
257  }
258  else if (problem == 3)
259  {
260  ls_coeff = new FunctionCoefficient(Gyroid);
261  smooth_steps = 0;
262  }
263  else if (problem == 4)
264  {
265  ls_coeff = new FunctionCoefficient(doughnut_cheese);
266  smooth_steps = 0;
267  }
268  else { MFEM_ABORT("Unrecognized -problem option."); }
269
270  const double dx = AvgElementSize(pmesh);
271  DistanceSolver *dist_solver = NULL;
272  if (solver_type == 0)
273  {
274  auto ds = new HeatDistanceSolver(t_param * dx * dx);
275  if (problem == 0)
276  {
277  ds->transform = false;
278  }
279  ds->smooth_steps = smooth_steps;
280  ds->vis_glvis = false;
281  dist_solver = ds;
282  }
283  else if (solver_type == 1)
284  {
285  const int p = 10;
286  const int newton_iter = 50;
287  auto ds = new PLapDistanceSolver(p, newton_iter);
288  dist_solver = ds;
289  }
290  else { MFEM_ABORT("Wrong solver option."); }
291  dist_solver->print_level = 1;
292
293  H1_FECollection fec(order, dim);
294  ParFiniteElementSpace pfes_s(&pmesh, &fec), pfes_v(&pmesh, &fec, dim);
295  ParGridFunction distance_s(&pfes_s), distance_v(&pfes_v);
296
297  // Smooth-out Gibbs oscillations from the input level set. The smoothing
298  // parameter here is specified to be mesh dependent with length scale dx.
299  ParGridFunction filt_gf(&pfes_s);
300  if (problem != 0)
301  {
302  PDEFilter filter(pmesh, 1.0 * dx);
303  filter.Filter(*ls_coeff, filt_gf);
304  }
305  else { filt_gf.ProjectCoefficient(*ls_coeff); }
306  delete ls_coeff;
307  GridFunctionCoefficient ls_filt_coeff(&filt_gf);
308
309  dist_solver->ComputeScalarDistance(ls_filt_coeff, distance_s);
310  dist_solver->ComputeVectorDistance(ls_filt_coeff, distance_v);
311
312  // Send the solution by socket to a GLVis server.
313  if (visualization)
314  {
315  int size = 500;
316  char vishost[] = "localhost";
317  int visport = 19916;
318
319  socketstream sol_sock_w;
320  common::VisualizeField(sol_sock_w, vishost, visport, filt_gf,
321  "Input Level Set", 0, 0, size, size);
322
323  MPI_Barrier(pmesh.GetComm());
324
325  socketstream sol_sock_ds;
326  common::VisualizeField(sol_sock_ds, vishost, visport, distance_s,
327  "Distance", size, 0, size, size,
328  "rRjmm********A");
329
330  MPI_Barrier(pmesh.GetComm());
331
332  socketstream sol_sock_dv;
333  common::VisualizeField(sol_sock_dv, vishost, visport, distance_v,
334  "Directions", 2*size, 0, size, size,
335  "rRjmm********vveA");
336  }
337
338  // ParaView output.
339  ParaViewDataCollection dacol("ParaViewDistance", &pmesh);
340  dacol.SetLevelsOfDetail(order);
341  dacol.RegisterField("filtered_level_set", &filt_gf);
342  dacol.RegisterField("distance", &distance_s);
343  dacol.SetTime(1.0);
344  dacol.SetCycle(1);
345  dacol.Save();
346
347  ConstantCoefficient zero(0.0);
348  const double s_norm = distance_s.ComputeL2Error(zero),
349  v_norm = distance_v.ComputeL2Error(zero);
350  if (myid == 0)
351  {
352  cout << fixed << setprecision(10) << "Norms: "
353  << s_norm << " " << v_norm << std::endl;
354  }
355
356  delete dist_solver;
357
358  return 0;
359 }
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:521
Helper class for ParaView visualization data.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
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:525
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
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
constexpr int visport
double sine_ls(const Vector &x)
Definition: distance.cpp:93
double Sph(const mfem::Vector &xx)
Definition: distance.cpp:136
double sphere_ls(const Vector &x)
Definition: distance.cpp:100
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetTime(double t)
Set physical time (for time-dependent simulations)
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
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
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
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
void DGyroid(const mfem::Vector &xx, mfem::Vector &vals)
Definition: distance.cpp:149
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:281
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
int dim
Definition: ex24.cpp:53
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:904
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
double doughnut_cheese(const Vector &coord)
Definition: sbm_aux.hpp:32
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
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:121
Class for parallel meshes.
Definition: pmesh.hpp:32
int main()
void Filter(ParGridFunction &func, ParGridFunction &ffield)
double AvgElementSize(ParMesh &pmesh)
Definition: dist_solver.cpp:42
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
virtual void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)=0