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