MFEM  v4.3.0
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-2021, 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 // mpirun -np 4 diffusion -m ../../data/inline-quad.mesh -rs 3 -o 1 -vis -lst 1
37 // mpirun -np 4 diffusion -m ../../data/inline-hex.mesh -rs 2 -o 2 -vis -lst 1 -ho 1 -alpha 10
38 //
39 // Problem 2: Circular hole of radius 0.2 at the center of the domain.
40 // Solves -nabla^2 u = f with inhomogeneous boundary conditions, and
41 // f is setup such that u = x^p + y^p, where p = 2 by default.
42 // This is a 2D convergence test.
43 // mpirun -np 4 diffusion -rs 2 -o 2 -vis -lst 2
44 //
45 // Problem 3: Domain is y = [0, 1] but mesh is shifted to [-1.e-4, 1].
46 // Solves -nabla^2 u = f with inhomogeneous boundary conditions,
47 // and f is setup such that u = sin(pi*x*y). This is a 2D
48 // convergence test. Second-order can be demonstrated by changing
49 // refinement level (-rs) for the sample run below. Higher-order
50 // convergence can be realized by also increasing the order (-o) of
51 // the finite element space and the number of high-order terms
52 // (-ho) to be included from the Taylor expansion used to enforce
53 // the boundary conditions.
54 // mpirun -np 4 diffusion -rs 2 -o 1 -vis -lst 3
55 //
56 // Problem 4: Complex 2D shape:
57 // Solves -nabla^2 u = 1 with homogeneous boundary conditions.
58 // mpirun -np 4 diffusion -rs 5 -lst 4 -alpha 2
59 
60 #include "mfem.hpp"
61 #include "../common/mfem-common.hpp"
62 #include "sbm_aux.hpp"
63 #include "sbm_solver.hpp"
64 #include "marking.hpp"
65 #include "dist_solver.hpp"
66 
67 using namespace mfem;
68 using namespace std;
69 
70 int main(int argc, char *argv[])
71 {
72 #ifdef HYPRE_USING_CUDA
73  cout << "\nAs of mfem-4.3 and hypre-2.22.0 (July 2021) this miniapp\n"
74  << "is NOT supported with the CUDA version of hypre.\n\n";
75  return 255;
76 #endif
77 
78  // Initialize MPI.
79  int num_procs, myid;
80  MPI_Init(&argc, &argv);
81  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
82  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
83 
84  // Parse command-line options.
85  const char *mesh_file = "../../data/inline-quad.mesh";
86  int order = 2;
87  bool visualization = true;
88  int ser_ref_levels = 0;
89  int level_set_type = 1;
90  int ho_terms = 0;
91  double alpha = 1;
92  bool include_cut_cell = false;
93 
94  OptionsParser args(argc, argv);
95  args.AddOption(&mesh_file, "-m", "--mesh",
96  "Mesh file to use.");
97  args.AddOption(&order, "-o", "--order",
98  "Finite element order (polynomial degree) or -1 for"
99  " isoparametric space.");
100  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
101  "--no-visualization",
102  "Enable or disable GLVis visualization.");
103  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
104  "Number of times to refine the mesh uniformly in serial.");
105  args.AddOption(&level_set_type, "-lst", "--level-set-type",
106  "level-set-type:");
107  args.AddOption(&ho_terms, "-ho", "--high-order",
108  "Additional high-order terms to include");
109  args.AddOption(&alpha, "-alpha", "--alpha",
110  "Nitsche penalty parameter (~1 for 2D, ~10 for 3D).");
111  args.AddOption(&include_cut_cell, "-cut", "--cut", "-no-cut-cell",
112  "--no-cut-cell",
113  "Include or not include elements cut by true boundary.");
114 
115  args.Parse();
116  if (!args.Good())
117  {
118  if (myid == 0)
119  {
120  args.PrintUsage(cout);
121  }
122  MPI_Finalize();
123  return 1;
124  }
125  if (myid == 0) { args.PrintOptions(cout); }
126 
127  // Enable hardware devices such as GPUs, and programming models such as CUDA,
128  // OCCA, RAJA and OpenMP based on command line options.
129  Device device("cpu");
130  if (myid == 0) { device.Print(); }
131 
132  // Refine the mesh.
133  Mesh mesh(mesh_file, 1, 1);
134  int dim = mesh.Dimension();
135  for (int lev = 0; lev < ser_ref_levels; lev++) { mesh.UniformRefinement(); }
136 
137  // MPI distribution.
138  ParMesh pmesh(MPI_COMM_WORLD, mesh);
139  mesh.Clear();
140 
141  // Define a finite element space on the mesh. Here we use continuous Lagrange
142  // finite elements of the specified order. If order < 1, we fix it to 1.
143  if (order < 1) { order = 1; }
144  H1_FECollection fec(order, dim);
145  ParFiniteElementSpace pfespace(&pmesh, &fec);
146 
147  Vector vxyz;
148 
149  // Set the nodal grid function for the mesh, and modify the nodal positions
150  // for level_set_type = 3 such that some of the mesh elements are intersected
151  // by the true boundary (y = 0).
152  ParFiniteElementSpace pfespace_mesh(&pmesh, &fec, dim);
153  pmesh.SetNodalFESpace(&pfespace_mesh);
154  ParGridFunction x_mesh(&pfespace_mesh);
155  pmesh.SetNodalGridFunction(&x_mesh);
156  vxyz = *pmesh.GetNodes();
157  int nodes_cnt = vxyz.Size()/dim;
158  if (level_set_type == 3)
159  {
160  for (int i = 0; i < nodes_cnt; i++)
161  {
162  // Shift the mesh from y = [0, 1] to [-1.e-4, 1].
163  vxyz(i+nodes_cnt) = (1.+1.e-4)*vxyz(i+nodes_cnt)-1.e-4;
164  }
165  }
166  pmesh.SetNodes(vxyz);
167  pfespace.ExchangeFaceNbrData();
168  const int gtsize = pfespace.GlobalTrueVSize();
169  if (myid == 0)
170  {
171  cout << "Number of finite element unknowns: " << gtsize << endl;
172  }
173 
174  // Define the solution vector x as a finite element grid function
175  // corresponding to pfespace.
176  ParGridFunction x(&pfespace);
177  // ParGridFunction for level_set_value.
178  ParGridFunction level_set_val(&pfespace);
179 
180  // Determine if each element in the ParMesh is inside the actual domain,
181  // partially cut by its boundary, or completely outside the domain.
182  Dist_Level_Set_Coefficient dist_fun_level_coef(level_set_type);
183  level_set_val.ProjectCoefficient(dist_fun_level_coef);
184  // Exchange information for ghost elements i.e. elements that share a face
185  // with element on the current processor, but belong to another processor.
186  level_set_val.ExchangeFaceNbrData();
187  // Setup the class to mark all elements based on whether they are located
188  // inside or outside the true domain, or intersected by the true boundary.
189  ShiftedFaceMarker marker(pmesh, level_set_val, pfespace, include_cut_cell);
190  Array<int> elem_marker;
191  marker.MarkElements(elem_marker);
192 
193  // Visualize the element markers.
194  if (visualization)
195  {
196  L2_FECollection fecl2 = L2_FECollection(0, dim);
197  ParFiniteElementSpace pfesl2(&pmesh, &fecl2);
198  ParGridFunction elem_marker_gf(&pfesl2);
199  for (int i = 0; i < elem_marker_gf.Size(); i++)
200  {
201  elem_marker_gf(i) = (double)elem_marker[i];
202  }
203  char vishost[] = "localhost";
204  int visport = 19916, s = 350;
205  socketstream sol_sock;
206  common::VisualizeField(sol_sock, vishost, visport, elem_marker_gf,
207  "Element Flags", 0, 0, s, s, "Rjmpc");
208  }
209 
210  // Get a list of dofs associated with shifted boundary (SB) faces.
211  Array<int> sb_dofs; // Array of dofs on SB faces
212  marker.ListShiftedFaceDofs(elem_marker, sb_dofs);
213 
214  // Visualize the shifted boundary face dofs.
215  if (visualization)
216  {
217  ParGridFunction face_dofs(&pfespace);
218  face_dofs = 0.0;
219  for (int i = 0; i < sb_dofs.Size(); i++)
220  {
221  face_dofs(sb_dofs[i]) = 1.0;
222  }
223  char vishost[] = "localhost";
224  int visport = 19916, s = 350;
225  socketstream sol_sock;
226  common::VisualizeField(sol_sock, vishost, visport, face_dofs,
227  "Shifted Face Dofs", 0, s, s, s, "Rjmp");
228  }
229 
230  // Make a list of inactive tdofs that will be eliminated from the system.
231  // The inactive tdofs are the dofs for the elements located outside the true
232  // domain (and optionally, for the elements cut by the true boundary, if
233  // include_cut_cell = false) minus the dofs that are located on the surrogate
234  // boundary.
235  Array<int> ess_tdof_list;
236  Array<int> ess_shift_bdr;
237  marker.ListEssentialTDofs(elem_marker, sb_dofs, ess_tdof_list,
238  ess_shift_bdr);
239 
240  // Compute distance vector to the actual boundary.
241  ParFiniteElementSpace distance_vec_space(&pmesh, &fec, dim);
242  ParGridFunction distance(&distance_vec_space);
243  VectorCoefficient *dist_vec = NULL;
244  // Compute the distance field using the HeatDistanceSolver for
245  // level_set_type == 4 or analytically for all other level set types.
246  if (level_set_type == 4)
247  {
248  // Discrete distance vector.
249  double dx = AvgElementSize(pmesh);
250  ParGridFunction filt_gf(&pfespace);
251  PDEFilter *filter = new PDEFilter(pmesh, 2.0 * dx);
252  filter->Filter(dist_fun_level_coef, filt_gf);
253  delete filter;
254  GridFunctionCoefficient ls_filt_coeff(&filt_gf);
255 
256  if (visualization)
257  {
258  char vishost[] = "localhost";
259  int visport = 19916, s = 350;
260  socketstream sol_sock;
261  common::VisualizeField(sol_sock, vishost, visport, filt_gf,
262  "Input Level Set", 0, 2*s, s, s, "Rjmm");
263  }
264 
265  HeatDistanceSolver dist_func(2.0 * dx* dx);
266  dist_func.print_level = 1;
267  dist_func.smooth_steps = 1;
268  dist_func.ComputeVectorDistance(ls_filt_coeff, distance);
269  dist_vec = new VectorGridFunctionCoefficient(&distance);
270  }
271  else
272  {
273  // Analytic distance vector.
274  dist_vec = new Dist_Vector_Coefficient(dim, level_set_type);
275  distance.ProjectDiscCoefficient(*dist_vec);
276  }
277 
278  // Visualize the distance vector.
279  if (visualization)
280  {
281  char vishost[] = "localhost";
282  int visport = 19916, s = 350;
283  socketstream sol_sock;
284  common::VisualizeField(sol_sock, vishost, visport, distance,
285  "Distance Vector", s, s, s, s, "Rjmmpcvv", 1);
286  }
287 
288  // Set up a list to indicate element attributes to be included in assembly,
289  // so that inactive elements are excluded.
290  const int max_elem_attr = pmesh.attributes.Max();
291  Array<int> ess_elem(max_elem_attr);
292  ess_elem = 1;
293  bool inactive_elements = false;
294  for (int i = 0; i < pmesh.GetNE(); i++)
295  {
296  if (!include_cut_cell &&
297  (elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE ||
298  elem_marker[i] == ShiftedFaceMarker::SBElementType::CUT))
299  {
300  pmesh.SetAttribute(i, max_elem_attr+1);
301  inactive_elements = true;
302  }
303  if (include_cut_cell &&
304  elem_marker[i] == ShiftedFaceMarker::SBElementType::OUTSIDE)
305  {
306  pmesh.SetAttribute(i, max_elem_attr+1);
307  inactive_elements = true;
308  }
309  }
310  bool inactive_elements_global;
311  MPI_Allreduce(&inactive_elements, &inactive_elements_global, 1, MPI_C_BOOL,
312  MPI_LOR, MPI_COMM_WORLD);
313  if (inactive_elements_global) { ess_elem.Append(0); }
314  pmesh.SetAttributes();
315 
316  // Set up the linear form b(.) which corresponds to the right-hand side of
317  // the FEM linear system.
318  ParLinearForm b(&pfespace);
319  FunctionCoefficient *rhs_f = NULL;
320  if (level_set_type == 1 || level_set_type == 4)
321  {
322  rhs_f = new FunctionCoefficient(rhs_fun_circle);
323  }
324  else if (level_set_type == 2)
325  {
327  }
328  else if (level_set_type == 3)
329  {
331  }
332  else { MFEM_ABORT("RHS function not set for level set type.\n"); }
333  b.AddDomainIntegrator(new DomainLFIntegrator(*rhs_f), ess_elem);
334 
335  // Dirichlet BC that must be imposed on the true boundary.
336  ShiftedFunctionCoefficient *dbcCoef = NULL;
337  if (level_set_type == 1 || level_set_type == 4)
338  {
340  }
341  else if (level_set_type == 2)
342  {
344  }
345  else if (level_set_type == 3)
346  {
348  }
349  else
350  {
351  MFEM_ABORT("Dirichlet velocity function not set for level set type.\n");
352  }
353  // Add integrators corresponding to the shifted boundary method (SBM).
355  alpha, *dist_vec,
356  elem_marker,
357  include_cut_cell,
358  ho_terms));
359  b.AddBdrFaceIntegrator(new SBM2DirichletLFIntegrator(&pmesh, *dbcCoef,
360  alpha, *dist_vec,
361  elem_marker,
362  include_cut_cell,
363  ho_terms), ess_shift_bdr);
364  b.Assemble();
365 
366  // Set up the bilinear form a(.,.) on the finite element space corresponding
367  // to the Laplacian operator -Delta, by adding the Diffusion domain
368  // integrator and SBM integrator.
369  ParBilinearForm a(&pfespace);
370  ConstantCoefficient one(1.);
371  a.AddDomainIntegrator(new DiffusionIntegrator(one), ess_elem);
373  *dist_vec,
374  elem_marker,
375  include_cut_cell,
376  ho_terms));
377  a.AddBdrFaceIntegrator(new SBM2DirichletIntegrator(&pmesh, alpha, *dist_vec,
378  elem_marker,
379  include_cut_cell,
380  ho_terms), ess_shift_bdr);
381 
382  // Assemble the bilinear form and the corresponding linear system,
383  // applying any necessary transformations.
384  a.Assemble();
385 
386  // Project the exact solution as an initial condition for Dirichlet boundary.
387  x.ProjectCoefficient(*dbcCoef);
388 
389  // Form the linear system and solve it.
390  OperatorPtr A;
391  Vector B, X;
392  a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
393 
394  Solver *prec = new HypreBoomerAMG;
395  BiCGSTABSolver *bicg = new BiCGSTABSolver(MPI_COMM_WORLD);
396  bicg->SetRelTol(1e-12);
397  bicg->SetMaxIter(2000);
398  bicg->SetPrintLevel(1);
399  bicg->SetPreconditioner(*prec);
400  bicg->SetOperator(*A);
401  bicg->Mult(B, X);
402 
403  // Recover the solution as a finite element grid function.
404  a.RecoverFEMSolution(X, b, x);
405 
406  // Save the mesh and the solution.
407  ofstream mesh_ofs("diffusion.mesh");
408  mesh_ofs.precision(8);
409  pmesh.PrintAsOne(mesh_ofs);
410  ofstream sol_ofs("diffusion.gf");
411  sol_ofs.precision(8);
412  x.SaveAsOne(sol_ofs);
413 
414  // Save the solution in ParaView format
415  if (visualization)
416  {
417  ParaViewDataCollection dacol("ParaViewDiffusion", &pmesh);
418  dacol.SetLevelsOfDetail(order);
419  dacol.RegisterField("distance", &distance);
420  dacol.RegisterField("solution", &x);
421  dacol.SetTime(1.0);
422  dacol.SetCycle(1);
423  dacol.Save();
424  }
425 
426 
427  // Send the solution by socket to a GLVis server.
428  if (visualization)
429  {
430  char vishost[] = "localhost";
431  int visport = 19916, s = 350;
432  socketstream sol_sock;
433  common::VisualizeField(sol_sock, vishost, visport, x,
434  "Solution", s, 0, s, s, "Rj");
435  }
436 
437  // Construct an error grid function if the exact solution is known.
438  if (level_set_type == 2 || level_set_type == 3)
439  {
440  ParGridFunction err(x);
441  Vector pxyz(dim);
442  pxyz(0) = 0.;
443  for (int i = 0; i < nodes_cnt; i++)
444  {
445  pxyz(0) = vxyz(i);
446  pxyz(1) = vxyz(i+nodes_cnt);
447  double exact_val = 0.;
448  if (level_set_type == 2)
449  {
450  exact_val = dirichlet_velocity_xy_exponent(pxyz);
451  }
452  else if (level_set_type == 3)
453  {
454  exact_val = dirichlet_velocity_xy_sinusoidal(pxyz);
455  }
456  err(i) = std::fabs(x(i) - exact_val);
457  }
458 
459  if (visualization)
460  {
461  char vishost[] = "localhost";
462  int visport = 19916, s = 350;
463  socketstream sol_sock;
464  common::VisualizeField(sol_sock, vishost, visport, err,
465  "Error", 2*s, 0, s, s, "Rj");
466  }
467 
468  const double global_error = x.ComputeL2Error(*dbcCoef);
469  if (myid == 0)
470  {
471  std::cout << "Global L2 error: " << global_error << endl;
472  }
473  }
474 
475  // Free the used memory.
476  delete prec;
477  delete bicg;
478  delete dbcCoef;
479  delete rhs_f;
480  delete dist_vec;
481 
482  MPI_Finalize();
483 
484  return 0;
485 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
double rhs_fun_circle(const Vector &x)
f for the Poisson problem (-nabla^2 u = f).
Definition: sbm_aux.hpp:149
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:78
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(...
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:275
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
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:514
void MarkElements(Array< int > &elem_marker) const
Mark all the elements in the mesh using the SBElementType.
Definition: marking.cpp:17
void AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Interior Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:97
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:493
virtual void SetAttributes()
Definition: pmesh.cpp:1574
double dirichlet_velocity_xy_exponent(const Vector &x)
Definition: sbm_aux.hpp:136
double dirichlet_velocity_circle(const Vector &x)
Boundary conditions.
Definition: sbm_aux.hpp:131
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
Class for parallel linear form.
Definition: plinearform.hpp:26
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection.
void SetPrintLevel(int print_lvl)
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:150
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:746
constexpr char vishost[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
void ListShiftedFaceDofs(const Array< int > &elem_marker, Array< int > &sface_dof_list) const
Definition: marking.cpp:84
constexpr int visport
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:180
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Distance vector to the zero level-set.
Definition: sbm_aux.hpp:96
Level set coefficient - +1 inside the domain, -1 outside, 0 at the boundary.
Definition: sbm_aux.hpp:76
void AddBdrFaceIntegrator(LinearFormIntegrator *lfi)
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:83
void SetNodalFESpace(FiniteElementSpace *nfes)
Definition: mesh.cpp:4843
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
void SetTime(double t)
Set physical time (for time-dependent simulations)
double rhs_fun_xy_sinusoidal(const Vector &x)
Definition: sbm_aux.hpp:169
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
void SetRelTol(double rtol)
Definition: solvers.hpp:98
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:1280
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:154
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:867
double a
Definition: lissajous.cpp:41
double dirichlet_velocity_xy_sinusoidal(const Vector &x)
Definition: sbm_aux.hpp:142
virtual double ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL) const
Definition: pgridfunc.hpp:273
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:409
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
Class for parallel bilinear form.
void SetLevelsOfDetail(int levels_of_detail_)
void Clear()
Clear the contents of the Mesh.
Definition: mesh.hpp:828
const double alpha
Definition: ex15.cpp:369
A general function coefficient.
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1200
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:94
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
Vector coefficient defined by a vector GridFunction.
void PrintAsOne(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4527
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
virtual void Save() override
RefCoord s[3]
BiCGSTAB method.
Definition: solvers.hpp:395
void SetNodes(const Vector &node_coord)
Definition: mesh.cpp:7355
Base class for solvers.
Definition: operator.hpp:648
Class for parallel grid function.
Definition: pgridfunc.hpp:32
void SetNodalGridFunction(GridFunction *nodes, bool make_owner=false)
Definition: mesh.cpp:4871
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
int main()
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:202
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:285
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