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