MFEM v2.0
ex3.cpp
Go to the documentation of this file.
```00001 //                                MFEM Example 3
00002 //
00003 // Compile with: make ex3
00004 //
00005 // Sample runs:  ex3 ../data/beam-tet.mesh
00006 //               ex3 ../data/beam-hex.mesh
00007 //               ex3 ../data/escher.mesh
00008 //               ex3 ../data/fichera.mesh
00009 //               ex3 ../data/fichera-q2.vtk
00010 //               ex3 ../data/fichera-q3.mesh
00011 //               ex3 ../data/beam-hex-nurbs.mesh
00012 //
00013 // Description:  This example code solves a simple 3D electromagnetic diffusion
00014 //               problem corresponding to the second order definite Maxwell
00015 //               equation curl curl E + E = f with boundary condition
00016 //               E x n = <given tangential field>. Here, we use a given exact
00017 //               solution E and compute the corresponding r.h.s. f.
00018 //               We discretize with Nedelec finite elements.
00019 //
00020 //               The example demonstrates the use of H(curl) finite element
00021 //               spaces with the curl-curl and the (vector finite element) mass
00022 //               bilinear form, as well as the computation of discretization
00023 //               error when the exact solution is known.
00024 //
00025 //               We recommend viewing examples 1-2 before viewing this example.
00026
00027 #include <fstream>
00028 #include "mfem.hpp"
00029
00030 // Exact solution, E, and r.h.s., f. See below for implementation.
00031 void E_exact(const Vector &, Vector &);
00032 void f_exact(const Vector &, Vector &);
00033
00034 int main (int argc, char *argv[])
00035 {
00036    Mesh *mesh;
00037
00038    if (argc == 1)
00039    {
00040       cout << "\nUsage: ex3 <mesh_file>\n" << endl;
00041       return 1;
00042    }
00043
00044    // 1. Read the mesh from the given mesh file. In this 3D example, we can
00045    //    handle tetrahedral or hexahedral meshes with the same code.
00046    ifstream imesh(argv[1]);
00047    if (!imesh)
00048    {
00049       cerr << "\nCan not open mesh file: " << argv[1] << '\n' << endl;
00050       return 2;
00051    }
00052    mesh = new Mesh(imesh, 1, 1);
00053    imesh.close();
00054    if (mesh -> Dimension() != 3)
00055    {
00056       cerr << "\nThis example requires a 3D mesh\n" << endl;
00057       return 3;
00058    }
00059
00060    // 2. Refine the mesh to increase the resolution. In this example we do
00061    //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
00062    //    largest number that gives a final mesh with no more than 50,000
00063    //    elements.
00064    {
00065       int ref_levels =
00066          (int)floor(log(50000./mesh->GetNE())/log(2.)/mesh->Dimension());
00067       for (int l = 0; l < ref_levels; l++)
00068          mesh->UniformRefinement();
00069    }
00070
00071    // 3. Define a finite element space on the mesh. Here we use the lowest order
00072    //    Nedelec finite elements, but we can easily swich to higher-order spaces
00073    //    by changing the value of p.
00074    int p = 1;
00075    FiniteElementCollection *fec = new ND_FECollection(p, mesh -> Dimension());
00076    FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
00077    cout << "Number of unknowns: " << fespace->GetVSize() << endl;
00078
00079    // 4. Set up the linear form b(.) which corresponds to the right-hand side
00080    //    of the FEM linear system, which in this case is (f,phi_i) where f is
00081    //    given by the function f_exact and phi_i are the basis functions in the
00082    //    finite element fespace.
00083    VectorFunctionCoefficient f(3, f_exact);
00084    LinearForm *b = new LinearForm(fespace);
00086    b->Assemble();
00087
00088    // 5. Define the solution vector x as a finite element grid function
00089    //    corresponding to fespace. Initialize x by projecting the exact
00090    //    solution. Note that only values from the boundary edges will be used
00091    //    when eliminating the non-homogenious boundary condition to modify the
00092    //    r.h.s. vector b.
00093    GridFunction x(fespace);
00094    VectorFunctionCoefficient E(3, E_exact);
00095    x.ProjectCoefficient(E);
00096
00097    // 6. Set up the bilinear form corresponding to the EM diffusion operator
00098    //    curl muinv curl + sigma I, by adding the curl-curl and the mass domain
00099    //    integrators and finally imposing the non-homogeneous Dirichlet boundary
00100    //    conditions. The boundary conditions are implemented by marking all the
00101    //    boundary attributes from the mesh as essential (Dirichlet). After
00102    //    assembly and finalizing we extract the corresponding sparse matrix A.
00103    Coefficient *muinv = new ConstantCoefficient(1.0);
00104    Coefficient *sigma = new ConstantCoefficient(1.0);
00105    BilinearForm *a = new BilinearForm(fespace);
00108    a->Assemble();
00109    Array<int> ess_bdr(mesh->bdr_attributes.Max());
00110    ess_bdr = 1;
00111    a->EliminateEssentialBC(ess_bdr, x, *b);
00112    a->Finalize();
00113    const SparseMatrix &A = a->SpMat();
00114
00115    // 7. Define a simple symmetric Gauss-Seidel preconditioner and use it to
00116    //    solve the system Ax=b with PCG.
00117    GSSmoother M(A);
00118    x = 0.0;
00119    PCG(A, M, *b, x, 1, 500, 1e-12, 0.0);
00120
00121    // 8. Compute and print the L^2 norm of the error.
00122    cout << "\n|| E_h - E ||_{L^2} = " << x.ComputeL2Error(E) << '\n' << endl;
00123
00124    // 9. Save the refined mesh and the solution. This output can be viewed
00125    //    later using GLVis: "glvis -m refined.mesh -g sol.gf".
00126    {
00127       ofstream mesh_ofs("refined.mesh");
00128       mesh_ofs.precision(8);
00129       mesh->Print(mesh_ofs);
00130       ofstream sol_ofs("sol.gf");
00131       sol_ofs.precision(8);
00132       x.Save(sol_ofs);
00133    }
00134
00135    // 10. (Optional) Send the solution by socket to a GLVis server.
00136    char vishost[] = "localhost";
00137    int  visport   = 19916;
00138    osockstream sol_sock(visport, vishost);
00139    sol_sock << "solution\n";
00140    sol_sock.precision(8);
00141    mesh->Print(sol_sock);
00142    x.Save(sol_sock);
00143    sol_sock.send();
00144
00145    // 11. Free the used memory.
00146    delete a;
00147    delete sigma;
00148    delete muinv;
00149    delete b;
00150    delete fespace;
00151    delete fec;
00152    delete mesh;
00153
00154    return 0;
00155 }
00156
00157 // A parameter for the exact solution.
00158 const double kappa = M_PI;
00159
00160 void E_exact(const Vector &x, Vector &E)
00161 {
00162    E(0) = sin(kappa * x(1));
00163    E(1) = sin(kappa * x(2));
00164    E(2) = sin(kappa * x(0));
00165 }
00166
00167 void f_exact(const Vector &x, Vector &f)
00168 {
00169    f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
00170    f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
00171    f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
00172 }
```