From 45b475f30d6d00822557ee81ea573dbf51a50e24 Mon Sep 17 00:00:00 2001 From: Aneesh Srinivas Date: Thu, 17 Apr 2025 22:45:25 -0700 Subject: [PATCH 1/4] Create biharmonic.cpp --- examples/mfem/biharmonic.cpp | 101 +++++++++++++++++++++++++++++++++++ 1 file changed, 101 insertions(+) create mode 100644 examples/mfem/biharmonic.cpp diff --git a/examples/mfem/biharmonic.cpp b/examples/mfem/biharmonic.cpp new file mode 100644 index 0000000000..60299fc484 --- /dev/null +++ b/examples/mfem/biharmonic.cpp @@ -0,0 +1,101 @@ +#include +#include +#include +#include +#include + +using namespace mfem; +using namespace std; + +// CEED QFunction to apply Laplacian operator +extern "C" CeedInt LaplacianApply(CeedInt Q, + const CeedScalar *const *in, + CeedScalar *const *out, + void *ctx) { + const CeedScalar *grad_u = in[1]; + CeedScalar *v = out[0]; + + const int dim = 2; // 2D example + for (CeedInt i = 0; i < Q; i++) { + v[i] = 0.0; + for (int d = 0; d < dim; d++) { + v[i] += grad_u[i * dim + d]; + } + } + return 0; +} + +// CEED QFunction to build operator coefficients +extern "C" CeedInt LaplacianBuild(CeedInt Q,const CeedScalar *const *in,CeedScalar *const *out, void *ctx) { + CeedScalar *qdata = out[0]; + for (CeedInt i = 0; i < Q; i++) { + qdata[i] = 1.0; // constant coefficient + } + return 0; +} + +int main(int argc, char *argv[]) { + // 1. Initialize MFEM and CEED + Device device("cpu"); + cout << "Running on device: " << device.GetDeviceMemoryType() << endl; + + Ceed ceed; + CeedInit("/cpu/self", &ceed); + + // 2. Mesh + int nx = 8, ny = 8, order = 2; + Mesh mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true); + mesh.EnsureNodes(); + int dim = mesh.Dimension(); + + // 3. Finite element space + H1_FECollection fec(order, dim); + FiniteElementSpace fespace(&mesh, &fec); + cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl; + + // 4. Define GridFunctions + GridFunction rhs(&fespace); + GridFunction w(&fespace); + GridFunction u(&fespace); + + rhs.ProjectCoefficient(ConstantCoefficient(1.0)); // f = 1 + + // 5. Define CEED QFunctions + CeedQFunction qf_build, qf_apply; + + CeedQFunctionCreateInterior(ceed, 1, LaplacianBuild, __FILE__ ":LaplacianBuild", &qf_build); + CeedQFunctionAddInput(qf_build, "dx", CEED_EVAL_GRAD, dim); + CeedQFunctionAddOutput(qf_build, "qdata", CEED_EVAL_NONE, 1); + + CeedQFunctionCreateInterior(ceed, 1, LaplacianApply, __FILE__ ":LaplacianApply", &qf_apply); + CeedQFunctionAddInput(qf_apply, "u", CEED_EVAL_INTERP, 1); + CeedQFunctionAddInput(qf_apply, "grad_u", CEED_EVAL_GRAD, dim); + CeedQFunctionAddInput(qf_apply, "qdata", CEED_EVAL_NONE, 1); + CeedQFunctionAddOutput(qf_apply, "v", CEED_EVAL_INTERP, 1); + + // 6. Create CEED Operator + mfem::ceed::CeedOperatorBuilder builder(&fespace, qf_apply, qf_build, ceed); + CeedOperator op_laplace = builder.GetCeedOperator(); + + // 7. Solver + mfem::ceed::CeedSolver solver(op_laplace); + + // 8. Solve -Δ w = f + cout << "Solving -Δ w = f..." << endl; + solver.Solve(rhs, w); + + // 9. Solve -Δ u = w + cout << "Solving -Δ u = w..." << endl; + solver.Solve(w, u); + + // 10. Output + ofstream ufile("u.gf"), wfile("w.gf"); + u.Save(ufile); + w.Save(wfile); + + cout << "Solutions saved to 'u.gf' and 'w.gf'" << endl; + + // Cleanup + CeedDestroy(&ceed); + return 0; +} \ No newline at end of file From 13fefdaccd8784858eec178e2814282e0d9135fb Mon Sep 17 00:00:00 2001 From: Aneesh Srinivas Date: Fri, 18 Apr 2025 19:36:42 -0700 Subject: [PATCH 2/4] Update biharmonic.cpp --- examples/mfem/biharmonic.cpp | 130 +++++++++++++++++------------------ 1 file changed, 65 insertions(+), 65 deletions(-) diff --git a/examples/mfem/biharmonic.cpp b/examples/mfem/biharmonic.cpp index 60299fc484..36040fb3b1 100644 --- a/examples/mfem/biharmonic.cpp +++ b/examples/mfem/biharmonic.cpp @@ -1,101 +1,101 @@ +// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. +// All Rights Reserved. See the top-level LICENSE and NOTICE files for details. +// +// SPDX-License-Identifier: BSD-2-Clause +// +// This file is part of CEED: http://github.com/ceed + +// libCEED + MFEM Example: BP6 +// +// This example illustrates a simple usage of libCEED with the MFEM (mfem.org) finite element library. +// +// It solves the biharmonic equation \Nabla^2 u =f on Omega, u(\vec{x})=0 on ∂(Omega), du/d(\vec{x})*n=0 on ∂(Omega) by solving -\Nabla w=f, -\Nabla u=w +// +// +// + + +/// @file +/// MFEM mass operator based on libCEED #include #include #include +#include +#include #include #include using namespace mfem; -using namespace std; - -// CEED QFunction to apply Laplacian operator -extern "C" CeedInt LaplacianApply(CeedInt Q, - const CeedScalar *const *in, - CeedScalar *const *out, - void *ctx) { - const CeedScalar *grad_u = in[1]; - CeedScalar *v = out[0]; - - const int dim = 2; // 2D example - for (CeedInt i = 0; i < Q; i++) { - v[i] = 0.0; - for (int d = 0; d < dim; d++) { - v[i] += grad_u[i * dim + d]; - } - } - return 0; -} - -// CEED QFunction to build operator coefficients -extern "C" CeedInt LaplacianBuild(CeedInt Q,const CeedScalar *const *in,CeedScalar *const *out, void *ctx) { - CeedScalar *qdata = out[0]; - for (CeedInt i = 0; i < Q; i++) { - qdata[i] = 1.0; // constant coefficient - } - return 0; -} - -int main(int argc, char *argv[]) { - // 1. Initialize MFEM and CEED - Device device("cpu"); - cout << "Running on device: " << device.GetDeviceMemoryType() << endl; - - Ceed ceed; - CeedInit("/cpu/self", &ceed); - - // 2. Mesh - int nx = 8, ny = 8, order = 2; - Mesh mesh = Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL, true); + +int main(int argc, char *argv[]) +{ + // 1. Initialize MFEM and CEED + Device device("cpu"); + std::cout << "\nRunning on device: " << device.GetDeviceMemoryTypeName() << "\n\n"; + + Ceed ceed; + CeedInit("/cpu/self", &ceed); + + // 2. Read mesh from file + const char *mesh_file = "mesh.mesh"; + Mesh mesh(mesh_file, 1, 1, true); // generate_edges = 1, refine = 1, fix_orientation = true mesh.EnsureNodes(); int dim = mesh.Dimension(); - // 3. Finite element space + + // 3. Set up finite element space H1_FECollection fec(order, dim); FiniteElementSpace fespace(&mesh, &fec); - cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl; - // 4. Define GridFunctions + std::cout << "Number of finite element unknowns: " << fespace.GetTrueVSize() << std::endl; + + // 4. Load source function f from file GridFunction rhs(&fespace); + std::ifstream fin("f.gf"); + if (!fin) + { + std::cerr << "Error: Unable to open file 'f.gf'.\n"; + return 1; + } + rhs.Load(fin); + + // 5. set up intermediate grid functions GridFunction w(&fespace); GridFunction u(&fespace); - rhs.ProjectCoefficient(ConstantCoefficient(1.0)); // f = 1 - - // 5. Define CEED QFunctions + // 6. Use CEED built-in QFunctions for Laplacian CeedQFunction qf_build, qf_apply; - - CeedQFunctionCreateInterior(ceed, 1, LaplacianBuild, __FILE__ ":LaplacianBuild", &qf_build); + CeedQFunctionCreateInterior(ceed, 1, f_build_diff, "f_build_diff", &qf_build); CeedQFunctionAddInput(qf_build, "dx", CEED_EVAL_GRAD, dim); - CeedQFunctionAddOutput(qf_build, "qdata", CEED_EVAL_NONE, 1); + CeedQFunctionAddOutput(qf_build, "qdata", CEED_EVAL_NONE, dim*(dim+1)/2); - CeedQFunctionCreateInterior(ceed, 1, LaplacianApply, __FILE__ ":LaplacianApply", &qf_apply); + CeedQFunctionCreateInterior(ceed, 1, f_apply_diff, "f_apply_diff", &qf_apply); CeedQFunctionAddInput(qf_apply, "u", CEED_EVAL_INTERP, 1); CeedQFunctionAddInput(qf_apply, "grad_u", CEED_EVAL_GRAD, dim); - CeedQFunctionAddInput(qf_apply, "qdata", CEED_EVAL_NONE, 1); + CeedQFunctionAddInput(qf_apply, "qdata", CEED_EVAL_NONE, dim*(dim+1)/2); CeedQFunctionAddOutput(qf_apply, "v", CEED_EVAL_INTERP, 1); - // 6. Create CEED Operator + // 7. Build CEED operator with MFEM wrapper mfem::ceed::CeedOperatorBuilder builder(&fespace, qf_apply, qf_build, ceed); - CeedOperator op_laplace = builder.GetCeedOperator(); + CeedOperator ceed_op = builder.GetCeedOperator(); - // 7. Solver - mfem::ceed::CeedSolver solver(op_laplace); + // 8. Create CEED linear solver + mfem::ceed::CeedSolver solver(ceed_op); - // 8. Solve -Δ w = f - cout << "Solving -Δ w = f..." << endl; + // 9. Solve -Δ w = f and -Δ u = w + std::cout << "Solving -Δ w = f...\n"; solver.Solve(rhs, w); - // 9. Solve -Δ u = w - cout << "Solving -Δ u = w..." << endl; + + std::cout << "Solving -Δ u = w...\n"; solver.Solve(w, u); - // 10. Output - ofstream ufile("u.gf"), wfile("w.gf"); + // 10. Save solutions + std::ofstream ufile("u.gf"); u.Save(ufile); - w.Save(wfile); - - cout << "Solutions saved to 'u.gf' and 'w.gf'" << endl; + std::cout << "Solutions saved to 'u.gf' .\n"; - // Cleanup + // 11. destroy CeedDestroy(&ceed); return 0; } \ No newline at end of file From 4e062d4262d51ea0710397fff9d93759e8dad88b Mon Sep 17 00:00:00 2001 From: Aneesh Srinivas Date: Fri, 2 May 2025 00:52:51 -0700 Subject: [PATCH 3/4] added the .h and .hpp files for biharmonic --- examples/mfem/biharmonic.h | 38 ++++++++++++++++++++++++++++++++++++ examples/mfem/biharmonic.hpp | 37 +++++++++++++++++++++++++++++++++++ 2 files changed, 75 insertions(+) create mode 100644 examples/mfem/biharmonic.h create mode 100644 examples/mfem/biharmonic.hpp diff --git a/examples/mfem/biharmonic.h b/examples/mfem/biharmonic.h new file mode 100644 index 0000000000..f6e42a3cc6 --- /dev/null +++ b/examples/mfem/biharmonic.h @@ -0,0 +1,38 @@ +#ifndef BIHARMONIC_H +#define BIHARMONIC_H + +#include +#include +#include + +namespace biharmonic +{ + // Initialize the MFEM device (e.g., CPU, GPU). + void InitializeDevice(); + + // Initialize the libCEED context with the specified resource. + void InitializeCeed(Ceed &ceed, const std::string &resource = "/cpu/self"); + + // Load the mesh from a file and perform uniform refinements. + mfem::Mesh *LoadMesh(const std::string &mesh_file, int ref_levels = 1); + + // Set up the finite element space with the given polynomial order. + mfem::FiniteElementSpace *SetupFESpace(mfem::Mesh *mesh, int order); + + // Load the right-hand side function 'f' from a file into a GridFunction. + mfem::GridFunction *LoadRHS(const std::string &rhs_file, mfem::FiniteElementSpace *fespace); + + // Set up the CEED QFunctions for the Laplacian operator. + void SetupQFunctions(Ceed ceed, CeedQFunction &qf_build, CeedQFunction &qf_apply, int dim); + + // Build the CEED operator using the MFEM wrapper. + CeedOperator BuildCeedOperator(mfem::FiniteElementSpace *fespace, CeedQFunction qf_apply, CeedQFunction qf_build, Ceed ceed); + + // Solve the equation using the CEED solver. + void Solve(CeedOperator ceed_op, mfem::GridFunction &rhs, mfem::GridFunction &solution); + + // Save the solution to a file. + void SaveSolution(const mfem::GridFunction &solution, const std::string &filename); +} + +#endif // BIHARMONIC_H \ No newline at end of file diff --git a/examples/mfem/biharmonic.hpp b/examples/mfem/biharmonic.hpp new file mode 100644 index 0000000000..0381478666 --- /dev/null +++ b/examples/mfem/biharmonic.hpp @@ -0,0 +1,37 @@ +#ifndef BIHARMONIC_HPP +#define BIHARMONIC_HPP + +#include +#include + +namespace biharmonic +{ + // Initializes the MFEM device. + void InitializeDevice(); + + // Initializes the libCEED context. + void InitializeCeed(Ceed &ceed, const std::string &resource = "/cpu/self"); + + // Loads the mesh from a file and prepares it. + mfem::Mesh *LoadMesh(const std::string &mesh_file, int ref_levels); + + // Sets up the finite element space. + mfem::FiniteElementSpace *SetupFESpace(mfem::Mesh *mesh, int order); + + // Loads the source function 'f' from a file. + mfem::GridFunction *LoadRHS(const std::string &rhs_file, mfem::FiniteElementSpace *fespace); + + // Sets up the CEED QFunctions for the Laplacian. + void SetupQFunctions(Ceed ceed, CeedQFunction &qf_build, CeedQFunction &qf_apply, int dim); + + // Builds the CEED operator using the MFEM wrapper. + CeedOperator BuildCeedOperator(mfem::FiniteElementSpace *fespace, CeedQFunction qf_apply, CeedQFunction qf_build, Ceed ceed); + + // Solves the equation using the CEED solver. + void Solve(CeedOperator ceed_op, mfem::GridFunction &rhs, mfem::GridFunction &solution); + + // Saves the solution to a file. + void SaveSolution(const mfem::GridFunction &solution, const std::string &filename); +} + +#endif // BIHARMONIC_HPP \ No newline at end of file From 038152f21d3b193807ba914fe1bcb8920674f942 Mon Sep 17 00:00:00 2001 From: Aneesh Srinivas Date: Thu, 22 May 2025 12:30:43 -0700 Subject: [PATCH 4/4] edited biharmonic.cpp to be in the style of the other boundary problems --- examples/mfem/Makefile | 2 +- examples/mfem/biharmonic.cpp | 202 ++++++++++++++++++++++------------- examples/mfem/biharmonic.h | 44 +++----- examples/mfem/biharmonic.hpp | 40 ++----- 4 files changed, 154 insertions(+), 134 deletions(-) diff --git a/examples/mfem/Makefile b/examples/mfem/Makefile index cb5abeba01..17cc7cee97 100644 --- a/examples/mfem/Makefile +++ b/examples/mfem/Makefile @@ -24,7 +24,7 @@ MFEM_DEF = -DMFEM_DIR="\"$(abspath $(MFEM_DIR))\"" MFEM_LIB_FILE = mfem_is_not_built -include $(wildcard $(CONFIG_MK)) -MFEM_EXAMPLES = bp1 bp3 +MFEM_EXAMPLES = bp1 bp3 biharmonic .SUFFIXES: .SUFFIXES: .cpp diff --git a/examples/mfem/biharmonic.cpp b/examples/mfem/biharmonic.cpp index 36040fb3b1..34f7ef5add 100644 --- a/examples/mfem/biharmonic.cpp +++ b/examples/mfem/biharmonic.cpp @@ -5,97 +5,155 @@ // // This file is part of CEED: http://github.com/ceed -// libCEED + MFEM Example: BP6 +// libCEED + MFEM Example: BP1 // // This example illustrates a simple usage of libCEED with the MFEM (mfem.org) finite element library. // -// It solves the biharmonic equation \Nabla^2 u =f on Omega, u(\vec{x})=0 on ∂(Omega), du/d(\vec{x})*n=0 on ∂(Omega) by solving -\Nabla w=f, -\Nabla u=w -// -// +// 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 +// 'solution'). The mass matrix required for performing the projection is expressed as a new class, CeedMassOperator, derived from mfem::Operator. +// Internally, CeedMassOperator uses a CeedOperator object constructed based on an mfem::FiniteElementSpace. +// All libCEED objects use a Ceed device object constructed based on a command line argument (-ceed). // - +// The mass matrix is inverted using a simple conjugate gradient algorithm corresponding to CEED BP1, see http://ceed.exascaleproject.org/bps. +// Arbitrary mesh and solution orders in 1D, 2D and 3D are supported from the same code. +// +// Build with: +// +// make biharmonic [MFEM_DIR=] [CEED_DIR=] +// +// Sample runs: +// +// ./biharmonic +// ./biharmonic -ceed /cpu/self +// ./biharmonic -ceed /gpu/cuda +// ./biharmonic -m ../../../mfem/data/fichera.mesh +// ./biharmonic -m ../../../mfem/data/star.vtk -o 3 +// ./biharmonic -m ../../../mfem/data/inline-segment.mesh -o 8 /// @file /// MFEM mass operator based on libCEED -#include -#include -#include -#include -#include -#include -#include -using namespace mfem; +#include "biharmonic.hpp" -int main(int argc, char *argv[]) -{ - // 1. Initialize MFEM and CEED - Device device("cpu"); - std::cout << "\nRunning on device: " << device.GetDeviceMemoryTypeName() << "\n\n"; - - Ceed ceed; - CeedInit("/cpu/self", &ceed); - - // 2. Read mesh from file - const char *mesh_file = "mesh.mesh"; - Mesh mesh(mesh_file, 1, 1, true); // generate_edges = 1, refine = 1, fix_orientation = true - mesh.EnsureNodes(); - int dim = mesh.Dimension(); - - - // 3. Set up finite element space - H1_FECollection fec(order, dim); - FiniteElementSpace fespace(&mesh, &fec); +#include - std::cout << "Number of finite element unknowns: " << fespace.GetTrueVSize() << std::endl; +#include - // 4. Load source function f from file - GridFunction rhs(&fespace); - std::ifstream fin("f.gf"); - if (!fin) - { - std::cerr << "Error: Unable to open file 'f.gf'.\n"; +/// Continuous function to project on the discrete FE space +double solution(const mfem::Vector &pt) { + return pt.Norml2(); // distance to the origin +} + +//TESTARGS -ceed {ceed_resource} -t -no-vis --size 2000 --order 4 +int main(int argc, char *argv[]) { + // 1. Parse command-line options. + const char *ceed_spec = "/cpu/self"; +#ifndef MFEM_DIR + const char *mesh_file = "../../../mfem/data/star.mesh"; +#else + const char *mesh_file = MFEM_DIR "/data/star.mesh"; +#endif + int order = 1; + bool visualization = true; + bool test = false; + double max_nnodes = 50000; + + mfem::OptionsParser args(argc, argv); + args.AddOption(&ceed_spec, "-c", "-ceed", "Ceed specification."); + args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use."); + args.AddOption(&order, "-o", "--order", "Finite element order (polynomial degree)."); + args.AddOption(&max_nnodes, "-s", "--size", "Maximum size (number of DoFs)"); + args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization."); + args.AddOption(&test, "-t", "--test", "-no-test", "--no-test", "Enable or disable test mode."); + args.Parse(); + if (!args.Good()) { + args.PrintUsage(std::cout); return 1; } - rhs.Load(fin); - - // 5. set up intermediate grid functions - GridFunction w(&fespace); - GridFunction u(&fespace); - - // 6. Use CEED built-in QFunctions for Laplacian - CeedQFunction qf_build, qf_apply; - CeedQFunctionCreateInterior(ceed, 1, f_build_diff, "f_build_diff", &qf_build); - CeedQFunctionAddInput(qf_build, "dx", CEED_EVAL_GRAD, dim); - CeedQFunctionAddOutput(qf_build, "qdata", CEED_EVAL_NONE, dim*(dim+1)/2); + if (!test) { + args.PrintOptions(std::cout); + } - CeedQFunctionCreateInterior(ceed, 1, f_apply_diff, "f_apply_diff", &qf_apply); - CeedQFunctionAddInput(qf_apply, "u", CEED_EVAL_INTERP, 1); - CeedQFunctionAddInput(qf_apply, "grad_u", CEED_EVAL_GRAD, dim); - CeedQFunctionAddInput(qf_apply, "qdata", CEED_EVAL_NONE, dim*(dim+1)/2); - CeedQFunctionAddOutput(qf_apply, "v", CEED_EVAL_INTERP, 1); + // 2. Initialize a Ceed device object using the given Ceed specification. + Ceed ceed; + CeedInit(ceed_spec, &ceed); - // 7. Build CEED operator with MFEM wrapper - mfem::ceed::CeedOperatorBuilder builder(&fespace, qf_apply, qf_build, ceed); - CeedOperator ceed_op = builder.GetCeedOperator(); + // 3. Read the mesh from the given mesh file. + mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1); + int dim = mesh->Dimension(); - // 8. Create CEED linear solver - mfem::ceed::CeedSolver solver(ceed_op); + // 4. Refine the mesh to increase the resolution. + // In this example we do 'ref_levels' of uniform refinement. + // We choose 'ref_levels' to be the largest number that gives a final system with no more than 50,000 unknowns, approximately. + { + int ref_levels = (int)floor((log(max_nnodes / mesh->GetNE()) - dim * log(order)) / log(2.) / dim); + for (int l = 0; l < ref_levels; l++) { + mesh->UniformRefinement(); + } + } + if (mesh->GetNodalFESpace() == NULL) { + mesh->SetCurvature(1, false, -1, mfem::Ordering::byNODES); + } + if (mesh->NURBSext) { + mesh->SetCurvature(order, false, -1, mfem::Ordering::byNODES); + } - // 9. Solve -Δ w = f and -Δ u = w - std::cout << "Solving -Δ w = f...\n"; - solver.Solve(rhs, w); + // 5. Define a finite element space on the mesh. + // Here we use continuous Lagrange finite elements of the specified order. + MFEM_VERIFY(order > 0, "invalid order"); + mfem::FiniteElementCollection *fec = new mfem::H1_FECollection(order, dim); + mfem::FiniteElementSpace *fespace = new mfem::FiniteElementSpace(mesh, fec); + if (!test) { + std::cout << "Number of finite element unknowns: " << fespace->GetTrueVSize() << std::endl; + } - - std::cout << "Solving -Δ u = w...\n"; - solver.Solve(w, u); + // 6. Construct a rhs vector using the linear form f(v) = (solution, v), where v is a test function. + mfem::LinearForm b(fespace); + mfem::FunctionCoefficient sol_coeff(solution); + b.AddDomainIntegrator(new mfem::DomainLFIntegrator(sol_coeff)); + b.Assemble(); + + // 7. Construct a CeedMassOperator utilizing the 'ceed' device and using the 'fespace' object to extract data needed by the Ceed objects. + CeedMassOperator mass(ceed, fespace); + + // 8. Solve the discrete system using the conjugate gradients (CG) method. + mfem::CGSolver cg; + cg.SetRelTol(1e-6); + cg.SetMaxIter(100); + if (test) { + cg.SetPrintLevel(0); + } else { + cg.SetPrintLevel(3); + } + cg.SetOperator(mass); + + mfem::GridFunction sol(fespace); + sol = 0.0; + cg.Mult(b, sol); + + // 9. Compute and print the L2 projection error. + double err_l2 = sol.ComputeL2Error(sol_coeff); + if (!test) { + std::cout << "L2 projection error: " << err_l2 << std::endl; + } else { + if (fabs(sol.ComputeL2Error(sol_coeff)) > 2e-4) { + std::cout << "Error too large: " << err_l2 << std::endl; + } + } - // 10. Save solutions - std::ofstream ufile("u.gf"); - u.Save(ufile); - std::cout << "Solutions saved to 'u.gf' .\n"; + // 10. Open a socket connection to GLVis and send the mesh and solution for visualization. + if (visualization) { + char vishost[] = "localhost"; + int visport = 19916; + mfem::socketstream sol_sock(vishost, visport); + sol_sock.precision(8); + sol_sock << "solution\n" << *mesh << sol << std::flush; + } - // 11. destroy + // 11. Free memory and exit. + delete fespace; + delete fec; + delete mesh; CeedDestroy(&ceed); return 0; -} \ No newline at end of file +} diff --git a/examples/mfem/biharmonic.h b/examples/mfem/biharmonic.h index f6e42a3cc6..903a673f40 100644 --- a/examples/mfem/biharmonic.h +++ b/examples/mfem/biharmonic.h @@ -5,34 +5,16 @@ #include #include -namespace biharmonic -{ - // Initialize the MFEM device (e.g., CPU, GPU). - void InitializeDevice(); - - // Initialize the libCEED context with the specified resource. - void InitializeCeed(Ceed &ceed, const std::string &resource = "/cpu/self"); - - // Load the mesh from a file and perform uniform refinements. - mfem::Mesh *LoadMesh(const std::string &mesh_file, int ref_levels = 1); - - // Set up the finite element space with the given polynomial order. - mfem::FiniteElementSpace *SetupFESpace(mfem::Mesh *mesh, int order); - - // Load the right-hand side function 'f' from a file into a GridFunction. - mfem::GridFunction *LoadRHS(const std::string &rhs_file, mfem::FiniteElementSpace *fespace); - - // Set up the CEED QFunctions for the Laplacian operator. - void SetupQFunctions(Ceed ceed, CeedQFunction &qf_build, CeedQFunction &qf_apply, int dim); - - // Build the CEED operator using the MFEM wrapper. - CeedOperator BuildCeedOperator(mfem::FiniteElementSpace *fespace, CeedQFunction qf_apply, CeedQFunction qf_build, Ceed ceed); - - // Solve the equation using the CEED solver. - void Solve(CeedOperator ceed_op, mfem::GridFunction &rhs, mfem::GridFunction &solution); - - // Save the solution to a file. - void SaveSolution(const mfem::GridFunction &solution, const std::string &filename); -} - -#endif // BIHARMONIC_H \ No newline at end of file +namespace biharmonic { + void InitializeDevice(); + void InitializeCeed(Ceed &ceed, const std::string &resource = "/cpu/self"); + mfem::Mesh *LoadMesh(const std::string &mesh_file); + mfem::FiniteElementSpace *SetupFESpace(mfem::Mesh *mesh, int order); + mfem::GridFunction *LoadRHS(const std::string &rhs_file, mfem::FiniteElementSpace *fespace); + void SetupQFunctions(Ceed ceed, CeedQFunction &qf_build, CeedQFunction &qf_apply, int dim); + CeedOperator BuildCeedOperator(mfem::FiniteElementSpace *fespace, CeedQFunction qf_apply, CeedQFunction qf_build, Ceed ceed); + void Solve(CeedOperator ceed_op, mfem::GridFunction &rhs, mfem::GridFunction &solution); + void SaveSolution(const mfem::GridFunction &solution, const std::string &filename); +} // namespace biharmonic + +#endif // BIHARMONIC_H diff --git a/examples/mfem/biharmonic.hpp b/examples/mfem/biharmonic.hpp index 0381478666..eabb827029 100644 --- a/examples/mfem/biharmonic.hpp +++ b/examples/mfem/biharmonic.hpp @@ -1,37 +1,17 @@ #ifndef BIHARMONIC_HPP #define BIHARMONIC_HPP -#include +#include "biharmonic.h" #include -namespace biharmonic -{ - // Initializes the MFEM device. - void InitializeDevice(); +/// CEED QFunction for building quadrature data for diffusion +CEED_QFUNCTION(f_build_diff)(void *ctx, CeedInt Q, + const CeedScalar *const *in, + CeedScalar *const *out); - // Initializes the libCEED context. - void InitializeCeed(Ceed &ceed, const std::string &resource = "/cpu/self"); +/// CEED QFunction for applying diffusion operator +CEED_QFUNCTION(f_apply_diff)(void *ctx, CeedInt Q, + const CeedScalar *const *in, + CeedScalar *const *out); - // Loads the mesh from a file and prepares it. - mfem::Mesh *LoadMesh(const std::string &mesh_file, int ref_levels); - - // Sets up the finite element space. - mfem::FiniteElementSpace *SetupFESpace(mfem::Mesh *mesh, int order); - - // Loads the source function 'f' from a file. - mfem::GridFunction *LoadRHS(const std::string &rhs_file, mfem::FiniteElementSpace *fespace); - - // Sets up the CEED QFunctions for the Laplacian. - void SetupQFunctions(Ceed ceed, CeedQFunction &qf_build, CeedQFunction &qf_apply, int dim); - - // Builds the CEED operator using the MFEM wrapper. - CeedOperator BuildCeedOperator(mfem::FiniteElementSpace *fespace, CeedQFunction qf_apply, CeedQFunction qf_build, Ceed ceed); - - // Solves the equation using the CEED solver. - void Solve(CeedOperator ceed_op, mfem::GridFunction &rhs, mfem::GridFunction &solution); - - // Saves the solution to a file. - void SaveSolution(const mfem::GridFunction &solution, const std::string &filename); -} - -#endif // BIHARMONIC_HPP \ No newline at end of file +#endif // BIHARMONIC_HPP