MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
diffusion.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// 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 bool include_cut_cell = false;
108
109 OptionsParser args(argc, argv);
111 "Mesh file to use.");
113 "Finite element order (polynomial degree) or -1 for"
114 " isoparametric space.");
116 "--no-visualization",
117 "Enable or disable GLVis visualization.");
119 "Number of times to refine the mesh uniformly in serial.");
121 "level-set-type.");
123 "neumann-level-set-type.");
125 "no-dc", "--no-dcombo",
126 "Combination of two Dirichlet level sets.");
128 "Additional high-order terms to include");
130 "Nitsche penalty parameter (~1 for 2D, ~10 for 3D).");
132 "--no-cut-cell",
133 "Include or not include elements cut by true boundary.");
134 args.Parse();
135 if (!args.Good())
136 {
137 if (myid == 0)
138 {
139 args.PrintUsage(cout);
140 }
141 return 1;
142 }
143 if (myid == 0) { args.PrintOptions(cout); }
144
145 // Use Dirichlet level set if no level sets are specified.
146 if (dirichlet_level_set_type <= 0 && neumann_level_set_type <= 0)
147 {
148 dirichlet_level_set_type = 1;
149 }
150 MFEM_VERIFY((neumann_level_set_type > 0 && ho_terms < 1) == false,
151 "Shifted Neumann BC requires extra terms, i.e., -ho >= 1.");
152
153 // Enable hardware devices such as GPUs, and programming models such as CUDA,
154 // OCCA, RAJA and OpenMP based on command line options.
155 Device device("cpu");
156 if (myid == 0) { device.Print(); }
157
158 // Refine the mesh.
159 Mesh mesh(mesh_file, 1, 1);
160 int dim = mesh.Dimension();
161 for (int lev = 0; lev < ser_ref_levels; lev++) { mesh.UniformRefinement(); }
162 if (myid == 0)
163 {
164 std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
165 }
166 MFEM_VERIFY(mesh.Conforming(), "AMR capability is not implemented yet!");
167
168 // MPI distribution.
169 ParMesh pmesh(MPI_COMM_WORLD, mesh);
170 mesh.Clear();
171
172 // Define a finite element space on the mesh. Here we use continuous Lagrange
173 // finite elements of the specified order. If order < 1, we fix it to 1.
174 if (order < 1) { order = 1; }
175 H1_FECollection fec(order, dim);
176 ParFiniteElementSpace pfespace(&pmesh, &fec);
177
178 Vector vxyz;
179
180 // Set the nodal grid function for the mesh, and modify the nodal positions
181 // for dirichlet_level_set_type = 3 such that some of the mesh elements are
182 // intersected by the true boundary (y = 0).
183 ParFiniteElementSpace pfespace_mesh(&pmesh, &fec, dim);
184 pmesh.SetNodalFESpace(&pfespace_mesh);
185 ParGridFunction x_mesh(&pfespace_mesh);
186 pmesh.SetNodalGridFunction(&x_mesh);
187 vxyz = *pmesh.GetNodes();
188 int nodes_cnt = vxyz.Size()/dim;
189 if (dirichlet_level_set_type == 3)
190 {
191 for (int i = 0; i < nodes_cnt; i++)
192 {
193 // Shift the mesh from y = [0, 1] to [-1.e-4, 1].
194 vxyz(i+nodes_cnt) = (1.+1.e-4)*vxyz(i+nodes_cnt)-1.e-4;
195 }
196 }
197 pmesh.SetNodes(vxyz);
198 pfespace.ExchangeFaceNbrData();
199 const int gtsize = pfespace.GlobalTrueVSize();
200 if (myid == 0)
201 {
202 cout << "Number of finite element unknowns: " << gtsize << endl;
203 }
204
205 // Define the solution vector x as a finite element grid function
206 // corresponding to pfespace.
207 ParGridFunction x(&pfespace);
208
209 // Determine if each element in the ParMesh is inside the actual domain,
210 // partially cut by its boundary, or completely outside the domain.
211 // Setup the level-set coefficients, and mark the elements.
212 Dist_Level_Set_Coefficient *dirichlet_dist_coef = NULL;
213 Dist_Level_Set_Coefficient *dirichlet_dist_coef_2 = NULL;
214 Dist_Level_Set_Coefficient *neumann_dist_coef = NULL;
215 Combo_Level_Set_Coefficient combo_dist_coef;
216
217 ParGridFunction level_set_gf(&pfespace);
218 ShiftedFaceMarker marker(pmesh, pfespace, include_cut_cell);
219 Array<int> elem_marker;
220
221 // Dirichlet level-set.
222 if (dirichlet_level_set_type > 0)
223 {
224 dirichlet_dist_coef = new Dist_Level_Set_Coefficient(dirichlet_level_set_type);
225 const real_t dx = AvgElementSize(pmesh);
226 PDEFilter filter(pmesh, dx);
227 filter.Filter(*dirichlet_dist_coef, level_set_gf);
228 //level_set_gf.ProjectCoefficient(*dirichlet_dist_coef);
229 // Exchange information for ghost elements i.e. elements that share a face
230 // with element on the current processor, but belong to another processor.
231 level_set_gf.ExchangeFaceNbrData();
232 // Setup the class to mark all elements based on whether they are located
233 // inside or outside the true domain, or intersected by the true boundary.
234 marker.MarkElements(level_set_gf, elem_marker);
236 }
237
238 // Second Dirichlet level-set.
239 if (dirichlet_combo)
240 {
241 MFEM_VERIFY(dirichlet_level_set_type == 5,
242 "The combo level set example has been only set for"
243 " dirichlet_level_set_type == 5.");
244 dirichlet_dist_coef_2 = new Dist_Level_Set_Coefficient(6);
245 level_set_gf.ProjectCoefficient(*dirichlet_dist_coef_2);
246 level_set_gf.ExchangeFaceNbrData();
247 marker.MarkElements(level_set_gf, elem_marker);
249 }
250
251 // Neumann level-set.
252 if (neumann_level_set_type > 0)
253 {
254 neumann_dist_coef = new Dist_Level_Set_Coefficient(neumann_level_set_type);
255 level_set_gf.ProjectCoefficient(*neumann_dist_coef);
256 level_set_gf.ExchangeFaceNbrData();
257 marker.MarkElements(level_set_gf, elem_marker);
259 }
260
261 // Visualize the element markers.
262 if (visualization)
263 {
265 ParFiniteElementSpace pfesl2(&pmesh, &fecl2);
266 ParGridFunction elem_marker_gf(&pfesl2);
267 for (int i = 0; i < elem_marker_gf.Size(); i++)
268 {
269 elem_marker_gf(i) = (real_t)elem_marker[i];
270 }
271 char vishost[] = "localhost";
272 int visport = 19916, s = 350;
273 socketstream sol_sock;
274 common::VisualizeField(sol_sock, vishost, visport, elem_marker_gf,
275 "Element Flags", 0, 0, s, s, "Rjmpc");
276 }
277
278 // Get a list of dofs associated with shifted boundary (SB) faces.
279 Array<int> sb_dofs; // Array of dofs on SB faces
280 marker.ListShiftedFaceDofs(elem_marker, sb_dofs);
281
282 // Visualize the shifted boundary face dofs.
283 if (visualization)
284 {
285 ParGridFunction face_dofs(&pfespace);
286 face_dofs = 0.0;
287 for (int i = 0; i < sb_dofs.Size(); i++)
288 {
289 face_dofs(sb_dofs[i]) = 1.0;
290 }
291 char vishost[] = "localhost";
292 int visport = 19916, s = 350;
293 socketstream sol_sock;
294 common::VisualizeField(sol_sock, vishost, visport, face_dofs,
295 "Shifted Face Dofs", 0, s, s, s, "Rjmplo");
296 }
297
298 // Make a list of inactive tdofs that will be eliminated from the system.
299 // The inactive tdofs are the dofs for the elements located outside the true
300 // domain (and optionally, for the elements cut by the true boundary, if
301 // include_cut_cell = false) minus the dofs that are located on the surrogate
302 // boundary.
303 Array<int> ess_tdof_list;
304 Array<int> ess_shift_bdr;
305 marker.ListEssentialTDofs(elem_marker, sb_dofs, ess_tdof_list,
306 ess_shift_bdr);
307
308 // Compute distance vector to the actual boundary.
309 ParFiniteElementSpace distance_vec_space(&pmesh, &fec, dim);
310 ParGridFunction distance(&distance_vec_space);
311 VectorCoefficient *dist_vec = NULL;
312 // Compute the distance field analytically or using the HeatDistanceSolver.
313 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type == 2 ||
314 dirichlet_level_set_type == 3)
315 {
316 // Analytic distance vector.
317 dist_vec = new Dist_Vector_Coefficient(dim, dirichlet_level_set_type);
318 distance.ProjectDiscCoefficient(*dist_vec);
319 }
320 else
321 {
322 // Discrete distance vector.
323 real_t dx = AvgElementSize(pmesh);
324 ParGridFunction filt_gf(&pfespace);
325 PDEFilter filter(pmesh, 2.0 * dx);
326 filter.Filter(combo_dist_coef, filt_gf);
327 GridFunctionCoefficient ls_filt_coeff(&filt_gf);
328
329 if (visualization)
330 {
331 char vishost[] = "localhost";
332 int visport = 19916, s = 350;
333 socketstream sol_sock;
334 common::VisualizeField(sol_sock, vishost, visport, filt_gf,
335 "Input Level Set", 0, 2*s, s, s, "Rjmm");
336 }
337
338 HeatDistanceSolver dist_func(2.0 * dx * dx);
339 dist_func.print_level.FirstAndLast().Summary();
340 dist_func.smooth_steps = 1;
341 dist_func.ComputeVectorDistance(ls_filt_coeff, distance);
342 dist_vec = new VectorGridFunctionCoefficient(&distance);
343 }
344
345 // Visualize the distance vector.
346 if (visualization)
347 {
348 char vishost[] = "localhost";
349 int visport = 19916, s = 350;
350 socketstream sol_sock;
351 common::VisualizeField(sol_sock, vishost, visport, distance,
352 "Distance Vector", s, s, s, s, "Rjmmpcvv", 1);
353 }
354
355 // Set up a list to indicate element attributes to be included in assembly,
356 // so that inactive elements are excluded.
357 const int max_elem_attr = pmesh.attributes.Max();
358 Array<int> ess_elem(max_elem_attr);
359 ess_elem = 1;
360 bool inactive_elements = false;
361 for (int i = 0; i < pmesh.GetNE(); i++)
362 {
363 if (!include_cut_cell &&
364 (elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE ||
365 elem_marker[i] >= ShiftedFaceMarker::SBElementType::CUT))
366 {
367 pmesh.SetAttribute(i, max_elem_attr+1);
368 inactive_elements = true;
369 }
370 if (include_cut_cell &&
371 elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE)
372 {
373 pmesh.SetAttribute(i, max_elem_attr+1);
374 inactive_elements = true;
375 }
376 }
377 bool inactive_elements_global;
378 MPI_Allreduce(&inactive_elements, &inactive_elements_global, 1, MPI_C_BOOL,
379 MPI_LOR, MPI_COMM_WORLD);
380 if (inactive_elements_global) { ess_elem.Append(0); }
381 pmesh.SetAttributes();
382
383 // Set up the linear form b(.) which corresponds to the right-hand side of
384 // the FEM linear system.
385 ParLinearForm b(&pfespace);
386 FunctionCoefficient *rhs_f = NULL;
387 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type == 4 ||
388 dirichlet_level_set_type == 5 || dirichlet_level_set_type == 6 ||
389 dirichlet_level_set_type == 8 ||
390 neumann_level_set_type == 1 || neumann_level_set_type == 7)
391 {
393 }
394 else if (dirichlet_level_set_type == 2 || neumann_level_set_type == 2)
395 {
397 }
398 else if (dirichlet_level_set_type == 3)
399 {
401 }
402 else { MFEM_ABORT("RHS function not set for level set type.\n"); }
404
405 // Exact solution to project for Dirichlet boundaries
406 FunctionCoefficient *exactCoef = NULL;
407 // Dirichlet BC that must be imposed on the true boundary.
408 ShiftedFunctionCoefficient *dbcCoef = NULL;
409 if (dirichlet_level_set_type == 1 || dirichlet_level_set_type >= 4)
410 {
412 exactCoef = new FunctionCoefficient(homogeneous);
413 }
414 else if (dirichlet_level_set_type == 2)
415 {
418 }
419 else if (dirichlet_level_set_type == 3)
420 {
423 }
424
425 ShiftedFunctionCoefficient *dbcCoefCombo = NULL;
426 if (dirichlet_combo)
427 {
428 dbcCoefCombo = new ShiftedFunctionCoefficient(0.015);
429 }
430
431 // Homogeneous Neumann boundary condition coefficient
432 ShiftedFunctionCoefficient *nbcCoef = NULL;
433 ShiftedVectorFunctionCoefficient *normalbcCoef = NULL;
434 if (neumann_level_set_type == 1)
435 {
438 }
439 else if (neumann_level_set_type == 2)
440 {
444 }
445 else if (neumann_level_set_type == 7)
446 {
449 }
450 else if (neumann_level_set_type > 0)
451 {
452 MFEM_ABORT(" Normal vector coefficient not implemented for level set.");
453 }
454
455 // Add integrators corresponding to the shifted boundary method (SBM)
456 // for Dirichlet boundaries.
457 // For each LinearFormIntegrator, we indicate the marker that we have used
458 // for the cut-cell corresponding to the level-set.
459 int ls_cut_marker = ShiftedFaceMarker::SBElementType::CUT;
460 // For each BilinearFormIntegrators, we make a list of the markers
461 // corresponding to the cut-cell whose faces they will be applied to.
462 Array<int> bf_dirichlet_marker(0), bf_neumann_marker(0);
463
464 if (dirichlet_level_set_type > 0)
465 {
467 alpha, *dist_vec,
468 elem_marker,
469 include_cut_cell,
470 ho_terms,
471 ls_cut_marker));
473 alpha, *dist_vec,
474 elem_marker,
475 include_cut_cell,
476 ho_terms,
477 ls_cut_marker),
478 ess_shift_bdr);
479 bf_dirichlet_marker.Append(ls_cut_marker);
480 ls_cut_marker += 1;
481 }
482
483 if (dirichlet_combo)
484 {
486 alpha, *dist_vec,
487 elem_marker,
488 include_cut_cell,
489 ho_terms,
490 ls_cut_marker));
491 bf_dirichlet_marker.Append(ls_cut_marker);
492 ls_cut_marker += 1;
493 }
494
495 // Add integrators corresponding to the shifted boundary method (SBM)
496 // for Neumann boundaries.
497 if (neumann_level_set_type > 0)
498 {
499 MFEM_VERIFY(!include_cut_cell, "include_cut_cell option must be set to"
500 " false for Neumann boundary conditions.");
502 &pmesh, *nbcCoef, *dist_vec,
503 *normalbcCoef, elem_marker,
504 ho_terms, include_cut_cell,
505 ls_cut_marker));
506 bf_neumann_marker.Append(ls_cut_marker);
507 ls_cut_marker += 1;
508 }
509
510 b.Assemble();
511
512 // Set up the bilinear form a(.,.) on the finite element space corresponding
513 // to the Laplacian operator -Delta, by adding the Diffusion domain
514 // integrator and SBM integrator.
515 ParBilinearForm a(&pfespace);
516 ConstantCoefficient one(1.);
518 if (dirichlet_level_set_type > 0)
519 {
521 *dist_vec,
522 elem_marker,
523 bf_dirichlet_marker,
524 include_cut_cell,
525 ho_terms));
527 elem_marker,
528 bf_dirichlet_marker,
529 include_cut_cell,
530 ho_terms), ess_shift_bdr);
531 }
532
533 // Add neumann bilinearform integrator.
534 if (neumann_level_set_type > 0)
535 {
537 *dist_vec,
538 *normalbcCoef,
539 elem_marker,
540 bf_neumann_marker,
541 include_cut_cell,
542 ho_terms));
543 }
544
545 // Assemble the bilinear form and the corresponding linear system,
546 // applying any necessary transformations.
547 a.Assemble();
548
549 // Project the exact solution as an initial condition for Dirichlet boundary.
550 if (!exactCoef)
551 {
552 exactCoef = new FunctionCoefficient(homogeneous);
553 }
554 x.ProjectCoefficient(*exactCoef);
555
556 // Form the linear system and solve it.
557 OperatorPtr A;
558 Vector B, X;
559 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
560
561 Solver *prec = new HypreBoomerAMG;
562 BiCGSTABSolver bicg(MPI_COMM_WORLD);
563 bicg.SetRelTol(1e-12);
564 bicg.SetMaxIter(500);
565 bicg.SetPrintLevel(1);
566 bicg.SetPreconditioner(*prec);
567 bicg.SetOperator(*A);
568 bicg.Mult(B, X);
569
570 // Recover the solution as a finite element grid function.
571 a.RecoverFEMSolution(X, b, x);
572
573 // Save the mesh and the solution.
574 ofstream mesh_ofs("diffusion.mesh");
575 mesh_ofs.precision(8);
576 pmesh.PrintAsOne(mesh_ofs);
577 ofstream sol_ofs("diffusion.gf");
578 sol_ofs.precision(8);
579 x.SaveAsOne(sol_ofs);
580
581 if (visualization)
582 {
583 // Save the solution in ParaView format.
584 ParaViewDataCollection dacol("ParaViewDiffusion", &pmesh);
585 dacol.SetLevelsOfDetail(order);
586 dacol.RegisterField("distance", &distance);
587 dacol.RegisterField("level_set", &level_set_gf);
588 dacol.RegisterField("solution", &x);
589 dacol.SetTime(1.0);
590 dacol.SetCycle(1);
591 dacol.Save();
592
593 // Send the solution by socket to a GLVis server.
594 char vishost[] = "localhost";
595 int visport = 19916, s = 350;
596 socketstream sol_sock;
598 "Solution", s, 0, s, s, "Rj");
599 }
600
601 // Construct an error grid function if the exact solution is known.
602 if (dirichlet_level_set_type == 2 || dirichlet_level_set_type == 3 ||
603 (dirichlet_level_set_type == -1 && neumann_level_set_type == 2))
604 {
605 ParGridFunction error(x);
606 Vector pxyz(dim);
607 pxyz(0) = 0.;
608 for (int i = 0; i < nodes_cnt; i++)
609 {
610 pxyz(0) = vxyz(i);
611 pxyz(1) = vxyz(i+nodes_cnt);
612 real_t exact_val = 0.;
613 if (dirichlet_level_set_type == 2 || neumann_level_set_type == 2)
614 {
615 exact_val = dirichlet_velocity_xy_exponent(pxyz);
616 }
617 else if (dirichlet_level_set_type == 3)
618 {
619 exact_val = dirichlet_velocity_xy_sinusoidal(pxyz);
620 }
621 error(i) = std::fabs(x(i) - exact_val);
622 }
623
624 if (visualization)
625 {
626 char vishost[] = "localhost";
627 int visport = 19916, s = 350;
628 socketstream sol_sock;
629 common::VisualizeField(sol_sock, vishost, visport, error,
630 "Error", 2*s, 0, s, s, "Rj");
631 }
632
633 const real_t global_error = x.ComputeL2Error(*exactCoef);
634 if (myid == 0)
635 {
636 std::cout << "Global L2 error: " << global_error << endl;
637 }
638 }
639
640 const real_t norm = x.ComputeL1Error(one);
641 if (myid == 0) { std::cout << setprecision(10) << norm << std::endl; }
642
643 // Free the used memory.
644 delete prec;
645 delete normalbcCoef;
646 delete nbcCoef;
647 delete dbcCoefCombo;
648 delete dbcCoef;
649 delete exactCoef;
650 delete rhs_f;
651 delete dist_vec;
652 delete neumann_dist_coef;
653 delete dirichlet_dist_coef;
654 delete dirichlet_dist_coef_2;
655
656 return 0;
657}
Combination of level sets: +1 inside the true domain, -1 outside.
Definition sbm_aux.hpp:157
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:144
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
BiCGSTAB method.
Definition solvers.hpp:596
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the BiCGSTAB method.
Definition solvers.cpp:1380
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:609
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:286
Class for domain integration .
Definition lininteg.hpp:109
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
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
Mesh data type.
Definition mesh.hpp:56
bool Conforming() const
Definition mesh.hpp:2228
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition mesh.hpp:1336
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition mesh.cpp:6200
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void SetNodes(const Vector &node_coord)
Updates the vertex/node locations. Invokes NodesUpdated().
Definition mesh.cpp:8985
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
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:285
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
void ProjectDiscCoefficient(VectorCoefficient &coeff) override
Project a discontinuous vector coefficient as a grid function on a continuous finite element space....
real_t ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const override
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:1589
void SetNodalFESpace(FiniteElementSpace *nfes) override
Definition pmesh.cpp:2028
void PrintAsOne(std::ostream &out=mfem::out, const std::string &comments="") const
Definition pmesh.cpp:4968
Helper class for ParaView visualization data.
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
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:683
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:80
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
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)
const int visport
float real_t
Definition config.hpp:43
const char vishost[]
RefCoord s[3]
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