MFEM v2.0
ex2.cpp
Go to the documentation of this file.
```00001 //                                MFEM Example 2
00002 //
00003 // Compile with: make ex2
00004 //
00005 // Sample runs:  ex2 ../data/beam-tri.mesh
00007 //               ex2 ../data/beam-tet.mesh
00008 //               ex2 ../data/beam-hex.mesh
00010 //               ex2 ../data/beam-hex-nurbs.mesh
00011 //
00012 // Description:  This example code solves a simple linear elasticity problem
00013 //               describing a multi-material cantilever beam.
00014 //
00015 //               Specifically, we approximate the weak form of -div(sigma(u))=0
00017 //               tensor corresponding to displacement field u, and lambda and mu
00018 //               are the material Lame constants. The boundary conditions are
00019 //               u=0 on the fixed part of the boundary with attribute 1, and
00020 //               sigma(u).n=f on the remainder with f being a constant pull down
00021 //               vector on boundary elements with attribute 2, and zero
00022 //               otherwise. The geometry of the domain is assumed to be as
00023 //               follows:
00024 //
00025 //                                 +----------+----------+
00026 //                    boundary --->| material | material |<--- boundary
00027 //                    attribute 1  |    1     |    2     |     attribute 2
00028 //                    (fixed)      +----------+----------+     (pull down)
00029 //
00030 //               The example demonstrates the use of high-order and NURBS vector
00031 //               finite element spaces with the linear elasticity bilinear form,
00032 //               meshes with curved elements, and the definition of piece-wise
00033 //               constant and vector coefficient objects.
00034 //
00035 //               We recommend viewing example 1 before viewing this example.
00036
00037 #include <fstream>
00038 #include "mfem.hpp"
00039
00040 int main (int argc, char *argv[])
00041 {
00042    Mesh *mesh;
00043
00044    if (argc == 1)
00045    {
00046       cout << "\nUsage: ex2 <mesh_file>\n" << endl;
00047       return 1;
00048    }
00049
00050    // 1. Read the mesh from the given mesh file. We can handle triangular,
00051    //    quadrilateral, tetrahedral or hexahedral elements with the same code.
00052    ifstream imesh(argv[1]);
00053    if (!imesh)
00054    {
00055       cerr << "\nCan not open mesh file: " << argv[1] << '\n' << endl;
00056       return 2;
00057    }
00058    mesh = new Mesh(imesh, 1, 1);
00059    imesh.close();
00060
00061    if (mesh->attributes.Max() < 2 || mesh->bdr_attributes.Max() < 2)
00062    {
00063       cerr << "\nInput mesh should have at least two materials and "
00064            << "two boundary attributes! (See schematic in ex2.cpp)\n"
00065            << endl;
00066       return 3;
00067    }
00068
00069    int dim = mesh->Dimension();
00070
00071    // 2. Select the order of the finite element discretization space. For NURBS
00072    //    meshes, we increase the order by degree elevation.
00073    int p;
00074    cout << "Enter finite element space order --> " << flush;
00075    cin >> p;
00076
00077    if (mesh->NURBSext && p > mesh->NURBSext->GetOrder())
00078       mesh->DegreeElevate(p - mesh->NURBSext->GetOrder());
00079
00080    // 3. Refine the mesh to increase the resolution. In this example we do
00081    //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
00082    //    largest number that gives a final mesh with no more than 5,000
00083    //    elements.
00084    {
00085       int ref_levels =
00086          (int)floor(log(5000./mesh->GetNE())/log(2.)/dim);
00087       for (int l = 0; l < ref_levels; l++)
00088          mesh->UniformRefinement();
00089    }
00090
00091    // 4. Define a finite element space on the mesh. Here we use vector finite
00092    //    elements, i.e. dim copies of a scalar finite element space. The vector
00093    //    dimension is specified by the last argument of the FiniteElementSpace
00094    //    constructor. For NURBS meshes, we use the (degree elevated) NURBS space
00095    //    associated with the mesh nodes.
00096    FiniteElementCollection *fec;
00097    FiniteElementSpace *fespace;
00098    if (mesh->NURBSext)
00099    {
00100       fec = NULL;
00101       fespace = mesh->GetNodes()->FESpace();
00102    }
00103    else
00104    {
00105       fec = new H1_FECollection(p, dim);
00106       fespace = new FiniteElementSpace(mesh, fec, dim);
00107    }
00108    cout << "Number of unknowns: " << fespace->GetVSize() << endl
00109         << "Assembling: " << flush;
00110
00111    // 5. Set up the linear form b(.) which corresponds to the right-hand side of
00112    //    the FEM linear system. In this case, b_i equals the boundary integral
00113    //    of f*phi_i where f represents a "pull down" force on the Neumann part
00114    //    of the boundary and phi_i are the basis functions in the finite element
00115    //    fespace. The force is defined by the VectorArrayCoefficient object f,
00116    //    which is a vector of Coefficient objects. The fact that f is non-zero
00117    //    on boundary attribute 2 is indicated by the use of piece-wise constants
00118    //    coefficient for its last component.
00119    VectorArrayCoefficient f(dim);
00120    for (int i = 0; i < dim-1; i++)
00121       f.Set(i, new ConstantCoefficient(0.0));
00122    {
00123       Vector pull_force(mesh->bdr_attributes.Max());
00124       pull_force = 0.0;
00125       pull_force(1) = -1.0e-2;
00126       f.Set(dim-1, new PWConstCoefficient(pull_force));
00127    }
00128
00129    LinearForm *b = new LinearForm(fespace);
00131    cout << "r.h.s. ... " << flush;
00132    b->Assemble();
00133
00134    // 6. Define the solution vector x as a finite element grid function
00135    //    corresponding to fespace. Initialize x with initial guess of zero,
00136    //    which satisfies the boundary conditions.
00137    GridFunction x(fespace);
00138    x = 0.0;
00139
00140    // 7. Set up the bilinear form a(.,.) on the finite element space
00141    //    corresponding to the linear elasticity integrator with piece-wise
00142    //    constants coefficient lambda and mu. The boundary conditions are
00143    //    implemented by marking only boundary attribute 1 as essential. After
00144    //    assembly and finalizing we extract the corresponding sparse matrix A.
00145    Vector lambda(mesh->attributes.Max());
00146    lambda = 1.0;
00147    lambda(0) = lambda(1)*50;
00148    PWConstCoefficient lambda_func(lambda);
00149    Vector mu(mesh->attributes.Max());
00150    mu = 1.0;
00151    mu(0) = mu(1)*50;
00152    PWConstCoefficient mu_func(mu);
00153
00154    BilinearForm *a = new BilinearForm(fespace);
00156    cout << "matrix ... " << flush;
00157    a->Assemble();
00158    Array<int> ess_bdr(mesh->bdr_attributes.Max());
00159    ess_bdr = 0;
00160    ess_bdr[0] = 1;
00161    a->EliminateEssentialBC(ess_bdr, x, *b);
00162    a->Finalize();
00163    cout << "done." << endl;
00164    const SparseMatrix &A = a->SpMat();
00165
00166    // 8. Define a simple symmetric Gauss-Seidel preconditioner and use it to
00167    //    solve the system Ax=b with PCG.
00168    GSSmoother M(A);
00169    PCG(A, M, *b, x, 1, 500, 1e-8, 0.0);
00170
00171    // 9. For non-NURBS meshes, make the mesh curved based on the finite element
00172    //    space. This means that we define the mesh elements through a fespace
00173    //    based transformation of the reference element. This allows us to save
00174    //    the displaced mesh as a curved mesh when using high-order finite
00175    //    element displacement field. We assume that the initial mesh (read from
00176    //    the file) is not higher order curved mesh compared to the chosen FE
00177    //    space.
00178    if (!mesh->NURBSext)
00179       mesh->SetNodalFESpace(fespace);
00180
00181    // 10. Save the displaced mesh and the inverted solution (which gives the
00182    //     backward displacements to the original grid). This output can be
00183    //     viewed later using GLVis: "glvis -m displaced.mesh -g sol.gf".
00184    {
00185       GridFunction *nodes = mesh->GetNodes();
00186       *nodes += x;
00187       x *= -1;
00188       ofstream mesh_ofs("displaced.mesh");
00189       mesh_ofs.precision(8);
00190       mesh->Print(mesh_ofs);
00191       ofstream sol_ofs("sol.gf");
00192       sol_ofs.precision(8);
00193       x.Save(sol_ofs);
00194    }
00195
00196    // 11. (Optional) Send the above data by socket to a GLVis server. Use the
00197    //     "n" and "b" keys in GLVis to visualize the displacements.
00198    char vishost[] = "localhost";
00199    int  visport   = 19916;
00200    osockstream sol_sock(visport, vishost);
00201    sol_sock << "solution\n";
00202    sol_sock.precision(8);
00203    mesh->Print(sol_sock);
00204    x.Save(sol_sock);
00205    sol_sock.send();
00206
00207    // 12. Free the used memory.
00208    delete a;
00209    delete b;
00210    if (fec)
00211    {
00212       delete fespace;
00213       delete fec;
00214    }
00215    delete mesh;
00216
00217    return 0;
00218 }
```