MFEM  v4.3.0 Finite element discretization library
distance.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 //
76 // Problem 3: level set: Gyroid
77 // mpirun -np 4 distance -m ../../data/periodic-square.mesh -rs 5 -o 2 -t 1.0 -p 3
78 // mpirun -np 4 distance -m ../../data/periodic-cube.mesh -rs 3 -o 2 -t 1.0 -p 3
79
80 #include <fstream>
81 #include <iostream>
82 #include "dist_solver.hpp"
83 #include "../common/mfem-common.hpp"
84
85 using namespace std;
86 using namespace mfem;
87
88 double sine_ls(const Vector &x)
89 {
90  const double sine = 0.25 * std::sin(4 * M_PI * x(0)) +
91  0.05 * std::sin(16 * M_PI * x(0));
92  return (x(1) >= sine + 0.5) ? -1.0 : 1.0;
93 }
94
95 double sphere_ls(const Vector &x)
96 {
97  const int dim = x.Size();
98  if (dim == 2)
99  {
100  const double xc = x(0) - 0.5, yc = x(1) - 0.5;
101  const double r = sqrt(xc*xc + yc*yc);
102  return (r >= 0.4) ? -1.0 : 1.0;
103  }
104  else if (dim == 3)
105  {
106  const double xc = x(0) - 0.0, yc = x(1) - 0.0, zc = x(2) - 0.0;
107  const double r = sqrt(xc*xc + yc*yc + zc*zc);
108  return (r >= 0.8) ? -1.0 : 1.0;
109  }
110  else
111  {
112  return (x(0) >= 0.5) ? -1.0 : 1.0;
113  }
114 }
115
116 double Gyroid(const Vector &xx)
117 {
118  const double period = 2.0 * M_PI;
119  double x=xx[0]*period;
120  double y=xx[1]*period;
121  double z=0.0;
122  if (xx.Size()==3)
123  {
124  z=xx[2]*period;
125  }
126  return std::sin(x)*std::cos(y) +
127  std::sin(y)*std::cos(z) +
128  std::sin(z)*std::cos(x);
129 }
130
131 double Sph(const mfem::Vector &xx)
132 {
133  double R=0.4;
134  mfem::Vector lvec(3);
135  lvec=0.0;
136  for (int i=0; i<xx.Size(); i++)
137  {
138  lvec[i]=xx[i];
139  }
140
141  return lvec[0]*lvec[0]+lvec[1]*lvec[1]+lvec[2]*lvec[2]-R*R;
142 }
143
144 void DGyroid(const mfem::Vector &xx, mfem::Vector &vals)
145 {
146  vals.SetSize(xx.Size());
147  vals=0.0;
148
149  double pp=4*M_PI;
150
151  mfem::Vector lvec(3);
152  lvec=0.0;
153  for (int i=0; i<xx.Size(); i++)
154  {
155  lvec[i]=xx[i]*pp;
156  }
157
158  vals[0]=cos(lvec[0])*cos(lvec[1])-sin(lvec[2])*sin(lvec[0]);
159  vals[1]=-sin(lvec[0])*sin(lvec[1])+cos(lvec[1])*cos(lvec[2]);
160  if (xx.Size()>2)
161  {
162  vals[2]=-sin(lvec[1])*sin(lvec[2])+cos(lvec[2])*cos(lvec[0]);
163  }
164
165  vals*=pp;
166 }
167
168 int main(int argc, char *argv[])
169 {
170  // Initialize MPI.
171  int num_procs, myid;
172  MPI_Init(&argc, &argv);
173  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
174  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
175
176  // Parse command-line options.
177  const char *mesh_file = "../../data/inline-quad.mesh";
178  int solver_type = 0;
179  int problem = 1;
180  int rs_levels = 2;
181  int order = 2;
182  double t_param = 1.0;
183  const char *device_config = "cpu";
184  bool visualization = true;
185
186  OptionsParser args(argc, argv);
188  "Mesh file to use.");
190  "Solver type:\n\t"
191  "0: Heat\n\t"
192  "1: P-Laplacian");
194  "Problem type:\n\t"
195  "0: Point source\n\t"
196  "1: Circle / sphere level set in 2D / 3D\n\t"
197  "2: 2D sine-looking level set\n\t"
198  "3: Gyroid level set in 2D or 3D");
200  "Number of times to refine the mesh uniformly in serial.");
202  "Finite element order (polynomial degree) or -1 for"
203  " isoparametric space.");
205  "Diffusion time step (scaled internally scaled by dx*dx).");
207  "Device configuration string, see Device::Configure().");
209  "--no-visualization",
210  "Enable or disable GLVis visualization.");
211  args.Parse();
212  if (!args.Good())
213  {
214  if (myid == 0)
215  {
216  args.PrintUsage(cout);
217  }
218  MPI_Finalize();
219  return 1;
220  }
221  if (myid == 0) { args.PrintOptions(cout); }
222
223  // Enable hardware devices such as GPUs, and programming models such as CUDA,
224  // OCCA, RAJA and OpenMP based on command line options.
225  Device device(device_config);
226  if (myid == 0) { device.Print(); }
227
228  // Refine the mesh.
229  Mesh mesh(mesh_file, 1, 1);
230  const int dim = mesh.Dimension();
231  for (int lev = 0; lev < rs_levels; lev++) { mesh.UniformRefinement(); }
232
233  // MPI distribution.
234  ParMesh pmesh(MPI_COMM_WORLD, mesh);
235  mesh.Clear();
236
237  Coefficient *ls_coeff;
238  int smooth_steps;
239  if (problem == 0)
240  {
241  ls_coeff = new DeltaCoefficient(0.5, -0.5, 1000.0);
242  smooth_steps = 0;
243  }
244  else if (problem == 1)
245  {
246  ls_coeff = new FunctionCoefficient(sphere_ls);
247  smooth_steps = 0;
248  }
249  else if (problem == 2)
250  {
251  ls_coeff = new FunctionCoefficient(sine_ls);
252  smooth_steps = 0;
253  }
254  else
255  {
256  ls_coeff = new FunctionCoefficient(Gyroid);
257  smooth_steps = 0;
258  }
259
260  const double dx = AvgElementSize(pmesh);
261  DistanceSolver *dist_solver = NULL;
262  if (solver_type == 0)
263  {
264  auto ds = new HeatDistanceSolver(t_param * dx * dx);
265  if (problem == 0)
266  {
267  ds->transform = false;
268  }
269  ds->smooth_steps = smooth_steps;
270  ds->vis_glvis = false;
271  dist_solver = ds;
272  }
273  else if (solver_type == 1)
274  {
275  const int p = 10;
276  const int newton_iter = 50;
277  auto ds = new PLapDistanceSolver(p, newton_iter);
278  dist_solver = ds;
279  }
280  else { MFEM_ABORT("Wrong solver option."); }
281  dist_solver->print_level = 1;
282
283  H1_FECollection fec(order, dim);
284  ParFiniteElementSpace pfes_s(&pmesh, &fec), pfes_v(&pmesh, &fec, dim);
285  ParGridFunction distance_s(&pfes_s), distance_v(&pfes_v);
286
287  // Smooth-out Gibbs oscillations from the input level set. The smoothing
288  // parameter here is specified to be mesh dependent with length scale dx.
289  ParGridFunction filt_gf(&pfes_s);
290  PDEFilter *filter = new PDEFilter(pmesh, 1.0 * dx);
291  if (problem != 0)
292  {
293  filter->Filter(*ls_coeff, filt_gf);
294  }
295  else { filt_gf.ProjectCoefficient(*ls_coeff); }
296  delete filter;
297  delete ls_coeff;
298  GridFunctionCoefficient ls_filt_coeff(&filt_gf);
299
300  dist_solver->ComputeScalarDistance(ls_filt_coeff, distance_s);
301  dist_solver->ComputeVectorDistance(ls_filt_coeff, distance_v);
302
303  // Send the solution by socket to a GLVis server.
304  if (visualization)
305  {
306  int size = 500;
307  char vishost[] = "localhost";
308  int visport = 19916;
309
310  socketstream sol_sock_w;
311  common::VisualizeField(sol_sock_w, vishost, visport, filt_gf,
312  "Input Level Set", 0, 0, size, size);
313
314  MPI_Barrier(pmesh.GetComm());
315
316  socketstream sol_sock_ds;
317  common::VisualizeField(sol_sock_ds, vishost, visport, distance_s,
318  "Distance", size, 0, size, size,
319  "rRjmm********A");
320
321  MPI_Barrier(pmesh.GetComm());
322
323  socketstream sol_sock_dv;
324  common::VisualizeField(sol_sock_dv, vishost, visport, distance_v,
325  "Directions", 2*size, 0, size, size,
326  "rRjmm********vveA");
327  }
328
329  // ParaView output.
330  ParaViewDataCollection dacol("ParaViewDistance", &pmesh);
331  dacol.SetLevelsOfDetail(order);
332  dacol.RegisterField("filtered_level_set", &filt_gf);
333  dacol.RegisterField("distance", &distance_s);
334  dacol.SetTime(1.0);
335  dacol.SetCycle(1);
336  dacol.Save();
337
338  ConstantCoefficient zero(0.0);
339  const double d_norm = distance_s.ComputeL2Error(zero);
340  if (myid == 0)
341  {
342  cout << fixed << setprecision(10) << "Norm: " << d_norm << std::endl;
343  }
344
345  delete dist_solver;
346
347  MPI_Finalize();
348  return 0;
349 }
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:513
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:190
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:493
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:150
constexpr char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
double sine_ls(const Vector &x)
Definition: distance.cpp:88
double Sph(const mfem::Vector &xx)
Definition: distance.cpp:131
double sphere_ls(const Vector &x)
Definition: distance.cpp:95
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
void SetTime(double t)
Set physical time (for time-dependent simulations)
MPI_Comm GetComm() const
Definition: pmesh.hpp:276
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:144
int dim
Definition: ex24.cpp:53
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:828
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
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:116
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