MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
ex34.cpp
Go to the documentation of this file.
1// MFEM Example 34
2//
3// Compile with: make ex34
4//
5// Sample runs: ex34 -o 2
6// ex34 -o 2 -pa -hex
7//
8// Device sample runs:
9// ex34 -o 2 -pa -hex -d cuda
10// ex34 -o 2 -no-pa -d cuda
11//
12// Description: This example code solves a simple magnetostatic problem
13// curl curl A = J where the current density J is computed on a
14// subset of the domain as J = -sigma grad phi. We discretize the
15// vector potential with Nedelec finite elements, the scalar
16// potential with Lagrange finite elements, and the current
17// density with Raviart-Thomas finite elements.
18//
19// The example demonstrates the use of a SubMesh to compute the
20// scalar potential and its associated current density which is
21// then transferred to the original mesh and used as a source
22// function.
23//
24// Note that this example takes certain liberties with the
25// current density which is not necessarily divergence free
26// as it should be. This was done to focus on the use of the
27// SubMesh to transfer information between a full mesh and a
28// sub-domain. A more rigorous implementation might employ an
29// H(div) saddle point solver to obtain a divergence free J on
30// the SubMesh. It would then also need to ensure that the r.h.s.
31// of curl curl A = J does in fact lie in the range of the weak
32// curl operator by performing a divergence cleaning procedure
33// before the solve. After divergence cleaning the delta
34// parameter would probably not be needed.
35//
36// This example is designed to make use of a specific mesh which
37// has a known configuration of elements and boundary attributes.
38// Other meshes could be used but extra care would be required to
39// properly define the SubMesh and the necessary boundaries.
40//
41// We recommend viewing examples 1 and 3 before viewing this
42// example.
43
44#include "mfem.hpp"
45#include <fstream>
46#include <iostream>
47
48using namespace std;
49using namespace mfem;
50
51static bool pa_ = false;
52static bool algebraic_ceed_ = false;
53
55 bool visualization,
56 const Array<int> &phi0_attr,
57 const Array<int> &phi1_attr,
58 const Array<int> &jn_zero_attr,
59 GridFunction &j_cond);
60
61int main(int argc, char *argv[])
62{
63 // 1. Parse command-line options.
64 const char *mesh_file = "../data/fichera-mixed.mesh";
65 Array<int> cond_attr;
66 Array<int> submesh_elems;
67 Array<int> sym_plane_attr;
68 Array<int> phi0_attr;
69 Array<int> phi1_attr;
70 Array<int> jn_zero_attr;
71 int ref_levels = 1;
72 int order = 1;
73 real_t delta_const = 1e-6;
74 bool mixed = true;
75 bool static_cond = false;
76 const char *device_config = "cpu";
77 bool visualization = true;
78
79 OptionsParser args(argc, argv);
81 "Number of times to refine the mesh uniformly.");
83 "Finite element order (polynomial degree).");
85 "Magnetic Conductivity");
87 "--no-static-condensation", "Enable static condensation.");
89 "--hex-mesh", "Mixed mesh of hexahedral mesh.");
91 "--no-partial-assembly", "Enable Partial Assembly.");
93 "Device configuration string, see Device::Configure().");
94#ifdef MFEM_USE_CEED
95 args.AddOption(&algebraic_ceed_, "-a", "--algebraic", "-no-a", "--no-algebraic",
96 "Use algebraic Ceed solver");
97#endif
99 "--no-visualization",
100 "Enable or disable GLVis visualization.");
101
102 args.Parse();
103 if (!args.Good())
104 {
105 args.PrintUsage(cout);
106 return 1;
107 }
108 args.PrintOptions(cout);
109
110 if (!mixed || pa_)
111 {
112 mesh_file = "../data/fichera.mesh";
113 }
114
115 if (submesh_elems.Size() == 0)
116 {
117 if (strcmp(mesh_file, "../data/fichera-mixed.mesh") == 0)
118 {
119 submesh_elems.SetSize(5);
120 submesh_elems[0] = 0;
121 submesh_elems[1] = 2;
122 submesh_elems[2] = 3;
123 submesh_elems[3] = 4;
124 submesh_elems[4] = 9;
125 }
126 else if (strcmp(mesh_file, "../data/fichera.mesh") == 0)
127 {
128 submesh_elems.SetSize(7);
129 submesh_elems[0] = 10;
130 submesh_elems[1] = 14;
131 submesh_elems[2] = 34;
132 submesh_elems[3] = 36;
133 submesh_elems[4] = 37;
134 submesh_elems[5] = 38;
135 submesh_elems[6] = 39;
136 }
137 }
138 if (sym_plane_attr.Size() == 0)
139 {
140 if (strcmp(mesh_file, "../data/fichera-mixed.mesh") == 0 ||
141 strcmp(mesh_file, "../data/fichera.mesh") == 0)
142 {
143 sym_plane_attr.SetSize(8);
144 sym_plane_attr[0] = 9;
145 sym_plane_attr[1] = 10;
146 sym_plane_attr[2] = 11;
147 sym_plane_attr[3] = 12;
148 sym_plane_attr[4] = 13;
149 sym_plane_attr[5] = 14;
150 sym_plane_attr[6] = 15;
151 sym_plane_attr[7] = 16;
152 }
153 }
154 if (phi0_attr.Size() == 0)
155 {
156 if (strcmp(mesh_file, "../data/fichera-mixed.mesh") == 0 ||
157 strcmp(mesh_file, "../data/fichera.mesh") == 0)
158 {
159 phi0_attr.Append(2);
160 }
161 }
162 if (phi1_attr.Size() == 0)
163 {
164 if (strcmp(mesh_file, "../data/fichera-mixed.mesh") == 0 ||
165 strcmp(mesh_file, "../data/fichera.mesh") == 0)
166 {
167 phi1_attr.Append(23);
168 }
169 }
170 if (jn_zero_attr.Size() == 0)
171 {
172 if (strcmp(mesh_file, "../data/fichera-mixed.mesh") == 0 ||
173 strcmp(mesh_file, "../data/fichera.mesh") == 0)
174 {
175 jn_zero_attr.Append(25);
176 }
177 for (int i=0; i<sym_plane_attr.Size(); i++)
178 {
179 jn_zero_attr.Append(sym_plane_attr[i]);
180 }
181 }
182
183 // 2. Enable hardware devices such as GPUs, and programming models such as
184 // CUDA, OCCA, RAJA and OpenMP based on command line options.
185 Device device(device_config);
186 device.Print();
187
188 // 3. Read the (serial) mesh from the given mesh file on all processors. We
189 // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
190 // and volume meshes with the same code.
191 Mesh mesh(mesh_file, 1, 1);
192 int dim = mesh.Dimension();
193
194 if (!mixed || pa_)
195 {
196 mesh.UniformRefinement();
197
198 if (ref_levels > 0)
199 {
200 ref_levels--;
201 }
202 }
203
204 int submesh_attr = -1;
205 if (cond_attr.Size() == 0 && submesh_elems.Size() > 0)
206 {
207 int max_attr = mesh.attributes.Max();
208 submesh_attr = max_attr + 1;
209
210 for (int i=0; i<submesh_elems.Size(); i++)
211 {
212 mesh.SetAttribute(submesh_elems[i], submesh_attr);
213 }
214 mesh.SetAttributes();
215
216 if (cond_attr.Size() == 0)
217 {
218 cond_attr.Append(submesh_attr);
219 }
220 }
221
222 // 4. Refine the serial mesh on all processors to increase the resolution. In
223 // this example we do 'ref_levels' of uniform refinement.
224 {
225 for (int l = 0; l < ref_levels; l++)
226 {
227 mesh.UniformRefinement();
228 }
229 }
230
231 // 5b. Extract a submesh covering a portion of the domain
232 SubMesh mesh_cond(SubMesh::CreateFromDomain(mesh, cond_attr));
233
234 // 6. Define a suitable finite element space on the SubMesh and compute
235 // the current density as an H(div) field.
236 RT_FECollection fec_cond_rt(order - 1, dim);
237 FiniteElementSpace fes_cond_rt(&mesh_cond, &fec_cond_rt);
238 GridFunction j_cond(&fes_cond_rt);
239
240 ComputeCurrentDensityOnSubMesh(order, visualization,
241 phi0_attr, phi1_attr, jn_zero_attr, j_cond);
242
243 // 6a. Save the SubMesh and associated current density in parallel. This
244 // output can be viewed later using GLVis:
245 // "glvis -np <np> -m cond_mesh -g cond_j"
246 {
247 ostringstream mesh_name, cond_name;
248 mesh_name << "cond.mesh";
249 cond_name << "cond_j.gf";
250
251 ofstream mesh_ofs(mesh_name.str().c_str());
252 mesh_ofs.precision(8);
253 mesh_cond.Print(mesh_ofs);
254
255 ofstream cond_ofs(cond_name.str().c_str());
256 cond_ofs.precision(8);
257 j_cond.Save(cond_ofs);
258 }
259
260 // 6b. Send the current density, computed on the SubMesh, to a GLVis server.
261 if (visualization)
262 {
263 char vishost[] = "localhost";
264 int visport = 19916;
265 socketstream port_sock(vishost, visport);
266 port_sock.precision(8);
267 port_sock << "solution\n" << mesh_cond << j_cond
268 << "window_title 'Conductor J'"
269 << "window_geometry 400 0 400 350" << flush;
270 }
271
272 // 7. Define a parallel finite element space on the full mesh. Here we use
273 // the H(curl) finite elements for the vector potential and H(div) for the
274 // current density.
275 ND_FECollection fec_nd(order, dim);
276 RT_FECollection fec_rt(order - 1, dim);
277 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
278 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
279
280 GridFunction j_full(&fespace_rt);
281 j_full = 0.0;
282 mesh_cond.Transfer(j_cond, j_full);
283
284 // 7a. Send the transferred current density to a GLVis server.
285 if (visualization)
286 {
287 char vishost[] = "localhost";
288 int visport = 19916;
289 socketstream sol_sock(vishost, visport);
290 sol_sock.precision(8);
291 sol_sock << "solution\n" << mesh << j_full
292 << "window_title 'J Full'"
293 << "window_geometry 400 430 400 350" << flush;
294 }
295
296 // 8. Determine the list of true (i.e. parallel conforming) essential
297 // boundary dofs. In this example, the boundary conditions are defined by
298 // marking all the boundary attributes except for those on a symmetry
299 // plane as essential (Dirichlet) and converting them to a list of true
300 // dofs.
301 Array<int> ess_tdof_list;
302 Array<int> ess_bdr;
303 if (mesh.bdr_attributes.Size())
304 {
305 ess_bdr.SetSize(mesh.bdr_attributes.Max());
306 ess_bdr = 1;
307 for (int i=0; i<sym_plane_attr.Size(); i++)
308 {
309 ess_bdr[sym_plane_attr[i]-1] = 0;
310 }
311 fespace_nd.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
312 }
313
314 // 9. Set up the parallel linear form b(.) which corresponds to the
315 // right-hand side of the FEM linear system, which in this case is
316 // (J,W_i) where J is given by the function H(div) field transferred
317 // from the SubMesh and W_i are the basis functions in the finite
318 // element fespace.
319 VectorGridFunctionCoefficient jCoef(&j_full);
320 LinearForm b(&fespace_nd);
322 b.Assemble();
323
324 // 10. Define the solution vector x as a parallel finite element grid
325 // function corresponding to fespace. Initialize x to zero.
326 GridFunction x(&fespace_nd);
327 x = 0.0;
328
329 // 11. Set up the parallel bilinear form corresponding to the EM diffusion
330 // operator curl muinv curl + delta I, by adding the curl-curl and the
331 // mass domain integrators. For standard magnetostatics equations choose
332 // delta << 1. Larger values of delta should make the linear system
333 // easier to solve at the expense of resembling a diffusive quasistatic
334 // magnetic field. A reasonable balance must be found whenever the mesh
335 // or problem setup is altered.
336 ConstantCoefficient muinv(1.0);
337 ConstantCoefficient delta(delta_const);
338 BilinearForm a(&fespace_nd);
339 if (pa_) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
342
343 // 12. Assemble the parallel bilinear form and the corresponding linear
344 // system, applying any necessary transformations such as: parallel
345 // assembly, eliminating boundary conditions, applying conforming
346 // constraints for non-conforming AMR, static condensation, etc.
347 if (static_cond) { a.EnableStaticCondensation(); }
348 a.Assemble();
349
350 OperatorPtr A;
351 Vector B, X;
352 a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
353
354 // 13. Solve the system AX=B
355 if (pa_) // Jacobi preconditioning in partial assembly mode
356 {
357 cout << "\nSolving for magnetic vector potential "
358 << "using CG with a Jacobi preconditioner" << endl;
359
360 OperatorJacobiSmoother M(a, ess_tdof_list);
361 PCG(*A, M, B, X, 1, 1000, 1e-12, 0.0);
362 }
363 else
364 {
365#ifndef MFEM_USE_SUITESPARSE
366 cout << "\nSolving for magnetic vector potential "
367 << "using CG with a Gauss-Seidel preconditioner" << endl;
368
369 // 13a. Define a simple symmetric Gauss-Seidel preconditioner and use
370 // it to solve the system Ax=b with PCG.
371 GSSmoother M((SparseMatrix&)(*A));
372 PCG(*A, M, B, X, 1, 500, 1e-12, 0.0);
373#else
374 cout << "\nSolving for magnetic vector potential "
375 << "using UMFPack" << endl;
376
377 // 13a. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
378 // system.
379 UMFPackSolver umf_solver;
380 umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
381 umf_solver.SetOperator(*A);
382 umf_solver.Mult(B, X);
383#endif
384 }
385
386 // 14. Recover the parallel grid function corresponding to X. This is the
387 // local finite element solution on each processor.
388 a.RecoverFEMSolution(X, b, x);
389
390 // 15. Save the refined mesh and the solution in parallel. This output can
391 // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
392 {
393 ostringstream mesh_name, sol_name;
394 mesh_name << "refined.mesh";
395 sol_name << "sol.gf";
396
397 ofstream mesh_ofs(mesh_name.str().c_str());
398 mesh_ofs.precision(8);
399 mesh.Print(mesh_ofs);
400
401 ofstream sol_ofs(sol_name.str().c_str());
402 sol_ofs.precision(8);
403 x.Save(sol_ofs);
404 }
405
406 // 16. Send the solution by socket to a GLVis server.
407 if (visualization)
408 {
409 char vishost[] = "localhost";
410 int visport = 19916;
411 socketstream sol_sock(vishost, visport);
412 sol_sock.precision(8);
413 sol_sock << "solution\n" << mesh << x
414 << "window_title 'Vector Potential'"
415 << "window_geometry 800 0 400 350" << flush;
416 }
417
418 // 17. Compute the magnetic flux as the curl of the solution
419 DiscreteLinearOperator curl(&fespace_nd, &fespace_rt);
421 curl.Assemble();
422 curl.Finalize();
423
424 GridFunction dx(&fespace_rt);
425 curl.Mult(x, dx);
426
427 // 18. Save the curl of the solution in parallel. This output can be viewed
428 // later using GLVis: "glvis -np <np> -m mesh -g dsol".
429 {
430 ostringstream dsol_name;
431 dsol_name << "dsol.gf";
432
433 ofstream dsol_ofs(dsol_name.str().c_str());
434 dsol_ofs.precision(8);
435 dx.Save(dsol_ofs);
436 }
437
438 // 19. Send the curl of the solution by socket to a GLVis server.
439 if (visualization)
440 {
441 char vishost[] = "localhost";
442 int visport = 19916;
443 socketstream sol_sock(vishost, visport);
444 sol_sock.precision(8);
445 sol_sock << "solution\n" << mesh << dx
446 << "window_title 'Magnetic Flux'"
447 << "window_geometry 1200 0 400 350" << flush;
448 }
449
450 // 20. Clean exit
451 return 0;
452}
453
455 bool visualization,
456 const Array<int> &phi0_attr,
457 const Array<int> &phi1_attr,
458 const Array<int> &jn_zero_attr,
459 GridFunction &j_cond)
460{
461 // Extract the finite element space and mesh on which j_cond is defined
462 FiniteElementSpace &fes_cond_rt = *j_cond.FESpace();
463 Mesh &mesh_cond = *fes_cond_rt.GetMesh();
464 int dim = mesh_cond.Dimension();
465
466 // Define a parallel finite element space on the SubMesh. Here we use the H1
467 // finite elements for the electrostatic potential.
468 H1_FECollection fec_h1(order, dim);
469 FiniteElementSpace fes_cond_h1(&mesh_cond, &fec_h1);
470
471 // Define the conductivity coefficient and the boundaries associated with the
472 // fixed potentials phi0 and phi1 which will drive the current.
473 ConstantCoefficient sigmaCoef(1.0);
474 Array<int> ess_bdr_phi(mesh_cond.bdr_attributes.Max());
475 Array<int> ess_bdr_j(mesh_cond.bdr_attributes.Max());
476 Array<int> ess_bdr_tdof_phi;
477 ess_bdr_phi = 0;
478 ess_bdr_j = 0;
479 for (int i=0; i<phi0_attr.Size(); i++)
480 {
481 ess_bdr_phi[phi0_attr[i]-1] = 1;
482 }
483 for (int i=0; i<phi1_attr.Size(); i++)
484 {
485 ess_bdr_phi[phi1_attr[i]-1] = 1;
486 }
487 for (int i=0; i<jn_zero_attr.Size(); i++)
488 {
489 ess_bdr_j[jn_zero_attr[i]-1] = 1;
490 }
491 fes_cond_h1.GetEssentialTrueDofs(ess_bdr_phi, ess_bdr_tdof_phi);
492
493 // Setup the bilinear form corresponding to -Div(sigma Grad phi)
494 BilinearForm a_h1(&fes_cond_h1);
496 a_h1.Assemble();
497
498 // Set the r.h.s. to zero
499 LinearForm b_h1(&fes_cond_h1);
500 b_h1 = 0.0;
501
502 // Setup the boundary conditions on phi
503 ConstantCoefficient one(1.0);
504 ConstantCoefficient zero(0.0);
505 GridFunction phi_h1(&fes_cond_h1);
506 phi_h1 = 0.0;
507
508 Array<int> bdr0(mesh_cond.bdr_attributes.Max()); bdr0 = 0;
509 for (int i=0; i<phi0_attr.Size(); i++)
510 {
511 bdr0[phi0_attr[i]-1] = 1;
512 }
513 phi_h1.ProjectBdrCoefficient(zero, bdr0);
514
515 Array<int> bdr1(mesh_cond.bdr_attributes.Max()); bdr1 = 0;
516 for (int i=0; i<phi1_attr.Size(); i++)
517 {
518 bdr1[phi1_attr[i]-1] = 1;
519 }
520 phi_h1.ProjectBdrCoefficient(one, bdr1);
521
522 {
523 OperatorPtr A;
524 Vector B, X;
525 a_h1.FormLinearSystem(ess_bdr_tdof_phi, phi_h1, b_h1, A, X, B);
526
527 // Solve the linear system
528 if (!pa_)
529 {
530#ifndef MFEM_USE_SUITESPARSE
531 cout << "\nSolving for electric potential using PCG "
532 << "with a Gauss-Seidel preconditioner" << endl;
533
534 // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
535 GSSmoother M((SparseMatrix&)(*A));
536 PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
537#else
538 cout << "\nSolving for electric potential using UMFPack" << endl;
539
540 // If MFEM was compiled with SuiteSparse,
541 // use UMFPACK to solve the system.
542 UMFPackSolver umf_solver;
543 umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
544 umf_solver.SetOperator(*A);
545 umf_solver.Mult(B, X);
546#endif
547 }
548 else
549 {
550 cout << "\nSolving for electric potential using CG" << endl;
551
552 if (UsesTensorBasis(fes_cond_h1))
553 {
554 if (algebraic_ceed_)
555 {
556 ceed::AlgebraicSolver M(a_h1, ess_bdr_tdof_phi);
557 PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
558 }
559 else
560 {
561 OperatorJacobiSmoother M(a_h1, ess_bdr_tdof_phi);
562 PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
563 }
564 }
565 else
566 {
567 CG(*A, B, X, 1, 400, 1e-12, 0.0);
568 }
569 }
570 a_h1.RecoverFEMSolution(X, b_h1, phi_h1);
571 }
572
573 if (visualization)
574 {
575 char vishost[] = "localhost";
576 int visport = 19916;
577 socketstream port_sock(vishost, visport);
578 port_sock.precision(8);
579 port_sock << "solution\n" << mesh_cond << phi_h1
580 << "window_title 'Conductor Potential'"
581 << "window_geometry 0 0 400 350" << flush;
582 }
583
584 // Solve for the current density J = -sigma Grad phi with boundary conditions
585 // J.n = 0 on the walls of the conductor but not on the ports where phi=0 and
586 // phi=1.
587
588 // J will be computed in H(div) so we need an RT mass matrix
589 BilinearForm m_rt(&fes_cond_rt);
591 m_rt.Assemble();
592
593 // Assemble the (sigma Grad phi) operator
594 MixedBilinearForm d_h1(&fes_cond_h1, &fes_cond_rt);
596 d_h1.Assemble();
597
598 // Compute the r.h.s, b_rt = sigma E = -sigma Grad phi
599 LinearForm b_rt(&fes_cond_rt);
600 d_h1.Mult(phi_h1, b_rt);
601 b_rt *= -1.0;
602
603 // Apply the necessary boundary conditions and solve for J in H(div)
604 cout << "\nSolving for current density in H(Div) "
605 << "using diagonally scaled CG" << endl;
606 cout << "Size of linear system: "
607 << fes_cond_rt.GetTrueVSize() << endl;
608
609 Array<int> ess_bdr_tdof_rt;
610 OperatorPtr M;
611 Vector B, X;
612
613 fes_cond_rt.GetEssentialTrueDofs(ess_bdr_j, ess_bdr_tdof_rt);
614
615 j_cond = 0.0;
616 m_rt.FormLinearSystem(ess_bdr_tdof_rt, j_cond, b_rt, M, X, B);
617
618 CGSolver cg;
619 cg.SetRelTol(1e-12);
620 cg.SetMaxIter(2000);
621 cg.SetPrintLevel(1);
622 cg.SetOperator(*M);
623 cg.Mult(B, X);
624 m_rt.RecoverFEMSolution(X, b_rt, j_cond);
625}
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
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:769
A "square matrix" operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
Adds new Domain Integrator. Assumes ownership of bfi.
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
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(....
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.
Integrator for for Nedelec elements.
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition device.hpp:123
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition device.cpp:286
Adds a domain interpolator. Assumes ownership of di.
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
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:588
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
Data type for Gauss-Seidel smoother of sparse matrix.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
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:469
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
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
Vector with associated FE space and LinearFormIntegrators.
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 SetAttribute(int i, int attr)
Set the attribute of element i.
Definition mesh.hpp:1336
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition mesh.hpp:2288
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
virtual void SetAttributes()
Determine the sets of unique attribute values in domain and boundary elements.
Definition mesh.cpp:1604
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY.
void Assemble(int skip_zeros=1)
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
Adds a domain integrator. Assumes ownership of bfi.
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition fe_coll.hpp:465
Pointer to an Operator of a specified type.
Definition handle.hpp:34
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition solvers.hpp:313
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:386
Data type sparse matrix.
Definition sparsemat.hpp:51
Subdomain representation of a topological parent in another Mesh.
Definition submesh.hpp:43
static void Transfer(const GridFunction &src, GridFunction &dst)
Transfer the dofs of a GridFunction.
Definition submesh.cpp:198
static SubMesh CreateFromDomain(const Mesh &parent, Array< int > domain_attributes)
Create a domain SubMesh from its parent.
Definition submesh.cpp:19
Direct sparse solver using UMFPACK.
Definition solvers.hpp:1096
virtual void SetOperator(const Operator &op)
Factorize the given Operator op which must be a SparseMatrix.
Definition solvers.cpp:3102
real_t Control[UMFPACK_CONTROL]
Definition solvers.hpp:1106
virtual void Mult(const Vector &b, Vector &x) const
Direct solution of the linear system using UMFPACK.
Definition solvers.cpp:3197
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition lininteg.hpp:347
Vector coefficient defined by a vector GridFunction.
Vector data type.
Definition vector.hpp:80
Wrapper for AlgebraicMultigrid object.
int dim
Definition ex24.cpp:53
void ComputeCurrentDensityOnSubMesh(int order, bool visualization, const Array< int > &phi0_attr, const Array< int > &phi1_attr, const Array< int > &jn_zero_attr, GridFunction &j_cond)
Definition ex34.cpp:454
int main()
real_t b
Definition lissajous.cpp:42
real_t delta
Definition lissajous.cpp:43
real_t a
Definition lissajous.cpp:41
const int visport
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:913
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, real_t RTOLERANCE, real_t ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
Definition solvers.cpp:898
bool UsesTensorBasis(const FiniteElementSpace &fes)
Return true if the mesh contains only one topology and the elements are tensor elements.
Definition fespace.hpp:1345
float real_t
Definition config.hpp:43
const char vishost[]