MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
pdiffusion.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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// MFEM ultraweak DPG parallel example for diffusion
13//
14// Compile with: make pdiffusion
15//
16// Sample runs
17// mpirun -np 4 pdiffusion -m ../../data/inline-quad.mesh -o 3 -sref 1 -pref 2 -theta 0.0 -prob 0
18// mpirun -np 4 pdiffusion -m ../../data/inline-hex.mesh -o 2 -sref 0 -pref 1 -theta 0.0 -prob 0 -sc
19// mpirun -np 4 pdiffusion -m ../../data/beam-tet.mesh -o 3 -sref 0 -pref 2 -theta 0.0 -prob 0 -sc
20
21// L-shape runs
22// Note: uniform ref are expected to give sub-optimal rate for the L-shape problem (rate = 2/3)
23// mpirun -np 4 pdiffusion -o 2 -sref 1 -pref 5 -theta 0.0 -prob 1
24
25// L-shape AMR runs
26// mpirun -np 4 pdiffusion -o 1 -sref 1 -pref 10 -theta 0.8 -prob 1
27// mpirun -np 4 pdiffusion -o 2 -sref 1 -pref 8 -theta 0.75 -prob 1 -sc
28// mpirun -np 4 pdiffusion -o 3 -sref 1 -pref 6 -theta 0.75 -prob 1 -sc -do 2
29
30// Description:
31// This example code demonstrates the use of MFEM to define and solve
32// the "ultraweak" (UW) DPG formulation for the Poisson problem in parallel
33
34// - Δ u = f, in Ω
35// u = u₀, on ∂Ω
36//
37// It solves two kinds of problems
38// a) A manufactured solution problem where u_exact = sin(π * (x + y + z)).
39// This example computes and prints out convergence rates for the L2 error.
40// b) The L-shape benchmark problem with AMR. The AMR process is driven by the
41// DPG built-in residual indicator.
42
43// The DPG UW deals with the First Order System
44// ∇ u - σ = 0, in Ω
45// - ∇⋅σ = f, in Ω
46// u = u₀, in ∂Ω
47
48// Ultraweak-DPG is obtained by integration by parts of both equations and the
49// introduction of trace unknowns on the mesh skeleton
50
51// u ∈ L²(Ω), σ ∈ (L²(Ω))ᵈⁱᵐ
52// û ∈ H^1/2, σ̂ ∈ H^-1/2
53// -(u , ∇⋅τ) + < û, τ⋅n> - (σ , τ) = 0, ∀ τ ∈ H(div,Ω)
54// (σ , ∇ v) - < σ̂, v > = (f,v) ∀ v ∈ H¹(Ω)
55// û = u₀ on ∂Ω
56
57// Note:
58// û := u and σ̂ := -σ on the mesh skeleton
59
60// -------------------------------------------------------------
61// | | u | σ | û | σ̂ | RHS |
62// -------------------------------------------------------------
63// | τ | -(u,∇⋅τ) | -(σ,τ) | < û, τ⋅n> | | 0 |
64// | | | | | | |
65// | v | | (σ,∇ v) | | -<σ̂,v> | (f,v) |
66
67// where (τ,v) ∈ H(div,Ω) × H¹(Ω)
68
69// For more information see https://doi.org/10.1007/978-3-319-01818-8_6
70
71#include "mfem.hpp"
72#include "util/pweakform.hpp"
74#include <fstream>
75#include <iostream>
76
77using namespace std;
78using namespace mfem;
79using namespace mfem::common;
80
86
87static const char *enum_str[] =
88{
89 "manufactured",
90 "lshape"
91};
92
94
95real_t exact_u(const Vector & X);
96void exact_gradu(const Vector & X, Vector &gradu);
98void exact_sigma(const Vector & X, Vector & sigma);
99real_t exact_hatu(const Vector & X);
100void exact_hatsigma(const Vector & X, Vector & hatsigma);
101real_t f_exact(const Vector & X);
102
103int main(int argc, char *argv[])
104{
105 // 0. Initialize MPI and HYPRE.
106 Mpi::Init();
107 int myid = Mpi::WorldRank();
108 Hypre::Init();
109
110 // 1. Parse command-line options.
111 const char *mesh_file = "../../data/inline-quad.mesh";
112 int order = 1;
113 int delta_order = 1;
114 int sref = 0; // initial uniform mesh refinements
115 int pref = 0; // parallel mesh refinements for AMR
116 int iprob = 0;
117 bool static_cond = false;
118 real_t theta = 0.7;
119 bool visualization = true;
120 bool paraview = false;
121
122 OptionsParser args(argc, argv);
123 args.AddOption(&mesh_file, "-m", "--mesh",
124 "Mesh file to use.");
125 args.AddOption(&order, "-o", "--order",
126 "Finite element order (polynomial degree).");
127 args.AddOption(&delta_order, "-do", "--delta_order",
128 "Order enrichment for DPG test space.");
129 args.AddOption(&sref, "-sref", "--num-serial-refinements",
130 "Number of initial serial uniform refinements");
131 args.AddOption(&pref, "-pref", "--num-parallel-refinements",
132 "Number of AMR refinements");
133 args.AddOption(&theta, "-theta", "--theta-factor",
134 "Refinement factor (0 indicates uniform refinements) ");
135 args.AddOption(&iprob, "-prob", "--problem", "Problem case"
136 " 0: manufactured, 1: L-shape");
137 args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
138 "--no-static-condensation", "Enable static condensation.");
139 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
140 "--no-visualization",
141 "Enable or disable GLVis visualization.");
142 args.AddOption(&paraview, "-paraview", "--paraview", "-no-paraview",
143 "--no-paraview",
144 "Enable or disable ParaView visualization.");
145 args.Parse();
146 if (!args.Good())
147 {
148 if (myid == 0)
149 {
150 args.PrintUsage(cout);
151 }
152 return 1;
153 }
154
155 if (iprob > 1) { iprob = 1; }
156 prob = (prob_type)iprob;
157
158 if (prob == prob_type::lshape)
159 {
160 mesh_file = "../../data/l-shape.mesh";
161 }
162
163 if (myid == 0)
164 {
165 args.PrintOptions(cout);
166 }
167
168 Mesh mesh(mesh_file, 1, 1);
169 int dim = mesh.Dimension();
170 MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
171
172 if (prob == prob_type::lshape)
173 {
174 /** rotate mesh to be consistent with l-shape benchmark problem
175 See https://doi.org/10.1016/j.amc.2013.05.068 */
176 mesh.EnsureNodes();
177 GridFunction *nodes = mesh.GetNodes();
178 int size = nodes->Size()/2;
179 for (int i = 0; i<size; i++)
180 {
181 real_t x = (*nodes)[2*i];
182 (*nodes)[2*i] = 2*(*nodes)[2*i+1]-1;
183 (*nodes)[2*i+1] = -2*x+1;
184 }
185 }
186
187
188 for (int i = 0; i<sref; i++)
189 {
190 mesh.UniformRefinement();
191 }
192
193 mesh.EnsureNCMesh();
194
195 ParMesh pmesh(MPI_COMM_WORLD, mesh);
196 mesh.Clear();
197
198 // Define spaces
199 enum TrialSpace
200 {
201 u_space = 0,
202 sigma_space = 1,
203 hatu_space = 2,
204 hatsigma_space = 3
205 };
206 enum TestSpace
207 {
208 tau_space = 0,
209 v_space = 1
210 };
211 // L2 space for u
212 FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
213 ParFiniteElementSpace *u_fes = new ParFiniteElementSpace(&pmesh,u_fec);
214
215 // Vector L2 space for σ
216 FiniteElementCollection *sigma_fec = new L2_FECollection(order-1,dim);
217 ParFiniteElementSpace *sigma_fes = new ParFiniteElementSpace(&pmesh,sigma_fec,
218 dim);
219
220 // H^1/2 space for û
221 FiniteElementCollection * hatu_fec = new H1_Trace_FECollection(order,dim);
222 ParFiniteElementSpace *hatu_fes = new ParFiniteElementSpace(&pmesh,hatu_fec);
223
224 // H^-1/2 space for σ̂
225 FiniteElementCollection * hatsigma_fec = new RT_Trace_FECollection(order-1,dim);
226 ParFiniteElementSpace *hatsigma_fes = new ParFiniteElementSpace(&pmesh,
227 hatsigma_fec);
228
229 // testspace fe collections
230 int test_order = order+delta_order;
231 FiniteElementCollection * tau_fec = new RT_FECollection(test_order-1, dim);
232 FiniteElementCollection * v_fec = new H1_FECollection(test_order, dim);
233
236
237 trial_fes.Append(u_fes);
238 trial_fes.Append(sigma_fes);
239 trial_fes.Append(hatu_fes);
240 trial_fes.Append(hatsigma_fes);
241 test_fec.Append(tau_fec);
242 test_fec.Append(v_fec);
243
244 // Required coefficients for the weak formulation
245 ConstantCoefficient one(1.0);
246 ConstantCoefficient negone(-1.0);
247 FunctionCoefficient f(f_exact); // rhs for the manufactured solution problem
248
249 // Required coefficients for the exact solutions
253
254 ParDPGWeakForm * a = new ParDPGWeakForm(trial_fes,test_fec);
255 a->StoreMatrices(true); // this is needed for estimation of residual
256
257 // -(u,∇⋅τ)
258 a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(one),
259 TrialSpace::u_space,TestSpace::tau_space);
260
261 // -(σ,τ)
262 a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(
263 negone)), TrialSpace::sigma_space, TestSpace::tau_space);
264
265 // (σ,∇ v)
266 a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),
267 TrialSpace::sigma_space,TestSpace::v_space);
268
269 // <û,τ⋅n>
270 a->AddTrialIntegrator(new NormalTraceIntegrator,
271 TrialSpace::hatu_space,TestSpace::tau_space);
272
273 // -<σ̂,v> (sign is included in σ̂)
274 a->AddTrialIntegrator(new TraceIntegrator,
275 TrialSpace::hatsigma_space, TestSpace::v_space);
276
277 // test integrators (space-induced norm for H(div) × H1)
278 // (∇⋅τ,∇⋅δτ)
279 a->AddTestIntegrator(new DivDivIntegrator(one),
280 TestSpace::tau_space, TestSpace::tau_space);
281 // (τ,δτ)
282 a->AddTestIntegrator(new VectorFEMassIntegrator(one),
283 TestSpace::tau_space, TestSpace::tau_space);
284 // (∇v,∇δv)
285 a->AddTestIntegrator(new DiffusionIntegrator(one),
286 TestSpace::v_space, TestSpace::v_space);
287 // (v,δv)
288 a->AddTestIntegrator(new MassIntegrator(one),
289 TestSpace::v_space, TestSpace::v_space);
290
291 // RHS
292 if (prob == prob_type::manufactured)
293 {
294 a->AddDomainLFIntegrator(new DomainLFIntegrator(f),TestSpace::v_space);
295 }
296
297 // GridFunction for Dirichlet bdr data
298 ParGridFunction hatu_gf;
299
300 // Visualization streams
301 socketstream u_out;
302 socketstream sigma_out;
303
304 if (myid == 0)
305 {
306 std::cout << "\n Ref |"
307 << " Dofs |"
308 << " L2 Error |"
309 << " Rate |"
310 << " Residual |"
311 << " Rate |"
312 << " PCG it |" << endl;
313 std::cout << std::string(72,'-') << endl;
314 }
315
316 Array<int> elements_to_refine; // for AMR
317 real_t err0 = 0.;
318 int dof0=0.;
319 real_t res0=0.0;
320
321 ParGridFunction u_gf(u_fes);
322 ParGridFunction sigma_gf(sigma_fes);
323 u_gf = 0.0;
324 sigma_gf = 0.0;
325
326 ParaViewDataCollection * paraview_dc = nullptr;
327
328 if (paraview)
329 {
330 paraview_dc = new ParaViewDataCollection(enum_str[prob], &pmesh);
331 paraview_dc->SetPrefixPath("ParaView/Diffusion");
332 paraview_dc->SetLevelsOfDetail(order);
333 paraview_dc->SetCycle(0);
334 paraview_dc->SetDataFormat(VTKFormat::BINARY);
335 paraview_dc->SetHighOrderOutput(true);
336 paraview_dc->SetTime(0.0); // set the time
337 paraview_dc->RegisterField("u",&u_gf);
338 paraview_dc->RegisterField("sigma",&sigma_gf);
339 }
340
341 if (static_cond) { a->EnableStaticCondensation(); }
342 for (int it = 0; it<=pref; it++)
343 {
344 a->Assemble();
345
346 Array<int> ess_tdof_list;
347 Array<int> ess_bdr;
348 if (pmesh.bdr_attributes.Size())
349 {
350 ess_bdr.SetSize(pmesh.bdr_attributes.Max());
351 ess_bdr = 1;
352 hatu_fes->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
353 }
354
355 // shift the ess_tdofs
356 for (int i = 0; i < ess_tdof_list.Size(); i++)
357 {
358 ess_tdof_list[i] += u_fes->GetTrueVSize() + sigma_fes->GetTrueVSize();
359 }
360
361 Array<int> offsets(5);
362 offsets[0] = 0;
363 offsets[1] = u_fes->GetVSize();
364 offsets[2] = sigma_fes->GetVSize();
365 offsets[3] = hatu_fes->GetVSize();
366 offsets[4] = hatsigma_fes->GetVSize();
367 offsets.PartialSum();
368 BlockVector x(offsets);
369 x = 0.0;
370 hatu_gf.MakeRef(hatu_fes,x.GetBlock(2),0);
371 hatu_gf.ProjectBdrCoefficient(uex,ess_bdr);
372
373 Vector X,B;
374 OperatorPtr Ah;
375 a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
376
377 BlockOperator * A = Ah.As<BlockOperator>();
378
380 M.owns_blocks = 1;
381 int skip = 0;
382 if (!static_cond)
383 {
384 HypreBoomerAMG * amg0 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(0,0));
385 HypreBoomerAMG * amg1 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(1,1));
386 amg0->SetPrintLevel(0);
387 amg1->SetPrintLevel(0);
388 M.SetDiagonalBlock(0,amg0);
389 M.SetDiagonalBlock(1,amg1);
390 skip=2;
391 }
393 skip));
394 amg2->SetPrintLevel(0);
395 M.SetDiagonalBlock(skip,amg2);
396 HypreSolver * prec;
397 if (dim == 2)
398 {
399 // AMS preconditioner for 2D H(div) (trace) space
400 prec = new HypreAMS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatsigma_fes);
401 }
402 else
403 {
404 // ADS preconditioner for 3D H(div) (trace) space
405 prec = new HypreADS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatsigma_fes);
406 }
407 M.SetDiagonalBlock(skip+1,prec);
408
409 CGSolver cg(MPI_COMM_WORLD);
410 cg.SetRelTol(1e-12);
411 cg.SetMaxIter(2000);
412 cg.SetPrintLevel(0);
413 cg.SetPreconditioner(M);
414 cg.SetOperator(*A);
415 cg.Mult(B, X);
416
417 a->RecoverFEMSolution(X,x);
418
419 Vector & residuals = a->ComputeResidual(x);
420
421 real_t residual = residuals.Norml2();
422
423 real_t maxresidual = residuals.Max();
424 real_t globalresidual = residual * residual;
425
426 MPI_Allreduce(MPI_IN_PLACE, &maxresidual, 1, MPITypeMap<real_t>::mpi_type,
427 MPI_MAX, MPI_COMM_WORLD);
428 MPI_Allreduce(MPI_IN_PLACE, &globalresidual, 1,
429 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
430
431 globalresidual = sqrt(globalresidual);
432
433 u_gf.MakeRef(u_fes,x.GetBlock(0),0);
434 sigma_gf.MakeRef(sigma_fes,x.GetBlock(1),0);
435
436 int dofs = u_fes->GlobalTrueVSize() + sigma_fes->GlobalTrueVSize()
437 + hatu_fes->GlobalTrueVSize() + hatsigma_fes->GlobalTrueVSize();
438
439 real_t u_err = u_gf.ComputeL2Error(uex);
440 real_t sigma_err = sigma_gf.ComputeL2Error(sigmaex);
441 real_t L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
442 real_t rate_err = (it) ? dim*log(err0/L2Error)/log((real_t)dof0/dofs) : 0.0;
443 real_t rate_res = (it) ? dim*log(res0/globalresidual)/log((
444 real_t)dof0/dofs) : 0.0;
445 err0 = L2Error;
446 res0 = globalresidual;
447 dof0 = dofs;
448
449 if (myid == 0)
450 {
451 std::ios oldState(nullptr);
452 oldState.copyfmt(std::cout);
453 std::cout << std::right << std::setw(5) << it << " | "
454 << std::setw(10) << dof0 << " | "
455 << std::setprecision(3)
456 << std::setw(10) << std::scientific << err0 << " | "
457 << std::setprecision(2)
458 << std::setw(6) << std::fixed << rate_err << " | "
459 << std::setprecision(3)
460 << std::setw(10) << std::scientific << res0 << " | "
461 << std::setprecision(2)
462 << std::setw(6) << std::fixed << rate_res << " | "
463 << std::setw(6) << std::fixed << cg.GetNumIterations() << " | "
464 << std::endl;
465 std::cout.copyfmt(oldState);
466 }
467
468 if (visualization)
469 {
470 const char * keys = (it == 0 && dim == 2) ? "jRcm\n" : nullptr;
471 char vishost[] = "localhost";
472 int visport = 19916;
473
474 VisualizeField(u_out,vishost,visport,u_gf,
475 "Numerical u", 0,0,500,500,keys);
476 VisualizeField(sigma_out,vishost,visport,sigma_gf,
477 "Numerical flux", 500,0,500,500,keys);
478 }
479
480 if (paraview)
481 {
482 paraview_dc->SetCycle(it);
483 paraview_dc->SetTime((real_t)it);
484 paraview_dc->Save();
485 }
486
487 if (it == pref) { break; }
488
489 elements_to_refine.SetSize(0);
490 for (int iel = 0; iel<pmesh.GetNE(); iel++)
491 {
492 if (residuals[iel] >= theta * maxresidual)
493 {
494 elements_to_refine.Append(iel);
495 }
496 }
497
498 pmesh.GeneralRefinement(elements_to_refine);
499
500 for (int i =0; i<trial_fes.Size(); i++)
501 {
502 trial_fes[i]->Update(false);
503 }
504 a->Update();
505 }
506
507 if (paraview)
508 {
509 delete paraview_dc;
510 }
511
512 delete a;
513 delete tau_fec;
514 delete v_fec;
515 delete hatsigma_fes;
516 delete hatsigma_fec;
517 delete hatu_fes;
518 delete hatu_fec;
519 delete sigma_fec;
520 delete sigma_fes;
521 delete u_fec;
522 delete u_fes;
523
524 return 0;
525}
526
528{
529 switch (prob)
530 {
532 {
533 real_t x = X[0];
534 real_t y = X[1];
535 real_t r = sqrt(x*x + y*y);
536 real_t alpha = 2./3.;
537 real_t phi = atan2(y,x);
538 if (phi < 0) { phi += 2*M_PI; }
539 return pow(r,alpha) * sin(alpha * phi);
540 }
541 break;
542 default:
543 {
544 real_t alpha = M_PI * (X.Sum());
545 return sin(alpha);
546 }
547 break;
548 }
549}
550
551void exact_gradu(const Vector & X, Vector & du)
552{
553 du.SetSize(X.Size());
554 switch (prob)
555 {
557 {
558 real_t x = X[0];
559 real_t y = X[1];
560 real_t r = sqrt(x*x + y*y);
561 real_t alpha = 2./3.;
562 real_t phi = atan2(y,x);
563 if (phi < 0) { phi += 2*M_PI; }
564
565 real_t r_x = x/r;
566 real_t r_y = y/r;
567 real_t phi_x = - y / (r*r);
568 real_t phi_y = x / (r*r);
569 real_t beta = alpha * pow(r,alpha - 1.);
570 du[0] = beta*(r_x * sin(alpha*phi) + r * phi_x * cos(alpha*phi));
571 du[1] = beta*(r_y * sin(alpha*phi) + r * phi_y * cos(alpha*phi));
572 }
573 break;
574 default:
575 {
576 real_t alpha = M_PI * (X.Sum());
577 du.SetSize(X.Size());
578 for (int i = 0; i<du.Size(); i++)
579 {
580 du[i] = M_PI * cos(alpha);
581 }
582 }
583 break;
584 }
585}
586
588{
589 switch (prob)
590 {
591 case prob_type::manufactured:
592 {
593 real_t alpha = M_PI * (X.Sum());
594 real_t u = sin(alpha);
595 return - M_PI*M_PI * u * X.Size();
596 }
597 break;
598 default:
599 MFEM_ABORT("Should be unreachable");
600 return 1;
601 break;
602 }
603}
604
605void exact_sigma(const Vector & X, Vector & sigma)
606{
607 // σ = ∇ u
609}
610
612{
613 return exact_u(X);
614}
615
616void exact_hatsigma(const Vector & X, Vector & hatsigma)
617{
618 exact_sigma(X,hatsigma);
619 hatsigma *= -1.;
620}
621
623{
624 MFEM_VERIFY(prob!=prob_type::lshape,
625 "f_exact should not be called for l-shape benchmark problem, i.e., f = 0")
626 return -exact_laplacian_u(X);
627}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void PartialSum()
Fill the entries of the array with the cumulative sum of the entries.
Definition array.cpp:103
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
A class to handle Block diagonal preconditioners in a matrix-free implementation.
A class to handle Block systems in a matrix-free implementation.
Array< int > & RowOffsets()
Return the row offsets for block starts.
Operator & GetBlock(int i, int j)
Return a reference to block i,j.
A class to handle Vectors in a block fashion.
Vector & GetBlock(int i)
Get the i-th vector in the block.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
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)
void SetPrefixPath(const std::string &prefix)
Set the path where the DataCollection will be saved.
for Raviart-Thomas elements
Class for domain integration .
Definition lininteg.hpp:109
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
A general function coefficient.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:322
The Auxiliary-space Divergence Solver in hypre.
Definition hypre.hpp:1922
The Auxiliary-space Maxwell Solver in hypre.
Definition hypre.hpp:1845
The BoomerAMG solver in hypre.
Definition hypre.hpp:1691
void SetPrintLevel(int print_level)
Definition hypre.hpp:1771
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
Abstract class for hypre's solvers and preconditioners.
Definition hypre.hpp:1162
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
int GetNumIterations() const
Returns the number of iterations taken during the last call to Mult()
Definition solvers.hpp:260
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
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10558
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
Definition mesh.cpp:6159
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 GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
void EnsureNCMesh(bool simplices_nonconforming=false)
Definition mesh.cpp:10626
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
Pointer to an Operator of a specified type.
Definition handle.hpp:34
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get().
Definition handle.hpp:104
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, Operator *&A, Vector &X, Vector &B, int copy_interior=0)
Form a constrained linear system using a matrix-free approach.
Definition operator.cpp:114
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 representing the parallel weak formulation. (Convenient for DPG Equations)
Definition pweakform.hpp:26
Abstract parallel finite element space.
Definition pfespace.hpp:29
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
int GetTrueVSize() const override
Return the number of local vector true dofs.
Definition pfespace.hpp:289
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 ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, const Array< int > &attr)
void MakeRef(FiniteElementSpace *f, real_t *v) override
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Class for parallel meshes.
Definition pmesh.hpp:34
Helper class for ParaView visualization data.
void SetHighOrderOutput(bool high_order_output_)
void SetLevelsOfDetail(int levels_of_detail_)
virtual void Save() override
void SetDataFormat(VTKFormat fmt)
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:386
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:445
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:922
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
Vector beta
real_t sigma(const Vector &x)
Definition maxwell.cpp:102
const real_t alpha
Definition ex15.cpp:369
int dim
Definition ex24.cpp:53
prob_type
Definition ex25.cpp:149
@ lshape
Definition ex25.cpp:152
int main()
real_t a
Definition lissajous.cpp:41
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
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
float real_t
Definition config.hpp:43
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
const char vishost[]
void exact_gradu(const Vector &X, Vector &gradu)
real_t exact_hatu(const Vector &X)
real_t exact_laplacian_u(const Vector &X)
real_t exact_u(const Vector &X)
void exact_sigma(const Vector &X, Vector &sigma)
void exact_hatsigma(const Vector &X, Vector &hatsigma)
prob_type prob
prob_type
@ manufactured
@ lshape
real_t f_exact(const Vector &X)
Helper struct to convert a C++ type to an MPI type.