MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
convection-diffusion.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11//
12// MFEM Ultraweak DPG example for convection-diffusion
13//
14// Compile with: make convection-diffusion
15//
16// sample runs
17// convection-diffusion -m ../../data/star.mesh -o 2 -ref 2 -theta 0.0 -eps 1e-1 -beta '2 3'
18// convection-diffusion -m ../../data/beam-hex.mesh -o 2 -ref 2 -theta 0.0 -eps 1e0 -beta '1 0 2'
19// convection-diffusion -m ../../data/inline-tri.mesh -o 3 -ref 2 -theta 0.0 -eps 1e-2 -beta '4 2' -sc
20
21// AMR runs
22// convection-diffusion -o 3 -ref 5 -prob 1 -eps 1e-1 -theta 0.75
23// convection-diffusion -o 2 -ref 9 -prob 1 -eps 1e-2 -theta 0.75
24// convection-diffusion -o 3 -ref 9 -prob 1 -eps 1e-3 -theta 0.75 -sc
25
26// Description:
27// This example code demonstrates the use of MFEM to define and solve
28// the "ultraweak" (UW) DPG formulation for the convection-diffusion problem
29
30// - εΔu + ∇⋅(βu) = f, in Ω
31// u = u_0, on ∂Ω
32
33// It solves the following kinds of problems
34// (a) A manufactured solution where u_exact = sin(π * (x + y + z)).
35// (b) The 2D Erickson-Johnson problem
36
37// The DPG UW deals with the First Order System
38// - ∇⋅σ + ∇⋅(βu) = f, in Ω
39// 1/ε σ - ∇u = 0, in Ω
40// u = u_0, on ∂Ω
41
42// Ultraweak-DPG is obtained by integration by parts of both equations and the
43// introduction of trace unknowns on the mesh skeleton
44//
45// u ∈ L²(Ω), σ ∈ (L²(Ω))ᵈⁱᵐ
46// û ∈ H^1/2, σ̂ ∈ H^-1/2
47// -(βu , ∇v) + (σ , ∇v) + < f̂ , v > = (f,v), ∀ v ∈ H¹(Ω)
48// (u , ∇⋅τ) + 1/ε (σ , τ) + < û , τ⋅n > = 0, ∀ τ ∈ H(div,Ω)
49// û = u_0 on ∂Ω
50
51// Note:
52// f̂ := βu - σ, û := -u on the mesh skeleton
53
54// -------------------------------------------------------------
55// | | u | σ | û | f̂ | RHS |
56// -------------------------------------------------------------
57// | v |-(βu , ∇v) | (σ , ∇v) | | < f̂ ,v > | (f,v) |
58// | | | | | | |
59// | τ | (u ,∇⋅τ) | 1/ε(σ , τ)| <û,τ⋅n> | | 0 |
60
61// where (v,τ) ∈ H¹(Ωₕ) × H(div,Ωₕ)
62
63// For more information see https://doi.org/10.1016/j.camwa.2013.06.010
64
65#include "mfem.hpp"
66#include "util/weakform.hpp"
68#include <fstream>
69#include <iostream>
70
71
72using namespace mfem;
73using namespace mfem::common;
74
76{
78 EJ // see https://doi.org/10.1016/j.camwa.2013.06.010
79};
80
84
85real_t exact_u(const Vector & X);
86void exact_gradu(const Vector & X, Vector & du);
88void exact_sigma(const Vector & X, Vector & sigma);
89real_t exact_hatu(const Vector & X);
90void exact_hatf(const Vector & X, Vector & hatf);
91real_t f_exact(const Vector & X);
93
94int main(int argc, char *argv[])
95{
96 // 1. Parse command-line options.
97 const char *mesh_file = "../../data/inline-quad.mesh";
98 int order = 1;
99 int delta_order = 1;
100 int ref = 1;
101 bool visualization = true;
102 int iprob = 0;
103 real_t theta = 0.0;
104 bool static_cond = false;
105 int visport = 19916;
106 epsilon = 1e0;
107
108 OptionsParser args(argc, argv);
109 args.AddOption(&mesh_file, "-m", "--mesh",
110 "Mesh file to use.");
111 args.AddOption(&order, "-o", "--order",
112 "Finite element order (polynomial degree).");
113 args.AddOption(&delta_order, "-do", "--delta-order",
114 "Order enrichment for DPG test space.");
115 args.AddOption(&epsilon, "-eps", "--epsilon",
116 "Epsilon coefficient");
117 args.AddOption(&ref, "-ref", "--num-refinements",
118 "Number of uniform refinements");
119 args.AddOption(&theta, "-theta", "--theta",
120 "Theta parameter for AMR");
121 args.AddOption(&iprob, "-prob", "--problem", "Problem case"
122 " 0: manufactured, 1: Erickson-Johnson ");
123 args.AddOption(&beta, "-beta", "--beta",
124 "Vector Coefficient beta");
125 args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
126 "--no-static-condensation", "Enable static condensation.");
127 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
128 "--no-visualization",
129 "Enable or disable GLVis visualization.");
130 args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
131 args.Parse();
132 if (!args.Good())
133 {
134 args.PrintUsage(std::cout);
135 return 1;
136 }
137
138 if (iprob > 1) { iprob = 1; }
139 prob = (prob_type)iprob;
140
141 if (prob == prob_type::EJ)
142 {
143 mesh_file = "../../data/inline-quad.mesh";
144 }
145
146 Mesh mesh(mesh_file, 1, 1);
147 int dim = mesh.Dimension();
148 MFEM_VERIFY(dim > 1, "Dimension = 1 is not supported in this example");
149
150 if (beta.Size() == 0)
151 {
153 beta = 0.0;
154 beta[0] = 1.;
155 }
156
157 args.PrintOptions(std::cout);
158
159 // Define spaces
160 enum TrialSpace
161 {
162 u_space = 0,
163 sigma_space = 1,
164 hatu_space = 2,
165 hatf_space = 3
166 };
167 enum TestSpace
168 {
169 v_space = 0,
170 tau_space = 1
171 };
172 // L2 space for u
173 FiniteElementCollection *u_fec = new L2_FECollection(order-1,dim);
174 FiniteElementSpace *u_fes = new FiniteElementSpace(&mesh,u_fec);
175
176 // Vector L2 space for σ
177 FiniteElementCollection *sigma_fec = new L2_FECollection(order-1,dim);
178 FiniteElementSpace *sigma_fes = new FiniteElementSpace(&mesh,sigma_fec, dim);
179
180 // H^1/2 space for û
181 FiniteElementCollection * hatu_fec = new H1_Trace_FECollection(order,dim);
182 FiniteElementSpace *hatu_fes = new FiniteElementSpace(&mesh,hatu_fec);
183
184 // H^-1/2 space for σ̂
185 FiniteElementCollection * hatf_fec = new RT_Trace_FECollection(order-1,dim);
186 FiniteElementSpace *hatf_fes = new FiniteElementSpace(&mesh,hatf_fec);
187
188 // testspace fe collections
189 int test_order = order+delta_order;
190 FiniteElementCollection * v_fec = new H1_FECollection(test_order, dim);
191 FiniteElementCollection * tau_fec = new RT_FECollection(test_order-1, dim);
192
193 // Coefficients
194 ConstantCoefficient one(1.0);
195 ConstantCoefficient negone(-1.0);
198 ConstantCoefficient negeps1(-1./epsilon);
200
203 Vector negbeta = beta; negbeta.Neg();
204 DenseMatrix bbt(beta.Size());
205 MultVVt(beta, bbt);
206 MatrixConstantCoefficient bbtcoeff(bbt);
207 VectorConstantCoefficient negbetacoeff(negbeta);
208
211
212 trial_fes.Append(u_fes);
213 trial_fes.Append(sigma_fes);
214 trial_fes.Append(hatu_fes);
215 trial_fes.Append(hatf_fes);
216 test_fec.Append(v_fec);
217 test_fec.Append(tau_fec);
218
219 FiniteElementCollection *coeff_fec = new L2_FECollection(0,dim);
220 FiniteElementSpace *coeff_fes = new FiniteElementSpace(&mesh,coeff_fec);
221 GridFunction c1_gf, c2_gf;
222 GridFunctionCoefficient c1_coeff(&c1_gf);
223 GridFunctionCoefficient c2_coeff(&c2_gf);
224
225 DPGWeakForm * a = new DPGWeakForm(trial_fes,test_fec);
226 a->StoreMatrices(true); // needed for residual calculation
227
228 //-(βu , ∇v)
229 a->AddTrialIntegrator(new MixedScalarWeakDivergenceIntegrator(betacoeff),
230 TrialSpace::u_space, TestSpace::v_space);
231
232 // (σ,∇ v)
233 a->AddTrialIntegrator(new TransposeIntegrator(new GradientIntegrator(one)),
234 TrialSpace::sigma_space, TestSpace::v_space);
235
236 // (u ,∇⋅τ)
237 a->AddTrialIntegrator(new MixedScalarWeakGradientIntegrator(negone),
238 TrialSpace::u_space, TestSpace::tau_space);
239
240 // 1/ε (σ,τ)
241 a->AddTrialIntegrator(new TransposeIntegrator(new VectorFEMassIntegrator(eps1)),
242 TrialSpace::sigma_space, TestSpace::tau_space);
243
244 // <û,τ⋅n>
245 a->AddTrialIntegrator(new NormalTraceIntegrator,
246 TrialSpace::hatu_space, TestSpace::tau_space);
247
248 // <f̂ ,v>
249 a->AddTrialIntegrator(new TraceIntegrator,
250 TrialSpace::hatf_space, TestSpace::v_space);
251
252 // mesh dependent test norm
253 c1_gf.SetSpace(coeff_fes);
254 c2_gf.SetSpace(coeff_fes);
255 setup_test_norm_coeffs(c1_gf,c2_gf);
256
257 // c1 (v,δv)
258 a->AddTestIntegrator(new MassIntegrator(c1_coeff),
259 TestSpace::v_space, TestSpace::v_space);
260 // ε (∇v,∇δv)
261 a->AddTestIntegrator(new DiffusionIntegrator(eps),
262 TestSpace::v_space, TestSpace::v_space);
263 // (β⋅∇v, β⋅∇δv)
264 a->AddTestIntegrator(new DiffusionIntegrator(bbtcoeff),
265 TestSpace::v_space, TestSpace::v_space);
266 // c2 (τ,δτ)
267 a->AddTestIntegrator(new VectorFEMassIntegrator(c2_coeff),
268 TestSpace::tau_space, TestSpace::tau_space);
269 // (∇⋅τ,∇⋅δτ)
270 a->AddTestIntegrator(new DivDivIntegrator(one),
271 TestSpace::tau_space, TestSpace::tau_space);
272
274 a->AddDomainLFIntegrator(new DomainLFIntegrator(f),TestSpace::v_space);
275
278 Array<int> elements_to_refine;
281
282 GridFunction hatu_gf, hatf_gf;
283
284 socketstream u_out;
285 socketstream sigma_out;
286
287 real_t res0 = 0.;
288 real_t err0 = 0.;
289 int dof0 = 0; // init to suppress gcc warning
290 std::cout << "\n Ref |"
291 << " Dofs |"
292 << " L2 Error |"
293 << " Rate |"
294 << " Residual |"
295 << " Rate |" << std::endl;
296 std::cout << std::string(64,'-') << std::endl;
297
298 if (static_cond) { a->EnableStaticCondensation(); }
299 for (int it = 0; it<=ref; it++)
300 {
301 a->Assemble();
302
303 Array<int> ess_tdof_list_uhat;
304 Array<int> ess_tdof_list_fhat;
305 Array<int> ess_bdr_uhat;
306 Array<int> ess_bdr_fhat;
307 if (mesh.bdr_attributes.Size())
308 {
309 ess_bdr_uhat.SetSize(mesh.bdr_attributes.Max());
310 ess_bdr_fhat.SetSize(mesh.bdr_attributes.Max());
311 ess_bdr_uhat = 1; ess_bdr_fhat = 0;
312 if (prob == prob_type::EJ)
313 {
314 ess_bdr_uhat = 0;
315 ess_bdr_fhat = 1;
316 ess_bdr_uhat[1] = 1;
317 ess_bdr_fhat[1] = 0;
318 }
319 hatu_fes->GetEssentialTrueDofs(ess_bdr_uhat, ess_tdof_list_uhat);
320 hatf_fes->GetEssentialTrueDofs(ess_bdr_fhat, ess_tdof_list_fhat);
321 }
322
323 // shift the ess_tdofs
324 int n = ess_tdof_list_uhat.Size();
325 int m = ess_tdof_list_fhat.Size();
326 Array<int> ess_tdof_list(n+m);
327 for (int j = 0; j < n; j++)
328 {
329 ess_tdof_list[j] = ess_tdof_list_uhat[j]
330 + u_fes->GetTrueVSize()
331 + sigma_fes->GetTrueVSize();
332 }
333 for (int j = 0; j < m; j++)
334 {
335 ess_tdof_list[j+n] = ess_tdof_list_fhat[j]
336 + u_fes->GetTrueVSize()
337 + sigma_fes->GetTrueVSize()
338 + hatu_fes->GetTrueVSize();
339 }
340
341 Array<int> offsets(5);
342 offsets[0] = 0;
343 int dofs = 0;
344 for (int i = 0; i<trial_fes.Size(); i++)
345 {
346 offsets[i+1] = trial_fes[i]->GetVSize();
347 dofs += trial_fes[i]->GetTrueVSize();
348 }
349 offsets.PartialSum();
350
351 BlockVector x(offsets); x = 0.0;
352 hatu_gf.MakeRef(hatu_fes,x.GetBlock(2),0);
353 hatf_gf.MakeRef(hatf_fes,x.GetBlock(3),0);
354 hatu_gf.ProjectBdrCoefficient(hatuex,ess_bdr_uhat);
355 hatf_gf.ProjectBdrCoefficientNormal(hatfex,ess_bdr_fhat);
356
357 OperatorPtr Ah;
358 Vector X,B;
359 a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
360
361 BlockMatrix * A = Ah.As<BlockMatrix>();
362
364 M.owns_blocks = 1;
365 for (int i = 0 ; i < A->NumRowBlocks(); i++)
366 {
367 M.SetDiagonalBlock(i,new DSmoother(A->GetBlock(i,i)));
368 }
369
370 CGSolver cg;
371 cg.SetRelTol(1e-8);
372 cg.SetMaxIter(20000);
373 cg.SetPrintLevel(0);
374 cg.SetPreconditioner(M);
375 cg.SetOperator(*A);
376 cg.Mult(B, X);
377
378 a->RecoverFEMSolution(X,x);
379
380 GridFunction u_gf, sigma_gf;
381 u_gf.MakeRef(u_fes,x.GetBlock(0),0);
382 sigma_gf.MakeRef(sigma_fes,x.GetBlock(1),0);
383
384 real_t u_err = u_gf.ComputeL2Error(uex);
385 real_t sigma_err = sigma_gf.ComputeL2Error(sigmaex);
386 real_t L2Error = sqrt(u_err*u_err + sigma_err*sigma_err);
387
388 Vector & residuals = a->ComputeResidual(x);
389 real_t residual = residuals.Norml2();
390
391 real_t rate_err = (it) ? dim*log(err0/L2Error)/log((real_t)dof0/dofs) : 0.0;
392 real_t rate_res = (it) ? dim*log(res0/residual)/log((real_t)dof0/dofs) : 0.0;
393
394 err0 = L2Error;
395 res0 = residual;
396 dof0 = dofs;
397
398 std::ios oldState(nullptr);
399 oldState.copyfmt(std::cout);
400 std::cout << std::right << std::setw(5) << it << " | "
401 << std::setw(10) << dof0 << " | "
402 << std::setprecision(3)
403 << std::setw(10) << std::scientific << err0 << " | "
404 << std::setprecision(2)
405 << std::setw(6) << std::fixed << rate_err << " | "
406 << std::setprecision(3)
407 << std::setw(10) << std::scientific << res0 << " | "
408 << std::setprecision(2)
409 << std::setw(6) << std::fixed << rate_res << " | "
410 << std::endl;
411 std::cout.copyfmt(oldState);
412
413 if (visualization)
414 {
415 const char * keys = (it == 0 && dim == 2) ? "jRcm\n" : nullptr;
416 char vishost[] = "localhost";
417 VisualizeField(u_out,vishost, visport, u_gf,
418 "Numerical u", 0,0, 500, 500, keys);
419 VisualizeField(sigma_out,vishost, visport, sigma_gf,
420 "Numerical flux", 501,0,500, 500, keys);
421 }
422
423 if (it == ref)
424 {
425 break;
426 }
427
428 elements_to_refine.SetSize(0);
429 real_t max_resid = residuals.Max();
430 for (int iel = 0; iel<mesh.GetNE(); iel++)
431 {
432 if (residuals[iel] > theta * max_resid)
433 {
434 elements_to_refine.Append(iel);
435 }
436 }
437
438 mesh.GeneralRefinement(elements_to_refine,1,1);
439 for (int i =0; i<trial_fes.Size(); i++)
440 {
441 trial_fes[i]->Update(false);
442 }
443 a->Update();
444
445 coeff_fes->Update();
446 c1_gf.Update();
447 c2_gf.Update();
448 setup_test_norm_coeffs(c1_gf,c2_gf);
449 }
450
451 delete coeff_fes;
452 delete coeff_fec;
453 delete a;
454 delete tau_fec;
455 delete v_fec;
456 delete hatf_fes;
457 delete hatf_fec;
458 delete hatu_fes;
459 delete hatu_fec;
460 delete sigma_fes;
461 delete sigma_fec;
462 delete u_fec;
463 delete u_fes;
464
465 return 0;
466}
467
469{
470 real_t x = X[0];
471 real_t y = X[1];
472 real_t z = 0.;
473 if (X.Size() == 3) { z = X[2]; }
474 switch (prob)
475 {
476 case EJ:
477 {
478 real_t alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
479 real_t r1 = (1. + alpha) / (2.*epsilon);
480 real_t r2 = (1. - alpha) / (2.*epsilon);
481 real_t denom = exp(-r2) - exp(-r1);
482
483 real_t g1 = exp(r2*(x-1.));
484 real_t g2 = exp(r1*(x-1.));
485 real_t g = g1-g2;
486
487 return g * cos(M_PI * y)/denom;
488 }
489 break;
490 default:
491 {
492 real_t alpha = M_PI * (x + y + z);
493 return sin(alpha);
494 }
495 break;
496 }
497}
498
499void exact_gradu(const Vector & X, Vector & du)
500{
501 real_t x = X[0];
502 real_t y = X[1];
503 real_t z = 0.;
504 if (X.Size() == 3) { z = X[2]; }
505 du.SetSize(X.Size());
506 du = 0.;
507 switch (prob)
508 {
509 case EJ:
510 {
511 real_t alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
512 real_t r1 = (1. + alpha) / (2.*epsilon);
513 real_t r2 = (1. - alpha) / (2.*epsilon);
514 real_t denom = exp(-r2) - exp(-r1);
515
516 real_t g1 = exp(r2*(x-1.));
517 real_t g1_x = r2*g1;
518 real_t g2 = exp(r1*(x-1.));
519 real_t g2_x = r1*g2;
520 real_t g = g1-g2;
521 real_t g_x = g1_x - g2_x;
522
523 du[0] = g_x * cos(M_PI * y)/denom;
524 du[1] = -M_PI * g * sin(M_PI*y)/denom;
525 }
526 break;
527 default:
528 {
529 real_t alpha = M_PI * (x + y + z);
530 du.SetSize(X.Size());
531 for (int i = 0; i<du.Size(); i++)
532 {
533 du[i] = M_PI * cos(alpha);
534 }
535 }
536 break;
537 }
538}
539
541{
542 real_t x = X[0];
543 real_t y = X[1];
544 real_t z = 0.;
545 if (X.Size() == 3) { z = X[2]; }
546
547 switch (prob)
548 {
549 case EJ:
550 {
551 real_t alpha = sqrt(1. + 4. * epsilon * epsilon * M_PI * M_PI);
552 real_t r1 = (1. + alpha) / (2.*epsilon);
553 real_t r2 = (1. - alpha) / (2.*epsilon);
554 real_t denom = exp(-r2) - exp(-r1);
555
556 real_t g1 = exp(r2*(x-1.));
557 real_t g1_x = r2*g1;
558 real_t g1_xx = r2*g1_x;
559 real_t g2 = exp(r1*(x-1.));
560 real_t g2_x = r1*g2;
561 real_t g2_xx = r1*g2_x;
562 real_t g = g1-g2;
563 real_t g_xx = g1_xx - g2_xx;
564
565 real_t u = g * cos(M_PI * y)/denom;
566 real_t u_xx = g_xx * cos(M_PI * y)/denom;
567 real_t u_yy = -M_PI * M_PI * u;
568 return u_xx + u_yy;
569 }
570 break;
571 default:
572 {
573 real_t alpha = M_PI * (x + y + z);
574 real_t u = sin(alpha);
575 return -M_PI*M_PI * u * X.Size();
576 }
577 break;
578 }
579}
580
581void exact_sigma(const Vector & X, Vector & sigma)
582{
583 // σ = ε ∇ u
585 sigma *= epsilon;
586}
587
589{
590 return -exact_u(X);
591}
592
593void exact_hatf(const Vector & X, Vector & hatf)
594{
597 real_t u = exact_u(X);
598 hatf.SetSize(X.Size());
599 for (int i = 0; i<hatf.Size(); i++)
600 {
601 hatf[i] = beta[i] * u - sigma[i];
602 }
603}
604
606{
607 // f = - εΔu + ∇⋅(βu)
608 Vector du;
609 exact_gradu(X,du);
610 real_t d2u = exact_laplacian_u(X);
611
612 real_t s = 0;
613 for (int i = 0; i<du.Size(); i++)
614 {
615 s += beta[i] * du[i];
616 }
617 return -epsilon * d2u + s;
618}
619
621{
622 Array<int> vdofs;
623 FiniteElementSpace * fes = c1_gf.FESpace();
624 Mesh * mesh = fes->GetMesh();
625 for (int i = 0; i < mesh->GetNE(); i++)
626 {
627 real_t volume = mesh->GetElementVolume(i);
628 real_t c1 = std::min(epsilon/volume, (real_t) 1.);
629 real_t c2 = std::min(1./epsilon, 1./volume);
630 fes->GetElementDofs(i,vdofs);
631 c1_gf.SetSubVector(vdofs,c1);
632 c2_gf.SetSubVector(vdofs,c2);
633 }
634}
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:758
int Size() const
Return the logical size of the array.
Definition array.hpp:147
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:830
A class to handle Block diagonal preconditioners in a matrix-free implementation.
SparseMatrix & GetBlock(int i, int j)
Return a reference to block (i,j). Reference may be invalid if Aij(i,j) == NULL.
Array< int > & RowOffsets()
Return the row offsets for block starts.
int NumRowBlocks() const
Return the number of row blocks.
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:538
void Mult(const Vector &b, Vector &x) const override
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:751
void SetOperator(const Operator &op) override
Set/update the solver for the given operator.
Definition solvers.hpp:551
A coefficient that is constant across space and time.
Class representing the DPG weak formulation. Given the variational formulation a(u,...
Definition weakform.hpp:34
Data type for scaled Jacobi-type smoother of sparse matrix.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
for Raviart-Thomas elements
Class for domain integration .
Definition lininteg.hpp:106
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:851
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
Definition fespace.cpp:3513
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition fespace.cpp:658
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
Definition fespace.cpp:4157
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:679
A general function coefficient.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition gridfunc.cpp:167
virtual void MakeRef(FiniteElementSpace *f, real_t *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition gridfunc.cpp:233
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
Returns ||exsol - u_h||_L2 for scalar or vector H1 or L2 elements.
FiniteElementSpace * FESpace()
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, const Array< int > &bdr_attr)
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition gridfunc.cpp:225
void ProjectBdrCoefficient(Coefficient &coeff, const Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition gridfunc.hpp:481
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:275
Arbitrary order "H^{1/2}-conforming" trace finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:338
void SetRelTol(real_t rtol)
Definition solvers.hpp:229
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:174
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:72
void SetMaxIter(int max_it)
Definition solvers.hpp:231
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:346
A matrix coefficient that is constant in space and time.
Mesh data type.
Definition mesh.hpp:64
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:290
void GeneralRefinement(const Array< Refinement > &refinements, int nonconforming=-1, int nc_limit=0)
Definition mesh.cpp:10883
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
real_t GetElementVolume(int i)
Definition mesh.cpp:121
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 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.
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition fe_coll.hpp:403
Arbitrary order "H^{-1/2}-conforming" face finite elements defined on the interface between mesh elem...
Definition fe_coll.hpp:462
Vector coefficient that is constant in space and time.
A general vector function coefficient.
Vector data type.
Definition vector.hpp:82
void Neg()
(*this) = -(*this)
Definition vector.cpp:375
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:679
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:931
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:1164
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
Vector beta
void setup_test_norm_coeffs(GridFunction &c1_gf, GridFunction &c2_gf)
real_t exact_hatu(const Vector &X)
real_t exact_laplacian_u(const Vector &X)
void exact_gradu(const Vector &X, Vector &du)
real_t epsilon
void exact_hatf(const Vector &X, Vector &hatf)
real_t exact_u(const Vector &X)
void exact_sigma(const Vector &X, Vector &sigma)
prob_type prob
real_t f_exact(const Vector &X)
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
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)
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
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[]