MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
diffusion.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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// Shifted Diffusion Miniapp: Finite element immersed boundary solver
14// ------------------------------------------------------------------
15//
16// This miniapp solves the Poisson problem with prescribed boundary conditions
17// on a surrogate domain using a high-order extension of the shifted boundary
18// method [1]. Using a level-set prescribed to represent the true boundary, a
19// distance function is computed for the immersed background mesh. A surrogate
20// domain is also computed to solve the Poisson problem, based on intersection
21// of the zero level-set with the background mesh. The distance function is used
22// with a Taylor expansion to enforce Dirichlet boundary conditions on the
23// (non-aligned) mesh faces of the surrogate domain, therefore "shifting" the
24// location where boundary conditions are imposed.
25//
26// [1] Atallah, Nabil M., Claudio Canuto and Guglielmo Scovazzi. "The
27// second-generation Shifted Boundary Method and its numerical analysis."
28// Computer Methods in Applied Mechanics and Engineering 372 (2020): 113341.
29//
30// Compile with: make diffusion
31//
32// Sample runs:
33//
34// Problem 1: Circular hole of radius 0.2 at the center of the domain.
35// Solves -nabla^2 u = 1 with homogeneous boundary conditions.
36// Dirichlet boundary condition
37// mpirun -np 4 diffusion -rs 3 -o 1 -vis -lst 1
38// mpirun -np 4 diffusion -m ../../data/inline-hex.mesh -rs 2 -o 2 -vis -lst 1 -ho 1 -alpha 10
39// Neumann boundary condition
40// mpirun -np 4 diffusion -rs 3 -o 1 -vis -nlst 1 -ho 1
41//
42// Problem 2: Circular hole of radius 0.2 at the center of the domain.
43// Solves -nabla^2 u = f with inhomogeneous boundary conditions,
44// and f is setup such that u = x^p + y^p, where p = 2 by default.
45// This is a 2D convergence test.
46// Dirichlet BC
47// mpirun -np 4 diffusion -rs 2 -o 2 -vis -lst 2
48// Neumann BC (inhomogeneous condition derived using exact solution)
49// mpirun -np 4 diffusion -rs 2 -o 2 -vis -nlst 2 -ho 1
50//
51// Problem 3: Domain is y = [0, 1] but mesh is shifted to [-1.e-4, 1].
52// Solves -nabla^2 u = f with inhomogeneous boundary conditions,
53// and f is setup such that u = sin(pi*x*y). This is a 2D
54// convergence test. Second-order can be demonstrated by changing
55// refinement level (-rs) for the sample run below. Higher-order
56// convergence can be realized by also increasing the order (-o) of
57// the finite element space and the number of high-order terms
58// (-ho) to be included from the Taylor expansion used to enforce
59// the boundary conditions.
60// mpirun -np 4 diffusion -rs 2 -o 1 -vis -lst 3
61//
62// Problem 4: Complex 2D / 3D shapes:
63// Solves -nabla^2 u = 1 with homogeneous boundary conditions.
64// mpirun -np 4 diffusion -m ../../data/inline-quad.mesh -rs 4 -lst 4 -alpha 10
65// mpirun -np 4 diffusion -m ../../data/inline-tri.mesh -rs 4 -lst 4 -alpha 10
66// mpirun -np 4 diffusion -m ../../data/inline-hex.mesh -rs 3 -lst 8 -alpha 10
67// mpirun -np 4 diffusion -m ../../data/inline-tet.mesh -rs 3 -lst 8 -alpha 10
68//
69// Problem 5: Circular hole with homogeneous Neumann, triangular hole with
70// inhomogeneous Dirichlet, and a square hole with homogeneous
71// Dirichlet boundary condition.
72// mpirun -np 4 diffusion -rs 3 -o 1 -vis -lst 5 -ho 1 -nlst 7 -alpha 10.0 -dc
73
74#include "mfem.hpp"
76#include "sbm_aux.hpp"
77#include "sbm_solver.hpp"
78#include "marking.hpp"
79
80using namespace mfem;
81using namespace std;
82using namespace common;
83
84int main(int argc, char *argv[])
85{
86#ifdef HYPRE_USING_GPU
87 cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this miniapp\n"
88 << "is NOT supported with the GPU version of hypre.\n\n";
89 return MFEM_SKIP_RETURN_VALUE;
90#endif
91
92 // Initialize MPI and HYPRE.
93 Mpi::Init(argc, argv);
94 int myid = Mpi::WorldRank();
96
97 // Parse command-line options.
98 const char *mesh_file = "../../data/inline-quad.mesh";
99 int order = 2;
100 bool visualization = true;
101 int ser_ref_levels = 0;
102 int dirichlet_level_set_type = -1;
103 int neumann_level_set_type = -1;
104 bool dirichlet_combo = false;
105 int ho_terms = 0;
106 real_t alpha = 1;
107 int visport = 19916;
108 bool include_cut_cell = false;
109
110 OptionsParser args(argc, argv);
111 args.AddOption(&mesh_file, "-m", "--mesh",
112 "Mesh file to use.");
113 args.AddOption(&order, "-o", "--order",
114 "Finite element order (polynomial degree) or -1 for"
115 " isoparametric space.");
116 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
117 "--no-visualization",
118 "Enable or disable GLVis visualization.");
119 args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
120 "Number of times to refine the mesh uniformly in serial.");
121 args.AddOption(&dirichlet_level_set_type, "-lst", "--level-set-type",
122 "level-set-type.");
123 args.AddOption(&neumann_level_set_type, "-nlst", "--neumann-level-set-type",
124 "neumann-level-set-type.");
125 args.AddOption(&dirichlet_combo, "-dc", "--dcombo",
126 "no-dc", "--no-dcombo",
127 "Combination of two Dirichlet level sets.");
128 args.AddOption(&ho_terms, "-ho", "--high-order",
129 "Additional high-order terms to include");
130 args.AddOption(&alpha, "-alpha", "--alpha",
131 "Nitsche penalty parameter (~1 for 2D, ~10 for 3D).");
132 args.AddOption(&include_cut_cell, "-cut", "--cut", "-no-cut-cell",
133 "--no-cut-cell",
134 "Include or not include elements cut by true boundary.");
135 args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
136 args.Parse();
137 if (!args.Good())
138 {
139 if (myid == 0)
140 {
141 args.PrintUsage(cout);
142 }
143 return 1;
144 }
145 if (myid == 0) { args.PrintOptions(cout); }
146
147 // Use Dirichlet level set if no level sets are specified.
148 if (dirichlet_level_set_type <= 0 && neumann_level_set_type <= 0)
149 {
150 dirichlet_level_set_type = 1;
151 }
152 MFEM_VERIFY((neumann_level_set_type > 0 && ho_terms < 1) == false,
153 "Shifted Neumann BC requires extra terms, i.e., -ho >= 1.");
154
155 // Enable hardware devices such as GPUs, and programming models such as CUDA,
156 // OCCA, RAJA and OpenMP based on command line options.
157 Device device("cpu");
158 if (myid == 0) { device.Print(); }
159
160 // Refine the mesh.
161 Mesh mesh(mesh_file, 1, 1);
162 int dim = mesh.Dimension();
163 for (int lev = 0; lev < ser_ref_levels; lev++) { mesh.UniformRefinement(); }
164 if (myid == 0)
165 {
166 std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
167 }
168 MFEM_VERIFY(mesh.Conforming(), "AMR capability is not implemented yet!");
169
170 // MPI distribution.
171 ParMesh pmesh(MPI_COMM_WORLD, mesh);
172 mesh.Clear();
173
174 // Define a finite element space on the mesh. Here we use continuous Lagrange
175 // finite elements of the specified order. If order < 1, we fix it to 1.
176 if (order < 1) { order = 1; }
177 H1_FECollection fec(order, dim);
178 ParFiniteElementSpace pfespace(&pmesh, &fec);
179
180 Vector vxyz;
181
182 // Set the nodal grid function for the mesh, and modify the nodal positions
183 // for dirichlet_level_set_type = 3 such that some of the mesh elements are
184 // intersected by the true boundary (y = 0).
185 ParFiniteElementSpace pfespace_mesh(&pmesh, &fec, dim);
186 pmesh.SetNodalFESpace(&pfespace_mesh);
187 ParGridFunction x_mesh(&pfespace_mesh);
188 pmesh.SetNodalGridFunction(&x_mesh);
189 vxyz = *pmesh.GetNodes();
190 int nodes_cnt = vxyz.Size()/dim;
191 if (dirichlet_level_set_type == 3)
192 {
193 for (int i = 0; i < nodes_cnt; i++)
194 {
195 // Shift the mesh from y = [0, 1] to [-1.e-4, 1].
196 vxyz(i+nodes_cnt) = (1.+1.e-4)*vxyz(i+nodes_cnt)-1.e-4;
197 }
198 }
199 pmesh.SetNodes(vxyz);
200 pfespace.ExchangeFaceNbrData();
201 const int gtsize = pfespace.GlobalTrueVSize();
202 if (myid == 0)
203 {
204 cout << "Number of finite element unknowns: " << gtsize << endl;
205 }
206
207 // Define the solution vector x as a finite element grid function
208 // corresponding to pfespace.
209 ParGridFunction x(&pfespace);
210
211 // Determine if each element in the ParMesh is inside the actual domain,
212 // partially cut by its boundary, or completely outside the domain.
213 // Setup the level-set coefficients, and mark the elements.
214 Dist_Level_Set_Coefficient *dirichlet_dist_coef = NULL;
215 Dist_Level_Set_Coefficient *dirichlet_dist_coef_2 = NULL;
216 Dist_Level_Set_Coefficient *neumann_dist_coef = NULL;
217 Combo_Level_Set_Coefficient combo_dist_coef;
218
219 ParGridFunction level_set_gf(&pfespace);
220 ShiftedFaceMarker marker(pmesh, pfespace, include_cut_cell);
221 Array<int> elem_marker;
222
223 // Dirichlet level-set.
224 if (dirichlet_level_set_type > 0)
225 {
226 dirichlet_dist_coef = new Dist_Level_Set_Coefficient(dirichlet_level_set_type);
227 const real_t dx = AvgElementSize(pmesh);
228 PDEFilter filter(pmesh, dx);
229 filter.Filter(*dirichlet_dist_coef, level_set_gf);
230 //level_set_gf.ProjectCoefficient(*dirichlet_dist_coef);
231 // Exchange information for ghost elements i.e. elements that share a face
232 // with element on the current processor, but belong to another processor.
233 level_set_gf.ExchangeFaceNbrData();
234 // Setup the class to mark all elements based on whether they are located
235 // inside or outside the true domain, or intersected by the true boundary.
236 marker.MarkElements(level_set_gf, elem_marker);
237 combo_dist_coef.Add_Level_Set_Coefficient(*dirichlet_dist_coef);
238 }
239
240 // Second Dirichlet level-set.
241 if (dirichlet_combo)
242 {
243 MFEM_VERIFY(dirichlet_level_set_type == 5,
244 "The combo level set example has been only set for"
245 " dirichlet_level_set_type == 5.");
246 dirichlet_dist_coef_2 = new Dist_Level_Set_Coefficient(6);
247 level_set_gf.ProjectCoefficient(*dirichlet_dist_coef_2);
248 level_set_gf.ExchangeFaceNbrData();
249 marker.MarkElements(level_set_gf, elem_marker);
250 combo_dist_coef.Add_Level_Set_Coefficient(*dirichlet_dist_coef_2);
251 }
252
253 // Neumann level-set.
254 if (neumann_level_set_type > 0)
255 {
256 neumann_dist_coef = new Dist_Level_Set_Coefficient(neumann_level_set_type);
257 level_set_gf.ProjectCoefficient(*neumann_dist_coef);
258 level_set_gf.ExchangeFaceNbrData();
259 marker.MarkElements(level_set_gf, elem_marker);
260 combo_dist_coef.Add_Level_Set_Coefficient(*neumann_dist_coef);
261 }
262
263 // Visualize the element markers.
264 if (visualization)
265 {
267 ParFiniteElementSpace pfesl2(&pmesh, &fecl2);
268 ParGridFunction elem_marker_gf(&pfesl2);
269 for (int i = 0; i < elem_marker_gf.Size(); i++)
270 {
271 elem_marker_gf(i) = (real_t)elem_marker[i];
272 }
273 char vishost[] = "localhost";
274 int s = 350;
275 socketstream sol_sock;
276 common::VisualizeField(sol_sock, vishost, visport, elem_marker_gf,
277 "Element Flags", 0, 0, s, s, "Rjmpc");
278 }
279
280 // Get a list of dofs associated with shifted boundary (SB) faces.
281 Array<int> sb_dofs; // Array of dofs on SB faces
282 marker.ListShiftedFaceDofs(elem_marker, sb_dofs);
283
284 // Visualize the shifted boundary face dofs.
285 if (visualization)
286 {
287 ParGridFunction face_dofs(&pfespace);
288 face_dofs = 0.0;
289 for (int i = 0; i < sb_dofs.Size(); i++)
290 {
291 face_dofs(sb_dofs[i]) = 1.0;
292 }
293 char vishost[] = "localhost";
294 int s = 350;
295 socketstream sol_sock;
296 common::VisualizeField(sol_sock, vishost, visport, face_dofs,
297 "Shifted Face Dofs", 0, s, s, s, "Rjmplo");
298 }
299
300 // Make a list of inactive tdofs that will be eliminated from the system.
301 // The inactive tdofs are the dofs for the elements located outside the true
302 // domain (and optionally, for the elements cut by the true boundary, if
303 // include_cut_cell = false) minus the dofs that are located on the surrogate
304 // boundary.
305 Array<int> ess_tdof_list;
306 Array<int> ess_shift_bdr;
307 marker.ListEssentialTDofs(elem_marker, sb_dofs, ess_tdof_list,
308 ess_shift_bdr);
309
310 // Compute distance vector to the actual boundary.
311 ParFiniteElementSpace distance_vec_space(&pmesh, &fec, dim);
312 ParGridFunction distance(&distance_vec_space);
313 VectorCoefficient *dist_vec = NULL;
314 // Compute the distance field analytically or using the HeatDistanceSolver.
315 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type == 2 ||
316 dirichlet_level_set_type == 3)
317 {
318 // Analytic distance vector.
319 dist_vec = new Dist_Vector_Coefficient(dim, dirichlet_level_set_type);
320 distance.ProjectDiscCoefficient(*dist_vec);
321 }
322 else
323 {
324 // Discrete distance vector.
325 real_t dx = AvgElementSize(pmesh);
326 ParGridFunction filt_gf(&pfespace);
327 PDEFilter filter(pmesh, 2.0 * dx);
328 filter.Filter(combo_dist_coef, filt_gf);
329 GridFunctionCoefficient ls_filt_coeff(&filt_gf);
330
331 if (visualization)
332 {
333 char vishost[] = "localhost";
334 int s = 350;
335 socketstream sol_sock;
336 common::VisualizeField(sol_sock, vishost, visport, filt_gf,
337 "Input Level Set", 0, 2*s, s, s, "Rjmm");
338 }
339
340 HeatDistanceSolver dist_func(2.0 * dx * dx);
341 dist_func.print_level.FirstAndLast().Summary();
342 dist_func.smooth_steps = 1;
343 dist_func.ComputeVectorDistance(ls_filt_coeff, distance);
344 dist_vec = new VectorGridFunctionCoefficient(&distance);
345 }
346
347 // Visualize the distance vector.
348 if (visualization)
349 {
350 char vishost[] = "localhost";
351 int s = 350;
352 socketstream sol_sock;
353 common::VisualizeField(sol_sock, vishost, visport, distance,
354 "Distance Vector", s, s, s, s, "Rjmmpcvv", 1);
355 }
356
357 // Set up a list to indicate element attributes to be included in assembly,
358 // so that inactive elements are excluded.
359 const int max_elem_attr = pmesh.attributes.Max();
360 Array<int> ess_elem(max_elem_attr);
361 ess_elem = 1;
362 bool inactive_elements = false;
363 for (int i = 0; i < pmesh.GetNE(); i++)
364 {
365 if (!include_cut_cell &&
366 (elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE ||
367 elem_marker[i] >= ShiftedFaceMarker::SBElementType::CUT))
368 {
369 pmesh.SetAttribute(i, max_elem_attr+1);
370 inactive_elements = true;
371 }
372 if (include_cut_cell &&
373 elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE)
374 {
375 pmesh.SetAttribute(i, max_elem_attr+1);
376 inactive_elements = true;
377 }
378 }
379 bool inactive_elements_global;
380 MPI_Allreduce(&inactive_elements, &inactive_elements_global, 1,
381 MFEM_MPI_CXX_BOOL, MPI_LOR, MPI_COMM_WORLD);
382 if (inactive_elements_global) { ess_elem.Append(0); }
383 pmesh.SetAttributes();
384
385 // Set up the linear form b(.) which corresponds to the right-hand side of
386 // the FEM linear system.
387 ParLinearForm b(&pfespace);
388 FunctionCoefficient *rhs_f = NULL;
389 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type == 4 ||
390 dirichlet_level_set_type == 5 || dirichlet_level_set_type == 6 ||
391 dirichlet_level_set_type == 8 ||
392 neumann_level_set_type == 1 || neumann_level_set_type == 7)
393 {
395 }
396 else if (dirichlet_level_set_type == 2 || neumann_level_set_type == 2)
397 {
399 }
400 else if (dirichlet_level_set_type == 3)
401 {
403 }
404 else { MFEM_ABORT("RHS function not set for level set type.\n"); }
405 b.AddDomainIntegrator(new DomainLFIntegrator(*rhs_f), ess_elem);
406
407 // Exact solution to project for Dirichlet boundaries
408 FunctionCoefficient *exactCoef = NULL;
409 // Dirichlet BC that must be imposed on the true boundary.
410 ShiftedFunctionCoefficient *dbcCoef = NULL;
411 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type >= 4)
412 {
414 exactCoef = new FunctionCoefficient(homogeneous);
415 }
416 else if (dirichlet_level_set_type == 2)
417 {
420 }
421 else if (dirichlet_level_set_type == 3)
422 {
425 }
426
427 ShiftedFunctionCoefficient *dbcCoefCombo = NULL;
428 if (dirichlet_combo)
429 {
430 dbcCoefCombo = new ShiftedFunctionCoefficient(0.015);
431 }
432
433 // Homogeneous Neumann boundary condition coefficient
434 ShiftedFunctionCoefficient *nbcCoef = NULL;
435 ShiftedVectorFunctionCoefficient *normalbcCoef = NULL;
436 if (neumann_level_set_type == 1)
437 {
440 }
441 else if (neumann_level_set_type == 2)
442 {
446 }
447 else if (neumann_level_set_type == 7)
448 {
451 }
452 else if (neumann_level_set_type > 0)
453 {
454 MFEM_ABORT(" Normal vector coefficient not implemented for level set.");
455 }
456
457 // Add integrators corresponding to the shifted boundary method (SBM)
458 // for Dirichlet boundaries.
459 // For each LinearFormIntegrator, we indicate the marker that we have used
460 // for the cut-cell corresponding to the level-set.
461 int ls_cut_marker = ShiftedFaceMarker::SBElementType::CUT;
462 // For each BilinearFormIntegrators, we make a list of the markers
463 // corresponding to the cut-cell whose faces they will be applied to.
464 Array<int> bf_dirichlet_marker(0), bf_neumann_marker(0);
465
466 if (dirichlet_level_set_type > 0)
467 {
468 b.AddInteriorFaceIntegrator(new SBM2DirichletLFIntegrator(&pmesh, *dbcCoef,
469 alpha, *dist_vec,
470 elem_marker,
471 include_cut_cell,
472 ho_terms,
473 ls_cut_marker));
474 b.AddBdrFaceIntegrator(new SBM2DirichletLFIntegrator(&pmesh, *dbcCoef,
475 alpha, *dist_vec,
476 elem_marker,
477 include_cut_cell,
478 ho_terms,
479 ls_cut_marker),
480 ess_shift_bdr);
481 bf_dirichlet_marker.Append(ls_cut_marker);
482 ls_cut_marker += 1;
483 }
484
485 if (dirichlet_combo)
486 {
487 b.AddInteriorFaceIntegrator(new SBM2DirichletLFIntegrator(&pmesh, *dbcCoefCombo,
488 alpha, *dist_vec,
489 elem_marker,
490 include_cut_cell,
491 ho_terms,
492 ls_cut_marker));
493 bf_dirichlet_marker.Append(ls_cut_marker);
494 ls_cut_marker += 1;
495 }
496
497 // Add integrators corresponding to the shifted boundary method (SBM)
498 // for Neumann boundaries.
499 if (neumann_level_set_type > 0)
500 {
501 MFEM_VERIFY(!include_cut_cell, "include_cut_cell option must be set to"
502 " false for Neumann boundary conditions.");
503 b.AddInteriorFaceIntegrator(new SBM2NeumannLFIntegrator(
504 &pmesh, *nbcCoef, *dist_vec,
505 *normalbcCoef, elem_marker,
506 ho_terms, include_cut_cell,
507 ls_cut_marker));
508 bf_neumann_marker.Append(ls_cut_marker);
509 ls_cut_marker += 1;
510 }
511
512 b.Assemble();
513
514 // Set up the bilinear form a(.,.) on the finite element space corresponding
515 // to the Laplacian operator -Delta, by adding the Diffusion domain
516 // integrator and SBM integrator.
517 ParBilinearForm a(&pfespace);
518 ConstantCoefficient one(1.);
519 a.AddDomainIntegrator(new DiffusionIntegrator(one), ess_elem);
520 if (dirichlet_level_set_type > 0)
521 {
522 a.AddInteriorFaceIntegrator(new SBM2DirichletIntegrator(&pmesh, alpha,
523 *dist_vec,
524 elem_marker,
525 bf_dirichlet_marker,
526 include_cut_cell,
527 ho_terms));
528 a.AddBdrFaceIntegrator(new SBM2DirichletIntegrator(&pmesh, alpha, *dist_vec,
529 elem_marker,
530 bf_dirichlet_marker,
531 include_cut_cell,
532 ho_terms), ess_shift_bdr);
533 }
534
535 // Add Neumann bilinear form integrator.
536 if (neumann_level_set_type > 0)
537 {
538 a.AddInteriorFaceIntegrator(new SBM2NeumannIntegrator(&pmesh,
539 *dist_vec,
540 *normalbcCoef,
541 elem_marker,
542 bf_neumann_marker,
543 include_cut_cell,
544 ho_terms));
545 }
546
547 // Assemble the bilinear form and the corresponding linear system,
548 // applying any necessary transformations.
549 a.Assemble();
550
551 // Project the exact solution as an initial condition for Dirichlet boundary.
552 if (!exactCoef)
553 {
554 exactCoef = new FunctionCoefficient(homogeneous);
555 }
556 x.ProjectCoefficient(*exactCoef);
557
558 // Form the linear system and solve it.
559 OperatorPtr A;
560 Vector B, X;
561 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
562
563 Solver *prec = new HypreBoomerAMG;
564 BiCGSTABSolver bicg(MPI_COMM_WORLD);
565 bicg.SetRelTol(1e-12);
566 bicg.SetMaxIter(500);
567 bicg.SetPrintLevel(1);
568 bicg.SetPreconditioner(*prec);
569 bicg.SetOperator(*A);
570 bicg.Mult(B, X);
571
572 // Recover the solution as a finite element grid function.
573 a.RecoverFEMSolution(X, b, x);
574
575 // Save the mesh and the solution.
576 ofstream mesh_ofs("diffusion.mesh");
577 mesh_ofs.precision(8);
578 pmesh.PrintAsOne(mesh_ofs);
579 ofstream sol_ofs("diffusion.gf");
580 sol_ofs.precision(8);
581 x.SaveAsOne(sol_ofs);
582
583 if (visualization)
584 {
585 // Save the solution in ParaView format.
586 ParaViewDataCollection dacol("ParaViewDiffusion", &pmesh);
587 dacol.SetLevelsOfDetail(order);
588 dacol.RegisterField("distance", &distance);
589 dacol.RegisterField("level_set", &level_set_gf);
590 dacol.RegisterField("solution", &x);
591 dacol.SetTime(1.0);
592 dacol.SetCycle(1);
593 dacol.Save();
594
595 // Send the solution by socket to a GLVis server.
596 char vishost[] = "localhost";
597 int s = 350;
598 socketstream sol_sock;
599 common::VisualizeField(sol_sock, vishost, visport, x,
600 "Solution", s, 0, s, s, "Rj");
601 }
602
603 // Construct an error grid function if the exact solution is known.
604 if (dirichlet_level_set_type == 2 || dirichlet_level_set_type == 3 ||
605 (dirichlet_level_set_type == -1 && neumann_level_set_type == 2))
606 {
607 ParGridFunction error(x);
608 Vector pxyz(dim);
609 pxyz(0) = 0.;
610 for (int i = 0; i < nodes_cnt; i++)
611 {
612 pxyz(0) = vxyz(i);
613 pxyz(1) = vxyz(i+nodes_cnt);
614 real_t exact_val = 0.;
615 if (dirichlet_level_set_type == 2 || neumann_level_set_type == 2)
616 {
617 exact_val = dirichlet_velocity_xy_exponent(pxyz);
618 }
619 else if (dirichlet_level_set_type == 3)
620 {
621 exact_val = dirichlet_velocity_xy_sinusoidal(pxyz);
622 }
623 error(i) = std::fabs(x(i) - exact_val);
624 }
625
626 if (visualization)
627 {
628 char vishost[] = "localhost";
629 int s = 350;
630 socketstream sol_sock;
631 common::VisualizeField(sol_sock, vishost, visport, error,
632 "Error", 2*s, 0, s, s, "Rj");
633 }
634
635 const real_t global_error = x.ComputeL2Error(*exactCoef);
636 if (myid == 0)
637 {
638 std::cout << "Global L2 error: " << global_error << endl;
639 }
640 }
641
642 const real_t norm = x.ComputeL1Error(one);
643 if (myid == 0) { std::cout << setprecision(10) << norm << std::endl; }
644
645 // Free the used memory.
646 delete prec;
647 delete normalbcCoef;
648 delete nbcCoef;
649 delete dbcCoefCombo;
650 delete dbcCoef;
651 delete exactCoef;
652 delete rhs_f;
653 delete dist_vec;
654 delete neumann_dist_coef;
655 delete dirichlet_dist_coef;
656 delete dirichlet_dist_coef_2;
657
658 return 0;
659}
Combination of level sets: +1 inside the true domain, -1 outside.
Definition sbm_aux.hpp:157
void Add_Level_Set_Coefficient(Dist_Level_Set_Coefficient &dls_)
Definition sbm_aux.hpp:164
Level set coefficient: +1 inside the true domain, -1 outside.
Definition sbm_aux.hpp:138
Distance vector to the zero level-set.
Definition sbm_aux.hpp:184
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
int Size() const
Return the logical size of the array.
Definition array.hpp:147
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
BiCGSTAB method.
Definition solvers.hpp:621
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the BiCGSTAB method.
Definition solvers.cpp:1456
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Definition solvers.hpp:634
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)
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:297
Class for domain integration .
Definition lininteg.hpp:106
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:275
The BoomerAMG solver in hypre.
Definition hypre.hpp:1717
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.cpp:33
void SetRelTol(real_t rtol)
Definition solvers.hpp:229
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:174
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:72
void SetMaxIter(int max_it)
Definition solvers.hpp:231
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:346
Mesh data type.
Definition mesh.hpp:64
bool Conforming() const
Return a bool indicating whether this mesh is conforming.
Definition mesh.hpp:2365
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition mesh.cpp:7629
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:761
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition mesh.cpp:6473
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:9294
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
Definition mesh.cpp:9306
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:11295
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:288
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).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
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.
Class for parallel bilinear form.
Abstract parallel finite element space.
Definition pfespace.hpp:29
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:346
Class for parallel grid function.
Definition pgridfunc.hpp:50
real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const override
Returns ||u_ex - u_h||_L2 in parallel for H1 or L2 elements.
void ProjectDiscCoefficient(VectorCoefficient &coeff) override
Project a discontinuous vector coefficient as a grid function on a continuous finite element space....
MFEM_DEPRECATED real_t ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const override
Returns ||u_ex - u_h||_L1 in parallel for H1 or L2 elements.
void ProjectCoefficient(Coefficient &coeff) override
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void SaveAsOne(const char *fname, int precision=16) const
Class for parallel linear form.
Class for parallel meshes.
Definition pmesh.hpp:34
void SetAttributes() override
Determine the sets of unique attribute values in domain and boundary elements.
Definition pmesh.cpp:1588
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition pmesh.cpp:2027
void PrintAsOne(std::ostream &out=mfem::out, const std::string &comments="") const
Definition pmesh.cpp:4967
Helper class for ParaView visualization data.
void SetLevelsOfDetail(int levels_of_detail_)
void ListShiftedFaceDofs(const Array< int > &elem_marker, Array< int > &sface_dof_list) const
Definition marking.cpp:106
void MarkElements(const ParGridFunction &ls_func, Array< int > &elem_marker)
Definition marking.cpp:17
void ListEssentialTDofs(const Array< int > &elem_marker, const Array< int > &sface_dof_list, Array< int > &ess_tdof_list, Array< int > &ess_shift_bdr) const
Definition marking.cpp:202
Base class for solvers.
Definition operator.hpp:780
Base class for vector Coefficients that optionally depend on time and space.
Vector coefficient defined by a vector GridFunction.
Vector data type.
Definition vector.hpp:82
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
IterativeSolver::PrintLevel print_level
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
void Filter(ParGridFunction &func, ParGridFunction &ffield)
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
int main()
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
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)
float real_t
Definition config.hpp:43
const char vishost[]
real_t rhs_fun_xy_exponent(const Vector &x)
Definition sbm_aux.hpp:273
real_t dirichlet_velocity_xy_exponent(const Vector &x)
Definition sbm_aux.hpp:223
real_t dirichlet_velocity_xy_sinusoidal(const Vector &x)
Definition sbm_aux.hpp:229
real_t rhs_fun_xy_sinusoidal(const Vector &x)
Definition sbm_aux.hpp:288
real_t homogeneous(const Vector &x)
Boundary conditions - Dirichlet.
Definition sbm_aux.hpp:218
real_t traction_xy_exponent(const Vector &x)
Neumann condition for exponent based solution.
Definition sbm_aux.hpp:256
void normal_vector_1(const Vector &x, Vector &p)
Definition sbm_aux.hpp:236
real_t rhs_fun_circle(const Vector &x)
f for the Poisson problem (-nabla^2 u = f).
Definition sbm_aux.hpp:268
void normal_vector_2(const Vector &x, Vector &p)
Normal vector for level_set_type = 7. Circle centered at [0.5 , 0.6].
Definition sbm_aux.hpp:246