MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
diffusion.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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"
75 #include "../common/mfem-common.hpp"
76 #include "sbm_aux.hpp"
77 #include "sbm_solver.hpp"
78 #include "marking.hpp"
79 #include "dist_solver.hpp"
80 
81 using namespace mfem;
82 using namespace std;
83 
84 int 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 242;
90 #endif
91 
92  // Initialize MPI and HYPRE.
93  Mpi::Init(argc, argv);
94  int myid = Mpi::WorldRank();
95  Hypre::Init();
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  double alpha = 1;
107  bool include_cut_cell = false;
108 
109  OptionsParser args(argc, argv);
110  args.AddOption(&mesh_file, "-m", "--mesh",
111  "Mesh file to use.");
112  args.AddOption(&order, "-o", "--order",
113  "Finite element order (polynomial degree) or -1 for"
114  " isoparametric space.");
115  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
116  "--no-visualization",
117  "Enable or disable GLVis visualization.");
118  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
119  "Number of times to refine the mesh uniformly in serial.");
120  args.AddOption(&dirichlet_level_set_type, "-lst", "--level-set-type",
121  "level-set-type.");
122  args.AddOption(&neumann_level_set_type, "-nlst", "--neumann-level-set-type",
123  "neumann-level-set-type.");
124  args.AddOption(&dirichlet_combo, "-dc", "--dcombo",
125  "no-dc", "--no-dcombo",
126  "Combination of two Dirichlet level sets.");
127  args.AddOption(&ho_terms, "-ho", "--high-order",
128  "Additional high-order terms to include");
129  args.AddOption(&alpha, "-alpha", "--alpha",
130  "Nitsche penalty parameter (~1 for 2D, ~10 for 3D).");
131  args.AddOption(&include_cut_cell, "-cut", "--cut", "-no-cut-cell",
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 double 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);
235  combo_dist_coef.Add_Level_Set_Coefficient(*dirichlet_dist_coef);
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);
248  combo_dist_coef.Add_Level_Set_Coefficient(*dirichlet_dist_coef_2);
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);
258  combo_dist_coef.Add_Level_Set_Coefficient(*neumann_dist_coef);
259  }
260 
261  // Visualize the element markers.
262  if (visualization)
263  {
264  L2_FECollection fecl2 = L2_FECollection(0, dim);
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) = (double)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  double 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  {
392  rhs_f = new FunctionCoefficient(rhs_fun_circle);
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"); }
403  b.AddDomainIntegrator(new DomainLFIntegrator(*rhs_f), ess_elem);
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  {
437  normalbcCoef = new ShiftedVectorFunctionCoefficient(dim, normal_vector_1);
438  }
439  else if (neumann_level_set_type == 2)
440  {
442  normalbcCoef = new ShiftedVectorFunctionCoefficient(dim, normal_vector_1);
444  }
445  else if (neumann_level_set_type == 7)
446  {
448  normalbcCoef = new ShiftedVectorFunctionCoefficient(dim, normal_vector_2);
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));
472  b.AddBdrFaceIntegrator(new SBM2DirichletLFIntegrator(&pmesh, *dbcCoef,
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  {
485  b.AddInteriorFaceIntegrator(new SBM2DirichletLFIntegrator(&pmesh, *dbcCoefCombo,
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.);
517  a.AddDomainIntegrator(new DiffusionIntegrator(one), ess_elem);
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));
526  a.AddBdrFaceIntegrator(new SBM2DirichletIntegrator(&pmesh, alpha, *dist_vec,
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;
597  common::VisualizeField(sol_sock, vishost, visport, x,
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  double 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 double global_error = x.ComputeL2Error(*exactCoef);
634  if (myid == 0)
635  {
636  std::cout << "Global L2 error: " << global_error << endl;
637  }
638  }
639 
640  const double 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 }
const char vishost[]
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(), hypre will be finalized automatically at program exit.
Definition: hypre.hpp:80
double rhs_fun_circle(const Vector &x)
f for the Poisson problem (-nabla^2 u = f).
Definition: sbm_aux.hpp:268
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Base class for vector Coefficients that optionally depend on time and space.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
bool Conforming() const
Definition: mesh.hpp:1650
Helper class for ParaView visualization data.
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
Project a discontinuous vector coefficient as a grid function on a continuous finite element space...
Definition: pgridfunc.cpp:546
void AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Interior Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:99
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
double dirichlet_velocity_xy_exponent(const Vector &x)
Definition: sbm_aux.hpp:223
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
Class for parallel linear form.
Definition: plinearform.hpp:26
double homogeneous(const Vector &x)
Boundary conditions - Dirichlet.
Definition: sbm_aux.hpp:218
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition: solvers.cpp:71
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:751
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9816
void ListShiftedFaceDofs(const Array< int > &elem_marker, Array< int > &sface_dof_list) const
Definition: marking.cpp:106
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
void SetMaxIter(int max_it)
Definition: solvers.hpp:200
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Definition: pgridfunc.hpp:281
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Distance vector to the zero level-set.
Definition: sbm_aux.hpp:183
Level set coefficient: +1 inside the true domain, -1 outside.
Definition: sbm_aux.hpp:137
static void Init()
Singleton creation with Mpi::Init();.
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:85
const int visport
int Dimension() const
Definition: mesh.hpp:1006
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
void SetTime(double t)
Set physical time (for time-dependent simulations)
double rhs_fun_xy_sinusoidal(const Vector &x)
Definition: sbm_aux.hpp:288
virtual double ComputeL1Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:266
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
IterativeSolver::PrintLevel print_level
Definition: dist_solver.hpp:29
void SetRelTol(double rtol)
Definition: solvers.hpp:198
Combination of level sets: +1 inside the true domain, -1 outside.
Definition: sbm_aux.hpp:156
void Add_Level_Set_Coefficient(Dist_Level_Set_Coefficient &dls_)
Definition: sbm_aux.hpp:164
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1373
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
virtual void ComputeVectorDistance(Coefficient &zero_level_set, ParGridFunction &distance)
double rhs_fun_xy_exponent(const Vector &x)
Definition: sbm_aux.hpp:273
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:899
double a
Definition: lissajous.cpp:41
double dirichlet_velocity_xy_sinusoidal(const Vector &x)
Definition: sbm_aux.hpp:229
void normal_vector_1(const Vector &x, Vector &p)
Definition: sbm_aux.hpp:236
int dim
Definition: ex24.cpp:53
void AddInteriorFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new interior Face Integrator. Assumes ownership of bfi.
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:558
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Class for parallel bilinear form.
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
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
void SetLevelsOfDetail(int levels_of_detail_)
double traction_xy_exponent(const Vector &x)
Neumann condition for exponent based solution.
Definition: sbm_aux.hpp:256
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:899
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
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)
Definition: fem_extras.cpp:92
Vector coefficient defined by a vector GridFunction.
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:173
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:220
virtual void Save() override
RefCoord s[3]
BiCGSTAB method.
Definition: solvers.hpp:544
Base class for solvers.
Definition: operator.hpp:655
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
Class for parallel meshes.
Definition: pmesh.hpp:32
void MarkElements(const ParGridFunction &ls_func, Array< int > &elem_marker)
Definition: marking.cpp:17
int main()
void Filter(ParGridFunction &func, ParGridFunction &ffield)
void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi)
Adds new boundary Face Integrator. Assumes ownership of bfi.
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:288
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