1 // libCEED + MFEM Example: BP1 2 // 3 // This example illustrates a simple usage of libCEED with the MFEM (mfem.org) 4 // finite element library. 5 // 6 // The example reads a mesh from a file and solves a simple linear system with a 7 // mass matrix (L2-projection of a given analytic function provided by 8 // 'solution'). The mass matrix required for performing the projection is 9 // expressed as a new class, CeedMassOperator, derived from mfem::Operator. 10 // Internally, CeedMassOperator uses a CeedOperator object constructed based on 11 // an mfem::FiniteElementSpace. All libCEED objects use a Ceed device object 12 // constructed based on a command line argument (-ceed). 13 // 14 // The mass matrix is inverted using a simple conjugate gradient algorithm 15 // corresponding to CEED BP1, see http://ceed.exascaleproject.org/bps. Arbitrary 16 // mesh and solution orders in 1D, 2D and 3D are supported from the same code. 17 // 18 // Build with: 19 // 20 // make bp1 [MFEM_DIR=</path/to/mfem>] [CEED_DIR=</path/to/libceed>] 21 // 22 // Sample runs: 23 // 24 // ./bp1 25 // ./bp1 -ceed /cpu/self 26 // ./bp1 -ceed /gpu/occa 27 // ./bp1 -ceed /cpu/occa 28 // ./bp1 -ceed /omp/occa 29 // ./bp1 -ceed /ocl/occa 30 // ./bp1 -m ../../../mfem/data/fichera.mesh 31 // ./bp1 -m ../../../mfem/data/star.vtk -o 3 32 // ./bp1 -m ../../../mfem/data/inline-segment.mesh -o 8 33 34 /// @file 35 /// MFEM mass operator based on libCEED 36 37 #include <ceed.h> 38 #include <mfem.hpp> 39 #include "bp1.hpp" 40 41 /// Continuous function to project on the discrete FE space 42 double solution(const mfem::Vector &pt) { 43 return pt.Norml2(); // distance to the origin 44 } 45 46 //TESTARGS -ceed {ceed_resource} -t -no-vis --size 2000 47 int main(int argc, char *argv[]) { 48 // 1. Parse command-line options. 49 const char *ceed_spec = "/cpu/self"; 50 #ifndef MFEM_DIR 51 const char *mesh_file = "../../../mfem/data/star.mesh"; 52 #else 53 const char *mesh_file = MFEM_DIR "/data/star.mesh"; 54 #endif 55 int order = 1; 56 bool visualization = true; 57 bool test = false; 58 double max_nnodes = 50000; 59 60 mfem::OptionsParser args(argc, argv); 61 args.AddOption(&ceed_spec, "-c", "-ceed", "Ceed specification."); 62 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use."); 63 args.AddOption(&order, "-o", "--order", 64 "Finite element order (polynomial degree)."); 65 args.AddOption(&max_nnodes, "-s", "--size", "Maximum size (number of DoFs)"); 66 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", 67 "--no-visualization", 68 "Enable or disable GLVis visualization."); 69 args.AddOption(&test, "-t", "--test", "-no-test", 70 "--no-test", 71 "Enable or disable test mode."); 72 args.Parse(); 73 if (!args.Good()) { 74 args.PrintUsage(std::cout); 75 return 1; 76 } 77 if (!test) { 78 args.PrintOptions(std::cout); 79 } 80 81 // 2. Initialize a Ceed device object using the given Ceed specification. 82 Ceed ceed; 83 CeedInit(ceed_spec, &ceed); 84 85 // 3. Read the mesh from the given mesh file. 86 mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1); 87 int dim = mesh->Dimension(); 88 89 // 4. Refine the mesh to increase the resolution. In this example we do 90 // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the 91 // largest number that gives a final system with no more than 50,000 92 // unknowns, approximately. 93 { 94 int ref_levels = 95 (int)floor((log(max_nnodes/mesh->GetNE())-dim*log(order))/log(2.)/dim); 96 for (int l = 0; l < ref_levels; l++) { 97 mesh->UniformRefinement(); 98 } 99 } 100 if (mesh->GetNodalFESpace() == NULL) { 101 mesh->SetCurvature(1, false, -1, mfem::Ordering::byNODES); 102 } 103 if (mesh->NURBSext) { 104 mesh->SetCurvature(order, false, -1, mfem::Ordering::byNODES); 105 } 106 107 // 5. Define a finite element space on the mesh. Here we use continuous 108 // Lagrange finite elements of the specified order. 109 MFEM_VERIFY(order > 0, "invalid order"); 110 mfem::FiniteElementCollection *fec = new mfem::H1_FECollection(order, dim); 111 mfem::FiniteElementSpace *fespace = new mfem::FiniteElementSpace(mesh, fec); 112 if (!test) { 113 std::cout << "Number of finite element unknowns: " 114 << fespace->GetTrueVSize() << std::endl; 115 } 116 117 // 6. Construct a rhs vector using the linear form f(v) = (solution, v), where 118 // v is a test function. 119 mfem::LinearForm b(fespace); 120 mfem::FunctionCoefficient sol_coeff(solution); 121 b.AddDomainIntegrator(new mfem::DomainLFIntegrator(sol_coeff)); 122 b.Assemble(); 123 124 // 7. Construct a CeedMassOperator utilizing the 'ceed' device and using the 125 // 'fespace' object to extract data needed by the Ceed objects. 126 CeedMassOperator mass(ceed, fespace); 127 128 // 8. Solve the discrete system using the conjugate gradients (CG) method. 129 mfem::CGSolver cg; 130 cg.SetRelTol(1e-6); 131 cg.SetMaxIter(100); 132 if (test) { 133 cg.SetPrintLevel(0); 134 } else { 135 cg.SetPrintLevel(3); 136 } 137 cg.SetOperator(mass); 138 139 mfem::GridFunction sol(fespace); 140 sol = 0.0; 141 cg.Mult(b, sol); 142 143 // 9. Compute and print the L2 projection error. 144 if (!test) { 145 std::cout << "L2 projection error: " << sol.ComputeL2Error(sol_coeff) 146 << std::endl; 147 } else { 148 if (fabs(sol.ComputeL2Error(sol_coeff))>2e-4) { 149 std::cout << "Error too large" << std::endl; 150 } 151 } 152 153 // 10. Open a socket connection to GLVis and send the mesh and solution for 154 // visualization. 155 if (visualization) { 156 char vishost[] = "localhost"; 157 int visport = 19916; 158 mfem::socketstream sol_sock(vishost, visport); 159 sol_sock.precision(8); 160 sol_sock << "solution\n" << *mesh << sol << std::flush; 161 } 162 163 // 11. Free memory and exit. 164 delete fespace; 165 delete fec; 166 delete mesh; 167 CeedDestroy(&ceed); 168 return 0; 169 } 170