MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
distance.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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 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>
94#include "sbm_aux.hpp"
95
96using namespace std;
97using namespace mfem;
98using namespace common;
99
101{
102 const real_t 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
108
110{
111 const int dim = x.Size();
112 const real_t xc = x(0) - 0.5;
113 const real_t yc = (dim > 1) ? x(1) - 0.5 : 0.0;
114 const real_t zc = (dim > 2) ? x(2) - 0.5 : 0.0;
115 const real_t r = sqrt(xc*xc + yc*yc + zc*zc);
116
117 return (r >= radius) ? -1.0 : 1.0;
118}
119
121{
122 const int dim = x.Size();
123 const real_t xc = x(0) - 0.5;
124 const real_t yc = (dim > 1) ? x(1) - 0.5 : 0.0;
125 const real_t zc = (dim > 2) ? x(2) - 0.5 : 0.0;
126 const real_t r = sqrt(xc*xc + yc*yc + zc*zc);
127
129}
130
131class ExactDistSphereLoc : public Coefficient
132{
133private:
134 ParGridFunction &dist;
135 const real_t dx;
136
137public:
138 ExactDistSphereLoc(ParGridFunction &d)
139 : dist(d), dx(dist.ParFESpace()->GetParMesh()->GetElementSize(0)) { }
140
141 virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
142 {
143 Vector pos(T.GetDimension());
144 T.Transform(ip, pos);
145 pos -= 0.5;
146 const real_t 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
156{
157 const real_t period = 2.0 * M_PI;
158 real_t x = xx[0]*period;
159 real_t y = xx[1]*period;
160 real_t 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
168{
169 real_t 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
180void DGyroid(const mfem::Vector &xx, mfem::Vector &vals)
181{
182 vals.SetSize(xx.Size());
183 vals=0.0;
184
185 real_t 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
204int 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 real_t t_param = 1.0;
218 const char *device_config = "cpu";
219 bool visualization = true;
220
221 OptionsParser args(argc, argv);
223 "Mesh file to use.");
225 "Solver type:\n\t"
226 "0: Heat\n\t"
227 "1: P-Laplacian\n\t"
228 "2: Rvachev scaling");
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.");
237 "Number of times to refine the mesh uniformly in serial.");
239 "Finite element order (polynomial degree) or -1 for"
240 " isoparametric space.");
242 "Diffusion time step (scaled internally scaled by dx*dx).");
244 "Device configuration string, see Device::Configure().");
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 real_t 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 real_t 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 real_t 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 real_t 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 real_t 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}
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
A coefficient that is constant across space and time.
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
void SetTime(real_t t)
Set physical time (for time-dependent simulations)
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
int GetDimension() const
Return the topological dimension of the reference element.
Definition eltrans.hpp:165
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
Class for integration point with weight.
Definition intrules.hpp:35
Mesh data type.
Definition mesh.hpp:56
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Abstract parallel finite element space.
Definition pfespace.hpp:29
Class for parallel grid function.
Definition pgridfunc.hpp:33
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
real_t GetValue(int i, const IntegrationPoint &ip, int vdim=1) const override
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Class for parallel meshes.
Definition pmesh.hpp:34
MPI_Comm GetComm() const
Definition pmesh.hpp:402
Helper class for ParaView visualization data.
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
Vector data type.
Definition vector.hpp:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
IterativeSolver::PrintLevel print_level
virtual void ComputeScalarDistance(Coefficient &zero_level_set, ParGridFunction &distance)=0
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
void Filter(ParGridFunction &func, ParGridFunction &ffield)
real_t sphere_ls(const Vector &x)
Definition distance.cpp:109
Definition distance.cpp:107
void DGyroid(const mfem::Vector &xx, mfem::Vector &vals)
Definition distance.cpp:180
real_t Sph(const mfem::Vector &xx)
Definition distance.cpp:167
real_t Gyroid(const Vector &xx)
Definition distance.cpp:155
real_t exact_dist_sphere(const Vector &x)
Definition distance.cpp:120
real_t sine_ls(const Vector &x)
Definition distance.cpp:100
int problem
Definition ex15.cpp:62
int dim
Definition ex24.cpp:53
int main()
real_t AvgElementSize(ParMesh &pmesh)
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, const char *keys, bool vec)
const int visport
float real_t
Definition config.hpp:43
const char vishost[]
real_t p(const Vector &x, real_t t)
real_t doughnut_cheese(const Vector &coord)
Definition sbm_aux.hpp:32