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
solution(const mfem::Vector & pt)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
main(int argc,char * argv[])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