bp3.cpp (5dab8e2af977e36360f79b2790b71cc75431f50b) bp3.cpp (dc00e230834c021cd7c5b556bbaba87dc4209676)
1// libCEED + MFEM Example: BP3
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 linear system with a
7// diffusion stiffness matrix (with a prescribed analytic solution, provided by
8// the function 'solution'). The diffusion matrix is expressed as a new class,

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

65 const char *ceed_spec = "/cpu/self";
66#ifndef MFEM_DIR
67 const char *mesh_file = "../../../mfem/data/star.mesh";
68#else
69 const char *mesh_file = MFEM_DIR "/data/star.mesh";
70#endif
71 int order = 2;
72 bool visualization = true;
1// libCEED + MFEM Example: BP3
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 linear system with a
7// diffusion stiffness matrix (with a prescribed analytic solution, provided by
8// the function 'solution'). The diffusion matrix is expressed as a new class,

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

65 const char *ceed_spec = "/cpu/self";
66#ifndef MFEM_DIR
67 const char *mesh_file = "../../../mfem/data/star.mesh";
68#else
69 const char *mesh_file = MFEM_DIR "/data/star.mesh";
70#endif
71 int order = 2;
72 bool visualization = true;
73 bool test = false;
73
74 mfem::OptionsParser args(argc, argv);
75 args.AddOption(&ceed_spec, "-c", "-ceed", "Ceed specification.");
76 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
77 args.AddOption(&order, "-o", "--order",
78 "Finite element order (polynomial degree).");
79 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
80 "--no-visualization",
81 "Enable or disable GLVis visualization.");
74
75 mfem::OptionsParser args(argc, argv);
76 args.AddOption(&ceed_spec, "-c", "-ceed", "Ceed specification.");
77 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
78 args.AddOption(&order, "-o", "--order",
79 "Finite element order (polynomial degree).");
80 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
81 "--no-visualization",
82 "Enable or disable GLVis visualization.");
83 args.AddOption(&test, "-t", "--test", "-no-test",
84 "--no-test",
85 "Enable or disable test mode.");
82 args.Parse();
83 if (!args.Good()) {
84 args.PrintUsage(std::cout);
85 return 1;
86 }
86 args.Parse();
87 if (!args.Good()) {
88 args.PrintUsage(std::cout);
89 return 1;
90 }
87 args.PrintOptions(std::cout);
91 if (!test) {
92 args.PrintOptions(std::cout);
93 }
88
89 // 2. Initialize a Ceed device object using the given Ceed specification.
90 Ceed ceed;
91 CeedInit(ceed_spec, &ceed);
92
93 // 3. Read the mesh from the given mesh file.
94 mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
95 int dim = mesh->Dimension();

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

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

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

119 mesh->SetCurvature(order, false, -1, mfem::Ordering::byNODES);
120 }
121
122 // 5. Define a finite element space on the mesh. Here we use continuous
123 // Lagrange finite elements of the specified order.
124 MFEM_VERIFY(order > 0, "invalid order");
125 mfem::FiniteElementCollection *fec = new mfem::H1_FECollection(order, dim);
126 mfem::FiniteElementSpace *fespace = new mfem::FiniteElementSpace(mesh, fec);
121 std::cout << "Number of finite element unknowns: "
122 << fespace->GetTrueVSize() << std::endl;
127 if (!test) {
128 std::cout << "Number of finite element unknowns: "
129 << fespace->GetTrueVSize() << std::endl;
130 }
123
124 mfem::FunctionCoefficient sol_coeff(solution);
125 mfem::Array<int> ess_tdof_list;
126 mfem::GridFunction sol(fespace);
127 if (mesh->bdr_attributes.Size()) {
128 mfem::Array<int> ess_bdr(mesh->bdr_attributes.Max());
129 ess_bdr = 1;
130 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);

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

145 mfem::Operator *D;
146 mfem::Vector X, B;
147 diff.FormLinearSystem(ess_tdof_list, sol, b, D, X, B);
148
149 // 8. Solve the discrete system using the conjugate gradients (CG) method.
150 mfem::CGSolver cg;
151 cg.SetRelTol(1e-6);
152 cg.SetMaxIter(1000);
131
132 mfem::FunctionCoefficient sol_coeff(solution);
133 mfem::Array<int> ess_tdof_list;
134 mfem::GridFunction sol(fespace);
135 if (mesh->bdr_attributes.Size()) {
136 mfem::Array<int> ess_bdr(mesh->bdr_attributes.Max());
137 ess_bdr = 1;
138 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);

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

153 mfem::Operator *D;
154 mfem::Vector X, B;
155 diff.FormLinearSystem(ess_tdof_list, sol, b, D, X, B);
156
157 // 8. Solve the discrete system using the conjugate gradients (CG) method.
158 mfem::CGSolver cg;
159 cg.SetRelTol(1e-6);
160 cg.SetMaxIter(1000);
153 cg.SetPrintLevel(3);
161 if (test) {
162 cg.SetPrintLevel(0);
163 } else {
164 cg.SetPrintLevel(3);
165 }
154 cg.SetOperator(*D);
155
156 cg.Mult(B, X);
157
158 // 9. Compute and print the L2 norm of the error.
166 cg.SetOperator(*D);
167
168 cg.Mult(B, X);
169
170 // 9. Compute and print the L2 norm of the error.
159 std::cout << "L2 norm of the error: " << sol.ComputeL2Error(sol_coeff)
160 << std::endl;
171 if (!test) {
172 std::cout << "L2 projection error: " << sol.ComputeL2Error(sol_coeff)
173 << std::endl;
174 } else {
175 if (fabs(sol.ComputeL2Error(sol_coeff))>1e-4) {
176 std::cout << "Error too large" << std::endl;
177 }
178 }
161
162 // 10. Open a socket connection to GLVis and send the mesh and solution for
163 // visualization.
164 if (visualization) {
165 char vishost[] = "localhost";
166 int visport = 19916;
167 mfem::socketstream sol_sock(vishost, visport);
168 sol_sock.precision(8);
169 sol_sock << "solution\n" << *mesh << sol << std::flush;
170 }
171
172 // 11. Free memory and exit.
173 delete fespace;
174 delete fec;
175 delete mesh;
176 CeedDestroy(&ceed);
177 return 0;
178}
179
180 // 10. Open a socket connection to GLVis and send the mesh and solution for
181 // visualization.
182 if (visualization) {
183 char vishost[] = "localhost";
184 int visport = 19916;
185 mfem::socketstream sol_sock(vishost, visport);
186 sol_sock.precision(8);
187 sol_sock << "solution\n" << *mesh << sol << std::flush;
188 }
189
190 // 11. Free memory and exit.
191 delete fespace;
192 delete fec;
193 delete mesh;
194 CeedDestroy(&ceed);
195 return 0;
196}