bp1.cpp (9b8727523add6a096aae88596ec6827f5594f30d) bp1.cpp (dc00e230834c021cd7c5b556bbaba87dc4209676)
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

--- 37 unchanged lines hidden (view full) ---

46 const char *ceed_spec = "/cpu/self";
47#ifndef MFEM_DIR
48 const char *mesh_file = "../../../mfem/data/star.mesh";
49#else
50 const char *mesh_file = MFEM_DIR "/data/star.mesh";
51#endif
52 int order = 1;
53 bool visualization = true;
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

--- 37 unchanged lines hidden (view full) ---

46 const char *ceed_spec = "/cpu/self";
47#ifndef MFEM_DIR
48 const char *mesh_file = "../../../mfem/data/star.mesh";
49#else
50 const char *mesh_file = MFEM_DIR "/data/star.mesh";
51#endif
52 int order = 1;
53 bool visualization = true;
54 bool test = false;
54
55 mfem::OptionsParser args(argc, argv);
56 args.AddOption(&ceed_spec, "-c", "-ceed", "Ceed specification.");
57 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
58 args.AddOption(&order, "-o", "--order",
59 "Finite element order (polynomial degree).");
60 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
61 "--no-visualization",
62 "Enable or disable GLVis visualization.");
55
56 mfem::OptionsParser args(argc, argv);
57 args.AddOption(&ceed_spec, "-c", "-ceed", "Ceed specification.");
58 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
59 args.AddOption(&order, "-o", "--order",
60 "Finite element order (polynomial degree).");
61 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
62 "--no-visualization",
63 "Enable or disable GLVis visualization.");
64 args.AddOption(&test, "-t", "--test", "-no-test",
65 "--no-test",
66 "Enable or disable test mode.");
63 args.Parse();
64 if (!args.Good()) {
65 args.PrintUsage(std::cout);
66 return 1;
67 }
67 args.Parse();
68 if (!args.Good()) {
69 args.PrintUsage(std::cout);
70 return 1;
71 }
68 args.PrintOptions(std::cout);
72 if (!test) {
73 args.PrintOptions(std::cout);
74 }
69
70 // 2. Initialize a Ceed device object using the given Ceed specification.
71 Ceed ceed;
72 CeedInit(ceed_spec, &ceed);
73
74 // 3. Read the mesh from the given mesh file.
75 mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
76 int dim = mesh->Dimension();

--- 17 unchanged lines hidden (view full) ---

94 mesh->SetCurvature(order, false, -1, mfem::Ordering::byNODES);
95 }
96
97 // 5. Define a finite element space on the mesh. Here we use continuous
98 // Lagrange finite elements of the specified order.
99 MFEM_VERIFY(order > 0, "invalid order");
100 mfem::FiniteElementCollection *fec = new mfem::H1_FECollection(order, dim);
101 mfem::FiniteElementSpace *fespace = new mfem::FiniteElementSpace(mesh, fec);
75
76 // 2. Initialize a Ceed device object using the given Ceed specification.
77 Ceed ceed;
78 CeedInit(ceed_spec, &ceed);
79
80 // 3. Read the mesh from the given mesh file.
81 mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
82 int dim = mesh->Dimension();

--- 17 unchanged lines hidden (view full) ---

100 mesh->SetCurvature(order, false, -1, mfem::Ordering::byNODES);
101 }
102
103 // 5. Define a finite element space on the mesh. Here we use continuous
104 // Lagrange finite elements of the specified order.
105 MFEM_VERIFY(order > 0, "invalid order");
106 mfem::FiniteElementCollection *fec = new mfem::H1_FECollection(order, dim);
107 mfem::FiniteElementSpace *fespace = new mfem::FiniteElementSpace(mesh, fec);
102 std::cout << "Number of finite element unknowns: "
103 << fespace->GetTrueVSize() << std::endl;
108 if (!test) {
109 std::cout << "Number of finite element unknowns: "
110 << fespace->GetTrueVSize() << std::endl;
111 }
104
105 // 6. Construct a rhs vector using the linear form f(v) = (solution, v), where
106 // v is a test function.
107 mfem::LinearForm b(fespace);
108 mfem::FunctionCoefficient sol_coeff(solution);
109 b.AddDomainIntegrator(new mfem::DomainLFIntegrator(sol_coeff));
110 b.Assemble();
111
112 // 7. Construct a CeedMassOperator utilizing the 'ceed' device and using the
113 // 'fespace' object to extract data needed by the Ceed objects.
114 CeedMassOperator mass(ceed, fespace);
115
116 // 8. Solve the discrete system using the conjugate gradients (CG) method.
117 mfem::CGSolver cg;
118 cg.SetRelTol(1e-6);
119 cg.SetMaxIter(100);
112
113 // 6. Construct a rhs vector using the linear form f(v) = (solution, v), where
114 // v is a test function.
115 mfem::LinearForm b(fespace);
116 mfem::FunctionCoefficient sol_coeff(solution);
117 b.AddDomainIntegrator(new mfem::DomainLFIntegrator(sol_coeff));
118 b.Assemble();
119
120 // 7. Construct a CeedMassOperator utilizing the 'ceed' device and using the
121 // 'fespace' object to extract data needed by the Ceed objects.
122 CeedMassOperator mass(ceed, fespace);
123
124 // 8. Solve the discrete system using the conjugate gradients (CG) method.
125 mfem::CGSolver cg;
126 cg.SetRelTol(1e-6);
127 cg.SetMaxIter(100);
120 cg.SetPrintLevel(3);
128 if (test) {
129 cg.SetPrintLevel(0);
130 } else {
131 cg.SetPrintLevel(3);
132 }
121 cg.SetOperator(mass);
122
123 mfem::GridFunction sol(fespace);
124 sol = 0.0;
125 cg.Mult(b, sol);
126
127 // 9. Compute and print the L2 projection error.
133 cg.SetOperator(mass);
134
135 mfem::GridFunction sol(fespace);
136 sol = 0.0;
137 cg.Mult(b, sol);
138
139 // 9. Compute and print the L2 projection error.
128 std::cout << "L2 projection error: " << sol.ComputeL2Error(sol_coeff)
129 << std::endl;
130 if (fabs(sol.ComputeL2Error(sol_coeff))>1e-4) {
131 std::cout << "Error too large" << std::endl;
132 return 1;
140 if (!test) {
141 std::cout << "L2 projection error: " << sol.ComputeL2Error(sol_coeff)
142 << std::endl;
143 } else {
144 if (fabs(sol.ComputeL2Error(sol_coeff))>1e-4) {
145 std::cout << "Error too large" << std::endl;
146 }
133 }
134
135 // 10. Open a socket connection to GLVis and send the mesh and solution for
136 // visualization.
137 if (visualization) {
138 char vishost[] = "localhost";
139 int visport = 19916;
140 mfem::socketstream sol_sock(vishost, visport);
141 sol_sock.precision(8);
142 sol_sock << "solution\n" << *mesh << sol << std::flush;
143 }
144
145 // 11. Free memory and exit.
146 delete fespace;
147 delete fec;
148 delete mesh;
149 CeedDestroy(&ceed);
150 return 0;
151}
147 }
148
149 // 10. Open a socket connection to GLVis and send the mesh and solution for
150 // visualization.
151 if (visualization) {
152 char vishost[] = "localhost";
153 int visport = 19916;
154 mfem::socketstream sol_sock(vishost, visport);
155 sol_sock.precision(8);
156 sol_sock << "solution\n" << *mesh << sol << std::flush;
157 }
158
159 // 11. Free memory and exit.
160 delete fespace;
161 delete fec;
162 delete mesh;
163 CeedDestroy(&ceed);
164 return 0;
165}