1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "mfem.hpp"
13 #include "unit_tests.hpp"
14
15 using namespace mfem;
16
17 namespace get_value
18 {
19
func_1D_lin(const Vector & x)20 double func_1D_lin(const Vector &x)
21 {
22 return x[0];
23 }
24
func_2D_lin(const Vector & x)25 double func_2D_lin(const Vector &x)
26 {
27 return x[0] + 2.0 * x[1];
28 }
29
func_3D_lin(const Vector & x)30 double func_3D_lin(const Vector &x)
31 {
32 return x[0] + 2.0 * x[1] + 3.0 * x[2];
33 }
34
Func_2D_lin(const Vector & x,Vector & v)35 void Func_2D_lin(const Vector &x, Vector &v)
36 {
37 v.SetSize(2);
38 v[0] = 1.234 * x[0] - 2.357 * x[1];
39 v[1] = 2.537 * x[0] + 4.321 * x[1];
40 }
41
Func_3D_lin(const Vector & x,Vector & v)42 void Func_3D_lin(const Vector &x, Vector &v)
43 {
44 v.SetSize(3);
45 v[0] = 1.234 * x[0] - 2.357 * x[1] + 3.572 * x[2];
46 v[1] = 2.537 * x[0] + 4.321 * x[1] - 1.234 * x[2];
47 v[2] = -2.572 * x[0] + 1.321 * x[1] + 3.234 * x[2];
48 }
49
func_1D_quad(const Vector & x)50 double func_1D_quad(const Vector &x)
51 {
52 return 2.0 * x[0] + x[0] * x[0];
53 }
54
dfunc_1D_quad(const Vector & x,Vector & v)55 void dfunc_1D_quad(const Vector &x, Vector &v)
56 {
57 v.SetSize(1);
58 v[0] = 2.0 + 2.0 * x[0];
59 }
60
func_2D_quad(const Vector & x)61 double func_2D_quad(const Vector &x)
62 {
63 return x[0] * x[0] + 2.0 * x[1] * x[1] + 3.0 * x[0] * x[1];
64 }
65
dfunc_2D_quad(const Vector & x,Vector & v)66 void dfunc_2D_quad(const Vector &x, Vector &v)
67 {
68 v.SetSize(2);
69 v[0] = 2.0 * x[0] + 3.0 * x[1];
70 v[1] = 4.0 * x[1] + 3.0 * x[0];
71 }
72
Func_2D_quad(const Vector & x,Vector & v)73 void Func_2D_quad(const Vector &x, Vector &v)
74 {
75 v.SetSize(2);
76 v[0] = 1.0 * x[0] * x[0] + 2.0 * x[1] * x[1] + 3.0 * x[0] * x[1];
77 v[1] = 2.0 * x[1] * x[1] + 3.0 * x[0] * x[0] + 1.0 * x[0] * x[1];
78 }
79
RotFunc_2D_quad(const Vector & x,Vector & v)80 void RotFunc_2D_quad(const Vector &x, Vector &v)
81 {
82 v.SetSize(1);
83 v[0] = 6.0 * x[0] + 1.0 * x[1] - 4.0 * x[1] - 3.0 * x[0];
84 }
85
DivFunc_2D_quad(const Vector & x)86 double DivFunc_2D_quad(const Vector &x)
87 {
88 return 3.0 * x[0] + 7.0 * x[1];
89 }
90
func_3D_quad(const Vector & x)91 double func_3D_quad(const Vector &x)
92 {
93 return x[0] * x[1] + 2.0 * x[1] * x[2] + 3.0 * x[2] * x[0];
94 }
95
dfunc_3D_quad(const Vector & x,Vector & v)96 void dfunc_3D_quad(const Vector &x, Vector &v)
97 {
98 v.SetSize(3);
99 v[0] = 1.0 * x[1] + 3.0 * x[2];
100 v[1] = 2.0 * x[2] + 1.0 * x[0];
101 v[2] = 3.0 * x[0] + 2.0 * x[1];
102 }
103
Func_3D_quad(const Vector & x,Vector & v)104 void Func_3D_quad(const Vector &x, Vector &v)
105 {
106 v.SetSize(3);
107 v[0] = 1.0 * x[0] * x[1] + 2.0 * x[1] * x[2] + 3.0 * x[2] * x[0];
108 v[1] = 2.0 * x[1] * x[2] + 3.0 * x[2] * x[0] + 1.0 * x[0] * x[1];
109 v[2] = 3.0 * x[2] * x[0] + 1.0 * x[0] * x[1] + 2.0 * x[1] * x[2];
110 }
111
CurlFunc_3D_quad(const Vector & x,Vector & v)112 void CurlFunc_3D_quad(const Vector &x, Vector &v)
113 {
114 v.SetSize(3);
115 v[0] = 1.0 * x[0] + 2.0 * x[2] - 2.0 * x[1] - 3.0 * x[0];
116 v[1] = 2.0 * x[1] + 3.0 * x[0] - 3.0 * x[2] - 1.0 * x[1];
117 v[2] = 3.0 * x[2] + 1.0 * x[1] - 1.0 * x[0] - 2.0 * x[2];
118 }
119
DivFunc_3D_quad(const Vector & x)120 double DivFunc_3D_quad(const Vector &x)
121 {
122 return 4.0 * x[0] + 3.0 * x[1] + 5.0 * x[2];
123 }
124
125 TEST_CASE("1D GetValue",
126 "[GridFunction]"
127 "[GridFunctionCoefficient]")
128 {
129 int log = 1;
130 int n = 1;
131 int dim = 1;
132 int order = 1;
133 int npts = 0;
134
135 double tol = 1e-6;
136
137 for (int type = (int)Element::SEGMENT;
138 type <= (int)Element::SEGMENT; type++)
139 {
140 Mesh mesh = Mesh::MakeCartesian1D(n, 2.0);
141
142 FunctionCoefficient funcCoef(func_1D_lin);
143
to_string(type)144 SECTION("1D GetValue tests for element type " + std::to_string(type))
145 {
146 H1_FECollection h1_fec(order, dim);
147 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
148 FiniteElement::VALUE);
149 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
150 FiniteElement::INTEGRAL);
151
152 FiniteElementSpace h1_fespace(&mesh, &h1_fec);
153 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
154 FiniteElementSpace dgi_fespace(&mesh, &dgi_fec);
155
156 GridFunction h1_x(&h1_fespace);
157 GridFunction dgv_x(&dgv_fespace);
158 GridFunction dgi_x(&dgi_fespace);
159
160 GridFunctionCoefficient h1_xCoef(&h1_x);
161 GridFunctionCoefficient dgv_xCoef(&dgv_x);
162 GridFunctionCoefficient dgi_xCoef(&dgi_x);
163
164 h1_x.ProjectCoefficient(funcCoef);
165 dgv_x.ProjectCoefficient(funcCoef);
166 dgi_x.ProjectCoefficient(funcCoef);
167
168 SECTION("Domain Evaluation 1D")
169 {
170 std::cout << "Domain Evaluation 1D" << std::endl;
171 for (int e = 0; e < mesh.GetNE(); e++)
172 {
173 ElementTransformation *T = mesh.GetElementTransformation(e);
174 const FiniteElement *fe = h1_fespace.GetFE(e);
175 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
176 2*order + 2);
177
178 double h1_gfc_err = 0.0;
179 double dgv_gfc_err = 0.0;
180 double dgi_gfc_err = 0.0;
181
182 double h1_gv_err = 0.0;
183 double dgv_gv_err = 0.0;
184 double dgi_gv_err = 0.0;
185
186 for (int j=0; j<ir.GetNPoints(); j++)
187 {
188 npts++;
189 const IntegrationPoint &ip = ir.IntPoint(j);
190 T->SetIntPoint(&ip);
191
192 double f_val = funcCoef.Eval(*T, ip);
193
194 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
195 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
196 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
197
198 double h1_gv_val = h1_x.GetValue(e, ip);
199 double dgv_gv_val = dgv_x.GetValue(e, ip);
200 double dgi_gv_val = dgi_x.GetValue(e, ip);
201
202 h1_gfc_err += fabs(f_val - h1_gfc_val);
203 dgv_gfc_err += fabs(f_val - dgv_gfc_val);
204 dgi_gfc_err += fabs(f_val - dgi_gfc_val);
205
206 h1_gv_err += fabs(f_val - h1_gv_val);
207 dgv_gv_err += fabs(f_val - dgv_gv_val);
208 dgi_gv_err += fabs(f_val - dgi_gv_val);
209
210 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
211 {
212 std::cout << e << ":" << j << " h1 gfc " << f_val << " "
213 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
214 << std::endl;
215 }
216 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
217 {
218 std::cout << e << ":" << j << " dgv gfc " << f_val << " "
219 << dgv_gfc_val << " "
220 << fabs(f_val - dgv_gfc_val)
221 << std::endl;
222 }
223 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
224 {
225 std::cout << e << ":" << j << " dgi gfc " << f_val << " "
226 << dgi_gfc_val << " "
227 << fabs(f_val - dgi_gfc_val)
228 << std::endl;
229 }
230 if (log > 0 && fabs(f_val - h1_gv_val) > tol)
231 {
232 std::cout << e << ":" << j << " h1 gv " << f_val << " "
233 << h1_gv_val << " " << fabs(f_val - h1_gv_val)
234 << std::endl;
235 }
236 if (log > 0 && fabs(f_val - dgv_gv_val) > tol)
237 {
238 std::cout << e << ":" << j << " dgv gv " << f_val << " "
239 << dgv_gv_val << " "
240 << fabs(f_val - dgv_gv_val)
241 << std::endl;
242 }
243 if (log > 0 && fabs(f_val - dgi_gv_val) > tol)
244 {
245 std::cout << e << ":" << j << " dgi gv " << f_val << " "
246 << dgi_gv_val << " "
247 << fabs(f_val - dgi_gv_val)
248 << std::endl;
249 }
250 }
251
252 h1_gfc_err /= ir.GetNPoints();
253 dgv_gfc_err /= ir.GetNPoints();
254 dgi_gfc_err /= ir.GetNPoints();
255
256 h1_gv_err /= ir.GetNPoints();
257 dgv_gv_err /= ir.GetNPoints();
258 dgi_gv_err /= ir.GetNPoints();
259
260 REQUIRE(h1_gfc_err == MFEM_Approx(0.0));
261 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
262 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
263
264 REQUIRE(h1_gv_err == MFEM_Approx(0.0));
265 REQUIRE(dgv_gv_err == MFEM_Approx(0.0));
266 REQUIRE(dgi_gv_err == MFEM_Approx(0.0));
267 }
268 }
269
270 SECTION("Boundary Evaluation 1D (H1 Context)")
271 {
272 std::cout << "Boundary Evaluation 1D (H1 Context)" << std::endl;
273 for (int be = 0; be < mesh.GetNBE(); be++)
274 {
275 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
276 const FiniteElement *fe = h1_fespace.GetBE(be);
277 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
278 2*order + 2);
279
280 double h1_err = 0.0;
281 double dgv_err = 0.0;
282 double dgi_err = 0.0;
283
284 for (int j=0; j<ir.GetNPoints(); j++)
285 {
286 npts++;
287 const IntegrationPoint &ip = ir.IntPoint(j);
288 T->SetIntPoint(&ip);
289
290 double f_val = funcCoef.Eval(*T, ip);
291 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
292 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
293 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
294
295 h1_err += fabs(f_val - h1_gfc_val);
296 dgv_err += fabs(f_val - dgv_gfc_val);
297 dgi_err += fabs(f_val - dgi_gfc_val);
298
299 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
300 {
301 std::cout << be << ":" << j << " h1 " << f_val << " "
302 << h1_gfc_val << " "
303 << fabs(f_val - h1_gfc_val)
304 << std::endl;
305 }
306 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
307 {
308 std::cout << be << ":" << j << " dgv " << f_val << " "
309 << dgv_gfc_val << " "
310 << fabs(f_val - dgv_gfc_val)
311 << std::endl;
312 }
313 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
314 {
315 std::cout << be << ":" << j << " dgi " << f_val << " "
316 << dgi_gfc_val << " "
317 << fabs(f_val - dgi_gfc_val)
318 << std::endl;
319 }
320 }
321 h1_err /= ir.GetNPoints();
322 dgv_err /= ir.GetNPoints();
323 dgi_err /= ir.GetNPoints();
324
325 REQUIRE(h1_err == MFEM_Approx(0.0));
326 REQUIRE(dgv_err == MFEM_Approx(0.0));
327 REQUIRE(dgi_err == MFEM_Approx(0.0));
328 }
329 }
330
331 SECTION("Boundary Evaluation 1D (DG Context)")
332 {
333 std::cout << "Boundary Evaluation 1D (DG Context)" << std::endl;
334 for (int be = 0; be < mesh.GetNBE(); be++)
335 {
336 FaceElementTransformations *T =
337 mesh.GetBdrFaceTransformations(be);
338 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
339 2*order + 2);
340
341 double h1_err = 0.0;
342 double dgv_err = 0.0;
343 double dgi_err = 0.0;
344
345 for (int j=0; j<ir.GetNPoints(); j++)
346 {
347 npts++;
348 const IntegrationPoint &ip = ir.IntPoint(j);
349
350 T->SetIntPoint(&ip);
351
352 double f_val = funcCoef.Eval(*T, ip);
353 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
354 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
355 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
356
357 h1_err += fabs(f_val - h1_gfc_val);
358 dgv_err += fabs(f_val - dgv_gfc_val);
359 dgi_err += fabs(f_val - dgi_gfc_val);
360
361 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
362 {
363 std::cout << be << ":" << j << " h1 " << f_val << " "
364 << h1_gfc_val << " "
365 << fabs(f_val - h1_gfc_val)
366 << std::endl;
367 }
368 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
369 {
370 std::cout << be << ":" << j << " dgv " << f_val << " "
371 << dgv_gfc_val << " "
372 << fabs(f_val - dgv_gfc_val)
373 << std::endl;
374 }
375 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
376 {
377 std::cout << be << ":" << j << " dgi " << f_val << " "
378 << dgi_gfc_val << " "
379 << fabs(f_val - dgi_gfc_val)
380 << std::endl;
381 }
382 }
383 h1_err /= ir.GetNPoints();
384 dgv_err /= ir.GetNPoints();
385 dgi_err /= ir.GetNPoints();
386
387 REQUIRE(h1_err == MFEM_Approx(0.0));
388 REQUIRE(dgv_err == MFEM_Approx(0.0));
389 REQUIRE(dgi_err == MFEM_Approx(0.0));
390 }
391 }
392 }
393 }
394 std::cout << "Checked GridFunction::GetValue at "
395 << npts << " 1D points" << std::endl;
396 }
397
398 #ifdef MFEM_USE_MPI
399 #
400 TEST_CASE("1D GetValue in Parallel",
401 "[ParGridFunction]"
402 "[GridFunctionCoefficient]"
403 "[Parallel]")
404 {
405 int num_procs;
406 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
407
408 int my_rank;
409 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
410
411 int log = 1;
412 int n = 2 * num_procs;
413 int dim = 1;
414 int order = 1;
415 int npts = 0;
416
417 double tol = 1e-6;
418
419 for (int type = (int)Element::SEGMENT;
420 type <= (int)Element::SEGMENT; type++)
421 {
422 Mesh mesh = Mesh::MakeCartesian1D(n, 2.0);
423 ParMesh pmesh(MPI_COMM_WORLD, mesh);
424 mesh.Clear();
425
426 FunctionCoefficient funcCoef(func_1D_lin);
427
to_string(type)428 SECTION("1D GetValue tests for element type " + std::to_string(type))
429 {
430 H1_FECollection h1_fec(order, dim);
431 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
432 FiniteElement::VALUE);
433 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
434 FiniteElement::INTEGRAL);
435
436 ParFiniteElementSpace h1_fespace(&pmesh, &h1_fec);
437 ParFiniteElementSpace dgv_fespace(&pmesh, &dgv_fec);
438 ParFiniteElementSpace dgi_fespace(&pmesh, &dgi_fec);
439
440 ParGridFunction h1_x(&h1_fespace);
441 ParGridFunction dgv_x(&dgv_fespace);
442 ParGridFunction dgi_x(&dgi_fespace);
443
444 GridFunctionCoefficient h1_xCoef(&h1_x);
445 GridFunctionCoefficient dgv_xCoef(&dgv_x);
446 GridFunctionCoefficient dgi_xCoef(&dgi_x);
447
448 h1_x.ProjectCoefficient(funcCoef);
449 dgv_x.ProjectCoefficient(funcCoef);
450 dgi_x.ProjectCoefficient(funcCoef);
451
452 h1_x.ExchangeFaceNbrData();
453 dgv_x.ExchangeFaceNbrData();
454 dgi_x.ExchangeFaceNbrData();
455
456 SECTION("Shared Face Evaluation 1D")
457 {
458 if (my_rank == 0)
459 {
460 std::cout << "Shared Face Evaluation 1D" << std::endl;
461 }
462 for (int sf = 0; sf < pmesh.GetNSharedFaces(); sf++)
463 {
464 FaceElementTransformations *FET =
465 pmesh.GetSharedFaceTransformations(sf);
466 ElementTransformation *T = &FET->GetElement2Transformation();
467 int e = FET->Elem2No;
468 int e_nbr = e - pmesh.GetNE();
469 const FiniteElement *fe = dgv_fespace.GetFaceNbrFE(e_nbr);
470 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
471 2*order + 2);
472
473 double h1_gfc_err = 0.0;
474 double dgv_gfc_err = 0.0;
475 double dgi_gfc_err = 0.0;
476
477 double h1_gv_err = 0.0;
478 double dgv_gv_err = 0.0;
479 double dgi_gv_err = 0.0;
480
481 for (int j=0; j<ir.GetNPoints(); j++)
482 {
483 npts++;
484 const IntegrationPoint &ip = ir.IntPoint(j);
485 T->SetIntPoint(&ip);
486
487 double f_val = funcCoef.Eval(*T, ip);
488
489 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
490 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
491 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
492
493 double h1_gv_val = h1_x.GetValue(e, ip);
494 double dgv_gv_val = dgv_x.GetValue(e, ip);
495 double dgi_gv_val = dgi_x.GetValue(e, ip);
496
497 h1_gfc_err += fabs(f_val - h1_gfc_val);
498 dgv_gfc_err += fabs(f_val - dgv_gfc_val);
499 dgi_gfc_err += fabs(f_val - dgi_gfc_val);
500
501 h1_gv_err += fabs(f_val - h1_gv_val);
502 dgv_gv_err += fabs(f_val - dgv_gv_val);
503 dgi_gv_err += fabs(f_val - dgi_gv_val);
504
505 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
506 {
507 std::cout << e << ":" << j << " h1 gfc " << f_val << " "
508 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
509 << std::endl;
510 }
511 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
512 {
513 std::cout << e << ":" << j << " dgv gfc " << f_val << " "
514 << dgv_gfc_val << " "
515 << fabs(f_val - dgv_gfc_val)
516 << std::endl;
517 }
518 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
519 {
520 std::cout << e << ":" << j << " dgi gfc " << f_val << " "
521 << dgi_gfc_val << " "
522 << fabs(f_val - dgi_gfc_val)
523 << std::endl;
524 }
525 if (log > 0 && fabs(f_val - h1_gv_val) > tol)
526 {
527 std::cout << e << ":" << j << " h1 gv " << f_val << " "
528 << h1_gv_val << " " << fabs(f_val - h1_gv_val)
529 << std::endl;
530 }
531 if (log > 0 && fabs(f_val - dgv_gv_val) > tol)
532 {
533 std::cout << e << ":" << j << " dgv gv " << f_val << " "
534 << dgv_gv_val << " "
535 << fabs(f_val - dgv_gv_val)
536 << std::endl;
537 }
538 if (log > 0 && fabs(f_val - dgi_gv_val) > tol)
539 {
540 std::cout << e << ":" << j << " dgi gv " << f_val << " "
541 << dgi_gv_val << " "
542 << fabs(f_val - dgi_gv_val)
543 << std::endl;
544 }
545 }
546
547 h1_gfc_err /= ir.GetNPoints();
548 dgv_gfc_err /= ir.GetNPoints();
549 dgi_gfc_err /= ir.GetNPoints();
550
551 h1_gv_err /= ir.GetNPoints();
552 dgv_gv_err /= ir.GetNPoints();
553 dgi_gv_err /= ir.GetNPoints();
554
555 REQUIRE(h1_gfc_err == MFEM_Approx(0.0));
556 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
557 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
558
559 REQUIRE(h1_gv_err == MFEM_Approx(0.0));
560 REQUIRE(dgv_gv_err == MFEM_Approx(0.0));
561 REQUIRE(dgi_gv_err == MFEM_Approx(0.0));
562 }
563 }
564 }
565 }
566 std::cout << my_rank << ": Checked GridFunction::GetValue at "
567 << npts << " 1D points" << std::endl;
568 }
569
570 #endif // MFEM_USE_MPI
571
572 TEST_CASE("2D GetValue",
573 "[GridFunction]"
574 "[GridFunctionCoefficient]")
575 {
576 int log = 1;
577 int n = 1;
578 int dim = 2;
579 int order = 1;
580 int npts = 0;
581
582 double tol = 1e-6;
583
584 for (int type = (int)Element::TRIANGLE;
585 type <= (int)Element::QUADRILATERAL; type++)
586 {
587 Mesh mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
588
589 FunctionCoefficient funcCoef(func_2D_lin);
590
to_string(type)591 SECTION("2D GetValue tests for element type " + std::to_string(type))
592 {
593 H1_FECollection h1_fec(order, dim);
594 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
595 FiniteElement::VALUE);
596 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
597 FiniteElement::INTEGRAL);
598
599 FiniteElementSpace h1_fespace(&mesh, &h1_fec);
600 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
601 FiniteElementSpace dgi_fespace(&mesh, &dgi_fec);
602
603 GridFunction h1_x(&h1_fespace);
604 GridFunction dgv_x(&dgv_fespace);
605 GridFunction dgi_x(&dgi_fespace);
606
607 GridFunctionCoefficient h1_xCoef(&h1_x);
608 GridFunctionCoefficient dgv_xCoef(&dgv_x);
609 GridFunctionCoefficient dgi_xCoef(&dgi_x);
610
611 h1_x.ProjectCoefficient(funcCoef);
612 dgv_x.ProjectCoefficient(funcCoef);
613 dgi_x.ProjectCoefficient(funcCoef);
614
615 SECTION("Domain Evaluation 2D")
616 {
617 std::cout << "Domain Evaluation 2D" << std::endl;
618 for (int e = 0; e < mesh.GetNE(); e++)
619 {
620 ElementTransformation *T = mesh.GetElementTransformation(e);
621 const FiniteElement *fe = h1_fespace.GetFE(e);
622 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
623 2*order + 2);
624
625 double h1_gfc_err = 0.0;
626 double dgv_gfc_err = 0.0;
627 double dgi_gfc_err = 0.0;
628
629 double h1_gv_err = 0.0;
630 double dgv_gv_err = 0.0;
631 double dgi_gv_err = 0.0;
632
633 for (int j=0; j<ir.GetNPoints(); j++)
634 {
635 npts++;
636 const IntegrationPoint &ip = ir.IntPoint(j);
637 T->SetIntPoint(&ip);
638
639 double f_val = funcCoef.Eval(*T, ip);
640
641 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
642 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
643 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
644
645 double h1_gv_val = h1_x.GetValue(e, ip);
646 double dgv_gv_val = dgv_x.GetValue(e, ip);
647 double dgi_gv_val = dgi_x.GetValue(e, ip);
648
649 h1_gfc_err += fabs(f_val - h1_gfc_val);
650 dgv_gfc_err += fabs(f_val - dgv_gfc_val);
651 dgi_gfc_err += fabs(f_val - dgi_gfc_val);
652
653 h1_gv_err += fabs(f_val - h1_gv_val);
654 dgv_gv_err += fabs(f_val - dgv_gv_val);
655 dgi_gv_err += fabs(f_val - dgi_gv_val);
656
657 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
658 {
659 std::cout << e << ":" << j << " h1 gfc " << f_val << " "
660 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
661 << std::endl;
662 }
663 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
664 {
665 std::cout << e << ":" << j << " dgv gfc " << f_val << " "
666 << dgv_gfc_val << " "
667 << fabs(f_val - dgv_gfc_val)
668 << std::endl;
669 }
670 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
671 {
672 std::cout << e << ":" << j << " dgi gfc " << f_val << " "
673 << dgi_gfc_val << " "
674 << fabs(f_val - dgi_gfc_val)
675 << std::endl;
676 }
677 if (log > 0 && fabs(f_val - h1_gv_val) > tol)
678 {
679 std::cout << e << ":" << j << " h1 gv " << f_val << " "
680 << h1_gv_val << " " << fabs(f_val - h1_gv_val)
681 << std::endl;
682 }
683 if (log > 0 && fabs(f_val - dgv_gv_val) > tol)
684 {
685 std::cout << e << ":" << j << " dgv gv " << f_val << " "
686 << dgv_gv_val << " "
687 << fabs(f_val - dgv_gv_val)
688 << std::endl;
689 }
690 if (log > 0 && fabs(f_val - dgi_gv_val) > tol)
691 {
692 std::cout << e << ":" << j << " dgi gv " << f_val << " "
693 << dgi_gv_val << " "
694 << fabs(f_val - dgi_gv_val)
695 << std::endl;
696 }
697 }
698
699 h1_gfc_err /= ir.GetNPoints();
700 dgv_gfc_err /= ir.GetNPoints();
701 dgi_gfc_err /= ir.GetNPoints();
702
703 h1_gv_err /= ir.GetNPoints();
704 dgv_gv_err /= ir.GetNPoints();
705 dgi_gv_err /= ir.GetNPoints();
706
707 REQUIRE(h1_gfc_err == MFEM_Approx(0.0));
708 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
709 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
710
711 REQUIRE(h1_gv_err == MFEM_Approx(0.0));
712 REQUIRE(dgv_gv_err == MFEM_Approx(0.0));
713 REQUIRE(dgi_gv_err == MFEM_Approx(0.0));
714 }
715 }
716
717 SECTION("Boundary Evaluation 2D (H1 Context)")
718 {
719 std::cout << "Boundary Evaluation 2D (H1 Context)" << std::endl;
720 for (int be = 0; be < mesh.GetNBE(); be++)
721 {
722 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
723 const FiniteElement *fe = h1_fespace.GetBE(be);
724 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
725 2*order + 2);
726
727 double h1_err = 0.0;
728 double dgv_err = 0.0;
729 double dgi_err = 0.0;
730
731 for (int j=0; j<ir.GetNPoints(); j++)
732 {
733 npts++;
734 const IntegrationPoint &ip = ir.IntPoint(j);
735 T->SetIntPoint(&ip);
736
737 double f_val = funcCoef.Eval(*T, ip);
738 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
739 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
740 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
741
742 h1_err += fabs(f_val - h1_gfc_val);
743 dgv_err += fabs(f_val - dgv_gfc_val);
744 dgi_err += fabs(f_val - dgi_gfc_val);
745
746 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
747 {
748 std::cout << be << ":" << j << " h1 " << f_val << " "
749 << h1_gfc_val << " "
750 << fabs(f_val - h1_gfc_val)
751 << std::endl;
752 }
753 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
754 {
755 std::cout << be << ":" << j << " dgv " << f_val << " "
756 << dgv_gfc_val << " "
757 << fabs(f_val - dgv_gfc_val)
758 << std::endl;
759 }
760 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
761 {
762 std::cout << be << ":" << j << " dgi " << f_val << " "
763 << dgi_gfc_val << " "
764 << fabs(f_val - dgi_gfc_val)
765 << std::endl;
766 }
767 }
768 h1_err /= ir.GetNPoints();
769 dgv_err /= ir.GetNPoints();
770 dgi_err /= ir.GetNPoints();
771
772 REQUIRE(h1_err == MFEM_Approx(0.0));
773 REQUIRE(dgv_err == MFEM_Approx(0.0));
774 REQUIRE(dgi_err == MFEM_Approx(0.0));
775 }
776 }
777
778 SECTION("Boundary Evaluation 2D (DG Context)")
779 {
780 std::cout << "Boundary Evaluation 2D (DG Context)" << std::endl;
781 for (int be = 0; be < mesh.GetNBE(); be++)
782 {
783 FaceElementTransformations *T =
784 mesh.GetBdrFaceTransformations(be);
785 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
786 2*order + 2);
787
788 double h1_err = 0.0;
789 double dgv_err = 0.0;
790 double dgi_err = 0.0;
791
792 for (int j=0; j<ir.GetNPoints(); j++)
793 {
794 npts++;
795 const IntegrationPoint &ip = ir.IntPoint(j);
796
797 T->SetIntPoint(&ip);
798
799 double f_val = funcCoef.Eval(*T, ip);
800 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
801 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
802 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
803
804 h1_err += fabs(f_val - h1_gfc_val);
805 dgv_err += fabs(f_val - dgv_gfc_val);
806 dgi_err += fabs(f_val - dgi_gfc_val);
807
808 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
809 {
810 std::cout << be << ":" << j << " h1 " << f_val << " "
811 << h1_gfc_val << " "
812 << fabs(f_val - h1_gfc_val)
813 << std::endl;
814 }
815 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
816 {
817 std::cout << be << ":" << j << " dgv " << f_val << " "
818 << dgv_gfc_val << " "
819 << fabs(f_val - dgv_gfc_val)
820 << std::endl;
821 }
822 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
823 {
824 std::cout << be << ":" << j << " dgi " << f_val << " "
825 << dgi_gfc_val << " "
826 << fabs(f_val - dgi_gfc_val)
827 << std::endl;
828 }
829 }
830 h1_err /= ir.GetNPoints();
831 dgv_err /= ir.GetNPoints();
832 dgi_err /= ir.GetNPoints();
833
834 REQUIRE(h1_err == MFEM_Approx(0.0));
835 REQUIRE(dgv_err == MFEM_Approx(0.0));
836 REQUIRE(dgi_err == MFEM_Approx(0.0));
837 }
838 }
839
840 SECTION("Edge Evaluation 2D (H1 Context)")
841 {
842 std::cout << "Edge Evaluation 2D (H1 Context)" << std::endl;
843 for (int e = 0; e < mesh.GetNEdges(); e++)
844 {
845 ElementTransformation *T = mesh.GetEdgeTransformation(e);
846 const FiniteElement *fe = h1_fespace.GetEdgeElement(e);
847 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
848 2*order + 2);
849
850 double h1_err = 0.0;
851
852 for (int j=0; j<ir.GetNPoints(); j++)
853 {
854 npts++;
855 const IntegrationPoint &ip = ir.IntPoint(j);
856 T->SetIntPoint(&ip);
857
858 double f_val = funcCoef.Eval(*T, ip);
859 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
860
861 h1_err += fabs(f_val - h1_gfc_val);
862
863 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
864 {
865 std::cout << e << ":" << j << " h1 " << f_val << " "
866 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
867 << std::endl;
868 }
869 }
870 h1_err /= ir.GetNPoints();
871
872 REQUIRE(h1_err == MFEM_Approx(0.0));
873 }
874 }
875 }
876 }
877 std::cout << "Checked GridFunction::GetValue at "
878 << npts << " 2D points" << std::endl;
879 }
880
881 #ifdef MFEM_USE_MPI
882 #
883 TEST_CASE("2D GetValue in Parallel",
884 "[ParGridFunction]"
885 "[GridFunctionCoefficient]"
886 "[Parallel]")
887 {
888 int num_procs;
889 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
890
891 int my_rank;
892 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
893
894 int log = 1;
895 int n = (int)ceil(sqrt(2*num_procs));
896 int dim = 2;
897 int order = 1;
898 int npts = 0;
899
900 double tol = 1e-6;
901
902 for (int type = (int)Element::TRIANGLE;
903 type <= (int)Element::QUADRILATERAL; type++)
904 {
905 Mesh mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
906 ParMesh pmesh(MPI_COMM_WORLD, mesh);
907 mesh.Clear();
908
909 FunctionCoefficient funcCoef(func_2D_lin);
910
to_string(type)911 SECTION("2D GetValue tests for element type " + std::to_string(type))
912 {
913 H1_FECollection h1_fec(order, dim);
914 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
915 FiniteElement::VALUE);
916 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
917 FiniteElement::INTEGRAL);
918
919 ParFiniteElementSpace h1_fespace(&pmesh, &h1_fec);
920 ParFiniteElementSpace dgv_fespace(&pmesh, &dgv_fec);
921 ParFiniteElementSpace dgi_fespace(&pmesh, &dgi_fec);
922
923 ParGridFunction h1_x(&h1_fespace);
924 ParGridFunction dgv_x(&dgv_fespace);
925 ParGridFunction dgi_x(&dgi_fespace);
926
927 GridFunctionCoefficient h1_xCoef(&h1_x);
928 GridFunctionCoefficient dgv_xCoef(&dgv_x);
929 GridFunctionCoefficient dgi_xCoef(&dgi_x);
930
931 h1_x.ProjectCoefficient(funcCoef);
932 dgv_x.ProjectCoefficient(funcCoef);
933 dgi_x.ProjectCoefficient(funcCoef);
934
935 h1_x.ExchangeFaceNbrData();
936 dgv_x.ExchangeFaceNbrData();
937 dgi_x.ExchangeFaceNbrData();
938
939 SECTION("Shared Face Evaluation 2D")
940 {
941 if (my_rank == 0)
942 {
943 std::cout << "Shared Face Evaluation 2D" << std::endl;
944 }
945 for (int sf = 0; sf < pmesh.GetNSharedFaces(); sf++)
946 {
947 FaceElementTransformations *FET =
948 pmesh.GetSharedFaceTransformations(sf);
949 ElementTransformation *T = &FET->GetElement2Transformation();
950 int e = FET->Elem2No;
951 int e_nbr = e - pmesh.GetNE();
952 const FiniteElement *fe = dgv_fespace.GetFaceNbrFE(e_nbr);
953 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
954 2*order + 2);
955
956 double h1_gfc_err = 0.0;
957 double dgv_gfc_err = 0.0;
958 double dgi_gfc_err = 0.0;
959
960 double h1_gv_err = 0.0;
961 double dgv_gv_err = 0.0;
962 double dgi_gv_err = 0.0;
963
964 for (int j=0; j<ir.GetNPoints(); j++)
965 {
966 npts++;
967 const IntegrationPoint &ip = ir.IntPoint(j);
968 T->SetIntPoint(&ip);
969
970 double f_val = funcCoef.Eval(*T, ip);
971 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
972 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
973 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
974
975 double h1_gv_val = h1_x.GetValue(e, ip);
976 double dgv_gv_val = dgv_x.GetValue(e, ip);
977 double dgi_gv_val = dgi_x.GetValue(e, ip);
978
979 h1_gfc_err += fabs(f_val - h1_gfc_val);
980 dgv_gfc_err += fabs(f_val - dgv_gfc_val);
981 dgi_gfc_err += fabs(f_val - dgi_gfc_val);
982
983 h1_gv_err += fabs(f_val - h1_gv_val);
984 dgv_gv_err += fabs(f_val - dgv_gv_val);
985 dgi_gv_err += fabs(f_val - dgi_gv_val);
986
987 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
988 {
989 std::cout << e << ":" << j << " h1 gfc " << f_val << " "
990 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
991 << std::endl;
992 }
993 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
994 {
995 std::cout << e << ":" << j << " dgv gfc " << f_val << " "
996 << dgv_gfc_val << " "
997 << fabs(f_val - dgv_gfc_val)
998 << std::endl;
999 }
1000 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
1001 {
1002 std::cout << e << ":" << j << " dgi gfc " << f_val << " "
1003 << dgi_gfc_val << " "
1004 << fabs(f_val - dgi_gfc_val)
1005 << std::endl;
1006 }
1007 if (log > 0 && fabs(f_val - h1_gv_val) > tol)
1008 {
1009 std::cout << e << ":" << j << " h1 gv " << f_val << " "
1010 << h1_gv_val << " " << fabs(f_val - h1_gv_val)
1011 << std::endl;
1012 }
1013 if (log > 0 && fabs(f_val - dgv_gv_val) > tol)
1014 {
1015 std::cout << e << ":" << j << " dgv gv " << f_val << " "
1016 << dgv_gv_val << " "
1017 << fabs(f_val - dgv_gv_val)
1018 << std::endl;
1019 }
1020 if (log > 0 && fabs(f_val - dgi_gv_val) > tol)
1021 {
1022 std::cout << e << ":" << j << " dgi gv " << f_val << " "
1023 << dgi_gv_val << " "
1024 << fabs(f_val - dgi_gv_val)
1025 << std::endl;
1026 }
1027 }
1028
1029 h1_gfc_err /= ir.GetNPoints();
1030 dgv_gfc_err /= ir.GetNPoints();
1031 dgi_gfc_err /= ir.GetNPoints();
1032
1033 h1_gv_err /= ir.GetNPoints();
1034 dgv_gv_err /= ir.GetNPoints();
1035 dgi_gv_err /= ir.GetNPoints();
1036
1037 REQUIRE(h1_gfc_err == MFEM_Approx(0.0));
1038 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
1039 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
1040
1041 REQUIRE(h1_gv_err == MFEM_Approx(0.0));
1042 REQUIRE(dgv_gv_err == MFEM_Approx(0.0));
1043 REQUIRE(dgi_gv_err == MFEM_Approx(0.0));
1044 }
1045 }
1046 }
1047 }
1048 std::cout << my_rank << ": Checked GridFunction::GetValue at "
1049 << npts << " 2D points" << std::endl;
1050 }
1051
1052 #endif // MFEM_USE_MPI
1053
1054 TEST_CASE("3D GetValue",
1055 "[GridFunction]"
1056 "[GridFunctionCoefficient]")
1057 {
1058 int log = 1;
1059 int n = 1;
1060 int dim = 3;
1061 int order = 1;
1062 int npts = 0;
1063
1064 double tol = 1e-6;
1065
1066 for (int type = (int)Element::TETRAHEDRON;
1067 type <= (int)Element::WEDGE; type++)
1068 {
1069 Mesh mesh = Mesh::MakeCartesian3D(
1070 n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
1071
1072 FunctionCoefficient funcCoef(func_3D_lin);
1073
to_string(type)1074 SECTION("3D GetValue tests for element type " + std::to_string(type))
1075 {
1076 H1_FECollection h1_fec(order, dim);
1077 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
1078 FiniteElement::VALUE);
1079 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
1080 FiniteElement::INTEGRAL);
1081
1082 FiniteElementSpace h1_fespace(&mesh, &h1_fec);
1083 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
1084 FiniteElementSpace dgi_fespace(&mesh, &dgi_fec);
1085
1086 GridFunction h1_x(&h1_fespace);
1087 GridFunction dgv_x(&dgv_fespace);
1088 GridFunction dgi_x(&dgi_fespace);
1089
1090 GridFunctionCoefficient h1_xCoef(&h1_x);
1091 GridFunctionCoefficient dgv_xCoef(&dgv_x);
1092 GridFunctionCoefficient dgi_xCoef(&dgi_x);
1093
1094 h1_x.ProjectCoefficient(funcCoef);
1095 dgv_x.ProjectCoefficient(funcCoef);
1096 dgi_x.ProjectCoefficient(funcCoef);
1097
1098 SECTION("Domain Evaluation 3D")
1099 {
1100 std::cout << "Domain Evaluation 3D" << std::endl;
1101 for (int e = 0; e < mesh.GetNE(); e++)
1102 {
1103 ElementTransformation *T = mesh.GetElementTransformation(e);
1104 const FiniteElement *fe = h1_fespace.GetFE(e);
1105 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
1106 2*order + 2);
1107
1108 double h1_gfc_err = 0.0;
1109 double dgv_gfc_err = 0.0;
1110 double dgi_gfc_err = 0.0;
1111
1112 double h1_gv_err = 0.0;
1113 double dgv_gv_err = 0.0;
1114 double dgi_gv_err = 0.0;
1115
1116 for (int j=0; j<ir.GetNPoints(); j++)
1117 {
1118 npts++;
1119 const IntegrationPoint &ip = ir.IntPoint(j);
1120 T->SetIntPoint(&ip);
1121
1122 double f_val = funcCoef.Eval(*T, ip);
1123
1124 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
1125 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
1126 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
1127
1128 double h1_gv_val = h1_x.GetValue(e, ip);
1129 double dgv_gv_val = dgv_x.GetValue(e, ip);
1130 double dgi_gv_val = dgi_x.GetValue(e, ip);
1131
1132 h1_gfc_err += fabs(f_val - h1_gfc_val);
1133 dgv_gfc_err += fabs(f_val - dgv_gfc_val);
1134 dgi_gfc_err += fabs(f_val - dgi_gfc_val);
1135
1136 h1_gv_err += fabs(f_val - h1_gv_val);
1137 dgv_gv_err += fabs(f_val - dgv_gv_val);
1138 dgi_gv_err += fabs(f_val - dgi_gv_val);
1139
1140 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
1141 {
1142 std::cout << e << ":" << j << " h1 gfc " << f_val << " "
1143 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
1144 << std::endl;
1145 }
1146 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
1147 {
1148 std::cout << e << ":" << j << " dgv gfc " << f_val << " "
1149 << dgv_gfc_val << " "
1150 << fabs(f_val - dgv_gfc_val)
1151 << std::endl;
1152 }
1153 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
1154 {
1155 std::cout << e << ":" << j << " dgi gfc " << f_val << " "
1156 << dgi_gfc_val << " "
1157 << fabs(f_val - dgi_gfc_val)
1158 << std::endl;
1159 }
1160 if (log > 0 && fabs(f_val - h1_gv_val) > tol)
1161 {
1162 std::cout << e << ":" << j << " h1 gv " << f_val << " "
1163 << h1_gv_val << " " << fabs(f_val - h1_gv_val)
1164 << std::endl;
1165 }
1166 if (log > 0 && fabs(f_val - dgv_gv_val) > tol)
1167 {
1168 std::cout << e << ":" << j << " dgv gv " << f_val << " "
1169 << dgv_gv_val << " "
1170 << fabs(f_val - dgv_gv_val)
1171 << std::endl;
1172 }
1173 if (log > 0 && fabs(f_val - dgi_gv_val) > tol)
1174 {
1175 std::cout << e << ":" << j << " dgi gv " << f_val << " "
1176 << dgi_gv_val << " "
1177 << fabs(f_val - dgi_gv_val)
1178 << std::endl;
1179 }
1180 }
1181
1182 h1_gfc_err /= ir.GetNPoints();
1183 dgv_gfc_err /= ir.GetNPoints();
1184 dgi_gfc_err /= ir.GetNPoints();
1185
1186 h1_gv_err /= ir.GetNPoints();
1187 dgv_gv_err /= ir.GetNPoints();
1188 dgi_gv_err /= ir.GetNPoints();
1189
1190 REQUIRE(h1_gfc_err == MFEM_Approx(0.0));
1191 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
1192 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
1193
1194 REQUIRE(h1_gv_err == MFEM_Approx(0.0));
1195 REQUIRE(dgv_gv_err == MFEM_Approx(0.0));
1196 REQUIRE(dgi_gv_err == MFEM_Approx(0.0));
1197 }
1198 }
1199
1200 SECTION("Boundary Evaluation 3D (H1 Context)")
1201 {
1202 std::cout << "Boundary Evaluation 3D (H1 Context)" << std::endl;
1203 for (int be = 0; be < mesh.GetNBE(); be++)
1204 {
1205 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
1206 const FiniteElement *fe = h1_fespace.GetBE(be);
1207 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
1208 2*order + 2);
1209
1210 double h1_err = 0.0;
1211 double dgv_err = 0.0;
1212 double dgi_err = 0.0;
1213
1214 for (int j=0; j<ir.GetNPoints(); j++)
1215 {
1216 npts++;
1217 const IntegrationPoint &ip = ir.IntPoint(j);
1218 T->SetIntPoint(&ip);
1219
1220 double f_val = funcCoef.Eval(*T, ip);
1221 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
1222 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
1223 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
1224
1225 h1_err += fabs(f_val - h1_gfc_val);
1226 dgv_err += fabs(f_val - dgv_gfc_val);
1227 dgi_err += fabs(f_val - dgi_gfc_val);
1228
1229 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
1230 {
1231 std::cout << be << ":" << j << " h1 " << f_val << " "
1232 << h1_gfc_val << " "
1233 << fabs(f_val - h1_gfc_val)
1234 << std::endl;
1235 }
1236 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
1237 {
1238 std::cout << be << ":" << j << " dgv " << f_val << " "
1239 << dgv_gfc_val << " "
1240 << fabs(f_val - dgv_gfc_val)
1241 << std::endl;
1242 }
1243 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
1244 {
1245 std::cout << be << ":" << j << " dgi " << f_val << " "
1246 << dgi_gfc_val << " "
1247 << fabs(f_val - dgi_gfc_val)
1248 << std::endl;
1249 }
1250 }
1251 h1_err /= ir.GetNPoints();
1252 dgv_err /= ir.GetNPoints();
1253 dgi_err /= ir.GetNPoints();
1254
1255 REQUIRE(h1_err == MFEM_Approx(0.0));
1256 REQUIRE(dgv_err == MFEM_Approx(0.0));
1257 REQUIRE(dgi_err == MFEM_Approx(0.0));
1258 }
1259 }
1260
1261 SECTION("Boundary Evaluation 3D (DG Context)")
1262 {
1263 std::cout << "Boundary Evaluation 3D (DG Context)" << std::endl;
1264 for (int be = 0; be < mesh.GetNBE(); be++)
1265 {
1266 FaceElementTransformations *T =
1267 mesh.GetBdrFaceTransformations(be);
1268 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
1269 2*order + 2);
1270
1271 double h1_err = 0.0;
1272 double dgv_err = 0.0;
1273 double dgi_err = 0.0;
1274
1275 for (int j=0; j<ir.GetNPoints(); j++)
1276 {
1277 npts++;
1278 const IntegrationPoint &ip = ir.IntPoint(j);
1279
1280 T->SetIntPoint(&ip);
1281
1282 double f_val = funcCoef.Eval(*T, ip);
1283 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
1284 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
1285 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
1286
1287 h1_err += fabs(f_val - h1_gfc_val);
1288 dgv_err += fabs(f_val - dgv_gfc_val);
1289 dgi_err += fabs(f_val - dgi_gfc_val);
1290
1291 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
1292 {
1293 std::cout << be << ":" << j << " h1 " << f_val << " "
1294 << h1_gfc_val << " "
1295 << fabs(f_val - h1_gfc_val)
1296 << std::endl;
1297 }
1298 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
1299 {
1300 std::cout << be << ":" << j << " dgv " << f_val << " "
1301 << dgv_gfc_val << " "
1302 << fabs(f_val - dgv_gfc_val)
1303 << std::endl;
1304 }
1305 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
1306 {
1307 std::cout << be << ":" << j << " dgi " << f_val << " "
1308 << dgi_gfc_val << " "
1309 << fabs(f_val - dgi_gfc_val)
1310 << std::endl;
1311 }
1312 }
1313 h1_err /= ir.GetNPoints();
1314 dgv_err /= ir.GetNPoints();
1315 dgi_err /= ir.GetNPoints();
1316
1317 REQUIRE(h1_err == MFEM_Approx(0.0));
1318 REQUIRE(dgv_err == MFEM_Approx(0.0));
1319 REQUIRE(dgi_err == MFEM_Approx(0.0));
1320 }
1321 }
1322
1323 SECTION("Edge Evaluation 3D (H1 Context)")
1324 {
1325 std::cout << "Edge Evaluation 3D (H1 Context)" << std::endl;
1326 for (int e = 0; e < mesh.GetNEdges(); e++)
1327 {
1328 ElementTransformation *T = mesh.GetEdgeTransformation(e);
1329 const FiniteElement *fe = h1_fespace.GetEdgeElement(e);
1330 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
1331 2*order + 2);
1332
1333 double h1_err = 0.0;
1334
1335 for (int j=0; j<ir.GetNPoints(); j++)
1336 {
1337 npts++;
1338 const IntegrationPoint &ip = ir.IntPoint(j);
1339 T->SetIntPoint(&ip);
1340
1341 double f_val = funcCoef.Eval(*T, ip);
1342 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
1343
1344 h1_err += fabs(f_val - h1_gfc_val);
1345
1346 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
1347 {
1348 std::cout << e << ":" << j << " h1 " << f_val << " "
1349 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
1350 << std::endl;
1351 }
1352 }
1353 h1_err /= ir.GetNPoints();
1354
1355 REQUIRE(h1_err == MFEM_Approx(0.0));
1356 }
1357 }
1358
1359 SECTION("Face Evaluation 3D (H1 Context)")
1360 {
1361 std::cout << "Face Evaluation 3D (H1 Context)" << std::endl;
1362 for (int f = 0; f < mesh.GetNFaces(); f++)
1363 {
1364 ElementTransformation *T = mesh.GetFaceTransformation(f);
1365 const FiniteElement *fe = h1_fespace.GetFaceElement(f);
1366 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
1367 2*order + 2);
1368
1369 double h1_err = 0.0;
1370
1371 for (int j=0; j<ir.GetNPoints(); j++)
1372 {
1373 npts++;
1374 const IntegrationPoint &ip = ir.IntPoint(j);
1375 T->SetIntPoint(&ip);
1376
1377 double f_val = funcCoef.Eval(*T, ip);
1378 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
1379
1380 h1_err += fabs(f_val - h1_gfc_val);
1381
1382 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
1383 {
1384 std::cout << f << ":" << j << " h1 " << f_val << " "
1385 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
1386 << std::endl;
1387 }
1388 }
1389 h1_err /= ir.GetNPoints();
1390
1391 REQUIRE(h1_err == MFEM_Approx(0.0));
1392 }
1393 }
1394 }
1395 }
1396 std::cout << "Checked GridFunction::GetValue at "
1397 << npts << " 3D points" << std::endl;
1398 }
1399
1400 #ifdef MFEM_USE_MPI
1401 #
1402 TEST_CASE("3D GetValue in Parallel",
1403 "[ParGridFunction]"
1404 "[GridFunctionCoefficient]"
1405 "[Parallel]")
1406 {
1407 int num_procs;
1408 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
1409
1410 int my_rank;
1411 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
1412
1413 int log = 1;
1414 int n = (int)ceil(pow(2*num_procs, 1.0 / 3.0));
1415 int dim = 3;
1416 int order = 1;
1417 int npts = 0;
1418
1419 double tol = 1e-6;
1420
1421 for (int type = (int)Element::TETRAHEDRON;
1422 type <= (int)Element::WEDGE; type++)
1423 {
1424 Mesh mesh = Mesh::MakeCartesian3D(
1425 n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
1426 ParMesh pmesh(MPI_COMM_WORLD, mesh);
1427 mesh.Clear();
1428
1429 FunctionCoefficient funcCoef(func_3D_lin);
1430
to_string(type)1431 SECTION("3D GetValue tests for element type " + std::to_string(type))
1432 {
1433 H1_FECollection h1_fec(order, dim);
1434 L2_FECollection l2_fec(order, dim);
1435 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
1436 FiniteElement::VALUE);
1437 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
1438 FiniteElement::INTEGRAL);
1439
1440 ParFiniteElementSpace h1_fespace(&pmesh, &h1_fec);
1441 ParFiniteElementSpace l2_fespace(&pmesh, &l2_fec);
1442 ParFiniteElementSpace dgv_fespace(&pmesh, &dgv_fec);
1443 ParFiniteElementSpace dgi_fespace(&pmesh, &dgi_fec);
1444
1445 ParGridFunction h1_x( &h1_fespace);
1446 ParGridFunction l2_x( &l2_fespace);
1447 ParGridFunction dgv_x(&dgv_fespace);
1448 ParGridFunction dgi_x(&dgi_fespace);
1449
1450 GridFunctionCoefficient h1_xCoef( &h1_x);
1451 GridFunctionCoefficient l2_xCoef( &l2_x);
1452 GridFunctionCoefficient dgv_xCoef(&dgv_x);
1453 GridFunctionCoefficient dgi_xCoef(&dgi_x);
1454
1455 h1_x.ProjectCoefficient(funcCoef);
1456 l2_x.ProjectCoefficient(funcCoef);
1457 dgv_x.ProjectCoefficient(funcCoef);
1458 dgi_x.ProjectCoefficient(funcCoef);
1459
1460 h1_x.ExchangeFaceNbrData();
1461 l2_x.ExchangeFaceNbrData();
1462 dgv_x.ExchangeFaceNbrData();
1463 dgi_x.ExchangeFaceNbrData();
1464
1465 SECTION("Shared Face Evaluation 3D")
1466 {
1467 if (my_rank == 0)
1468 {
1469 std::cout << "Shared Face Evaluation 3D" << std::endl;
1470 }
1471 for (int sf = 0; sf < pmesh.GetNSharedFaces(); sf++)
1472 {
1473 FaceElementTransformations *FET =
1474 pmesh.GetSharedFaceTransformations(sf);
1475 ElementTransformation *T = &FET->GetElement2Transformation();
1476 int e = FET->Elem2No;
1477 int e_nbr = e - pmesh.GetNE();
1478 const FiniteElement *fe = dgv_fespace.GetFaceNbrFE(e_nbr);
1479 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
1480 2*order + 2);
1481
1482 double h1_gfc_err = 0.0;
1483 double dgv_gfc_err = 0.0;
1484 double dgi_gfc_err = 0.0;
1485
1486 double h1_gv_err = 0.0;
1487 double dgv_gv_err = 0.0;
1488 double dgi_gv_err = 0.0;
1489
1490 for (int j=0; j<ir.GetNPoints(); j++)
1491 {
1492 npts++;
1493 const IntegrationPoint &ip = ir.IntPoint(j);
1494 T->SetIntPoint(&ip);
1495
1496 double f_val = funcCoef.Eval(*T, ip);
1497
1498 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
1499 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
1500 double dgi_gfc_val = dgi_xCoef.Eval(*T, ip);
1501
1502 double h1_gv_val = h1_x.GetValue(e, ip);
1503 double dgv_gv_val = dgv_x.GetValue(e, ip);
1504 double dgi_gv_val = dgi_x.GetValue(e, ip);
1505
1506 h1_gfc_err += fabs(f_val - h1_gfc_val);
1507 dgv_gfc_err += fabs(f_val - dgv_gfc_val);
1508 dgi_gfc_err += fabs(f_val - dgi_gfc_val);
1509
1510 h1_gv_err += fabs(f_val - h1_gv_val);
1511 dgv_gv_err += fabs(f_val - dgv_gv_val);
1512 dgi_gv_err += fabs(f_val - dgi_gv_val);
1513
1514 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
1515 {
1516 std::cout << e << ":" << j << " h1 gfc " << f_val << " "
1517 << h1_gfc_val << " " << fabs(f_val - h1_gfc_val)
1518 << std::endl;
1519 }
1520 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
1521 {
1522 std::cout << e << ":" << j << " dgv gfc " << f_val << " "
1523 << dgv_gfc_val << " "
1524 << fabs(f_val - dgv_gfc_val)
1525 << std::endl;
1526 }
1527 if (log > 0 && fabs(f_val - dgi_gfc_val) > tol)
1528 {
1529 std::cout << e << ":" << j << " dgi gfc " << f_val << " "
1530 << dgi_gfc_val << " "
1531 << fabs(f_val - dgi_gfc_val)
1532 << std::endl;
1533 }
1534 if (log > 0 && fabs(f_val - h1_gv_val) > tol)
1535 {
1536 std::cout << e << ":" << j << " h1 gv " << f_val << " "
1537 << h1_gv_val << " " << fabs(f_val - h1_gv_val)
1538 << std::endl;
1539 }
1540 if (log > 0 && fabs(f_val - dgv_gv_val) > tol)
1541 {
1542 std::cout << e << ":" << j << " dgv gv " << f_val << " "
1543 << dgv_gv_val << " "
1544 << fabs(f_val - dgv_gv_val)
1545 << std::endl;
1546 }
1547 if (log > 0 && fabs(f_val - dgi_gv_val) > tol)
1548 {
1549 std::cout << e << ":" << j << " dgi gv " << f_val << " "
1550 << dgi_gv_val << " "
1551 << fabs(f_val - dgi_gv_val)
1552 << std::endl;
1553 }
1554 }
1555
1556 h1_gfc_err /= ir.GetNPoints();
1557 dgv_gfc_err /= ir.GetNPoints();
1558 dgi_gfc_err /= ir.GetNPoints();
1559
1560 h1_gv_err /= ir.GetNPoints();
1561 dgv_gv_err /= ir.GetNPoints();
1562 dgi_gv_err /= ir.GetNPoints();
1563
1564 REQUIRE(h1_gfc_err == MFEM_Approx(0.0));
1565 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
1566 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
1567
1568 REQUIRE(h1_gv_err == MFEM_Approx(0.0));
1569 REQUIRE(dgv_gv_err == MFEM_Approx(0.0));
1570 REQUIRE(dgi_gv_err == MFEM_Approx(0.0));
1571 }
1572 }
1573 }
1574 }
1575 std::cout << my_rank << ": Checked GridFunction::GetValue at "
1576 << npts << " 3D points" << std::endl;
1577 }
1578
1579 #endif // MFEM_USE_MPI
1580
1581 TEST_CASE("2D GetVectorValue",
1582 "[GridFunction]"
1583 "[VectorGridFunctionCoefficient]")
1584 {
1585 int log = 1;
1586 int n = 1;
1587 int dim = 2;
1588 int order = 1;
1589 int npts = 0;
1590
1591 double tol = 1e-6;
1592
1593 for (int type = (int)Element::TRIANGLE;
1594 type <= (int)Element::QUADRILATERAL; type++)
1595 {
1596 Mesh mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
1597
1598 VectorFunctionCoefficient funcCoef(dim, Func_2D_lin);
1599
1600 SECTION("2D GetVectorValue tests for element type " +
to_string(type)1601 std::to_string(type))
1602 {
1603 H1_FECollection h1_fec(order, dim);
1604 ND_FECollection nd_fec(order+1, dim);
1605 RT_FECollection rt_fec(order+1, dim);
1606 L2_FECollection l2_fec(order, dim);
1607 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
1608 FiniteElement::VALUE);
1609 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
1610 FiniteElement::INTEGRAL);
1611
1612 FiniteElementSpace h1_fespace(&mesh, &h1_fec, dim);
1613 FiniteElementSpace nd_fespace(&mesh, &nd_fec);
1614 FiniteElementSpace rt_fespace(&mesh, &rt_fec);
1615 FiniteElementSpace l2_fespace(&mesh, &l2_fec, dim);
1616 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec, dim);
1617 FiniteElementSpace dgi_fespace(&mesh, &dgi_fec, dim);
1618
1619 GridFunction h1_x( &h1_fespace);
1620 GridFunction nd_x( &nd_fespace);
1621 GridFunction rt_x( &rt_fespace);
1622 GridFunction l2_x( &l2_fespace);
1623 GridFunction dgv_x(&dgv_fespace);
1624 GridFunction dgi_x(&dgi_fespace);
1625
1626 VectorGridFunctionCoefficient h1_xCoef( &h1_x);
1627 VectorGridFunctionCoefficient nd_xCoef( &nd_x);
1628 VectorGridFunctionCoefficient rt_xCoef( &rt_x);
1629 VectorGridFunctionCoefficient l2_xCoef( &l2_x);
1630 VectorGridFunctionCoefficient dgv_xCoef(&dgv_x);
1631 VectorGridFunctionCoefficient dgi_xCoef(&dgi_x);
1632
1633 h1_x.ProjectCoefficient(funcCoef);
1634 nd_x.ProjectCoefficient(funcCoef);
1635 rt_x.ProjectCoefficient(funcCoef);
1636 l2_x.ProjectCoefficient(funcCoef);
1637 dgv_x.ProjectCoefficient(funcCoef);
1638 dgi_x.ProjectCoefficient(funcCoef);
1639
1640 Vector f_val(dim); f_val = 0.0;
1641
1642 Vector h1_gfc_val(dim); h1_gfc_val = 0.0;
1643 Vector nd_gfc_val(dim); nd_gfc_val = 0.0;
1644 Vector rt_gfc_val(dim); rt_gfc_val = 0.0;
1645 Vector l2_gfc_val(dim); l2_gfc_val = 0.0;
1646 Vector dgv_gfc_val(dim); dgv_gfc_val = 0.0;
1647 Vector dgi_gfc_val(dim); dgi_gfc_val = 0.0;
1648
1649 Vector h1_gvv_val(dim); h1_gvv_val = 0.0;
1650 Vector nd_gvv_val(dim); nd_gvv_val = 0.0;
1651 Vector rt_gvv_val(dim); rt_gvv_val = 0.0;
1652 Vector l2_gvv_val(dim); l2_gvv_val = 0.0;
1653 Vector dgv_gvv_val(dim); dgv_gvv_val = 0.0;
1654 Vector dgi_gvv_val(dim); dgi_gvv_val = 0.0;
1655
1656 SECTION("Domain Evaluation 2D")
1657 {
1658 std::cout << "Domain Evaluation 2D" << std::endl;
1659 for (int e = 0; e < mesh.GetNE(); e++)
1660 {
1661 ElementTransformation *T = mesh.GetElementTransformation(e);
1662 const FiniteElement *fe = h1_fespace.GetFE(e);
1663 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
1664 2*order + 2);
1665
1666 double h1_gfc_err = 0.0;
1667 double nd_gfc_err = 0.0;
1668 double rt_gfc_err = 0.0;
1669 double l2_gfc_err = 0.0;
1670 double dgv_gfc_err = 0.0;
1671 double dgi_gfc_err = 0.0;
1672
1673 double h1_gvv_err = 0.0;
1674 double nd_gvv_err = 0.0;
1675 double rt_gvv_err = 0.0;
1676 double l2_gvv_err = 0.0;
1677 double dgv_gvv_err = 0.0;
1678 double dgi_gvv_err = 0.0;
1679
1680 for (int j=0; j<ir.GetNPoints(); j++)
1681 {
1682 npts++;
1683 const IntegrationPoint &ip = ir.IntPoint(j);
1684 T->SetIntPoint(&ip);
1685
1686 funcCoef.Eval(f_val, *T, ip);
1687
1688 h1_xCoef.Eval(h1_gfc_val, *T, ip);
1689 nd_xCoef.Eval(nd_gfc_val, *T, ip);
1690 rt_xCoef.Eval(rt_gfc_val, *T, ip);
1691 l2_xCoef.Eval(l2_gfc_val, *T, ip);
1692 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
1693 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
1694
1695 h1_x.GetVectorValue(e, ip, h1_gvv_val);
1696 nd_x.GetVectorValue(e, ip, nd_gvv_val);
1697 rt_x.GetVectorValue(e, ip, rt_gvv_val);
1698 l2_x.GetVectorValue(e, ip, l2_gvv_val);
1699 dgv_x.GetVectorValue(e, ip, dgv_gvv_val);
1700 dgi_x.GetVectorValue(e, ip, dgi_gvv_val);
1701
1702 double h1_gfc_dist = Distance(f_val, h1_gfc_val, dim);
1703 double nd_gfc_dist = Distance(f_val, nd_gfc_val, dim);
1704 double rt_gfc_dist = Distance(f_val, rt_gfc_val, dim);
1705 double l2_gfc_dist = Distance(f_val, l2_gfc_val, dim);
1706 double dgv_gfc_dist = Distance(f_val, dgv_gfc_val, dim);
1707 double dgi_gfc_dist = Distance(f_val, dgi_gfc_val, dim);
1708
1709 double h1_gvv_dist = Distance(f_val, h1_gvv_val, dim);
1710 double nd_gvv_dist = Distance(f_val, nd_gvv_val, dim);
1711 double rt_gvv_dist = Distance(f_val, rt_gvv_val, dim);
1712 double l2_gvv_dist = Distance(f_val, l2_gvv_val, dim);
1713 double dgv_gvv_dist = Distance(f_val, dgv_gvv_val, dim);
1714 double dgi_gvv_dist = Distance(f_val, dgi_gvv_val, dim);
1715
1716 h1_gfc_err += h1_gfc_dist;
1717 nd_gfc_err += nd_gfc_dist;
1718 rt_gfc_err += rt_gfc_dist;
1719 l2_gfc_err += l2_gfc_dist;
1720 dgv_gfc_err += dgv_gfc_dist;
1721 dgi_gfc_err += dgi_gfc_dist;
1722
1723 h1_gvv_err += h1_gvv_dist;
1724 nd_gvv_err += nd_gvv_dist;
1725 rt_gvv_err += rt_gvv_dist;
1726 l2_gvv_err += l2_gvv_dist;
1727 dgv_gvv_err += dgv_gvv_dist;
1728 dgi_gvv_err += dgi_gvv_dist;
1729
1730 if (log > 0 && h1_gfc_dist > tol)
1731 {
1732 std::cout << e << ":" << j << " h1 gfc ("
1733 << f_val[0] << "," << f_val[1] << ") vs. ("
1734 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ") "
1735 << h1_gfc_dist << std::endl;
1736 }
1737 if (log > 0 && nd_gfc_dist > tol)
1738 {
1739 std::cout << e << ":" << j << " nd gfc ("
1740 << f_val[0] << "," << f_val[1] << ") vs. ("
1741 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ") "
1742 << nd_gfc_dist << std::endl;
1743 }
1744 if (log > 0 && rt_gfc_dist > tol)
1745 {
1746 std::cout << e << ":" << j << " rt gfc ("
1747 << f_val[0] << "," << f_val[1] << ") vs. ("
1748 << rt_gfc_val[0] << "," << rt_gfc_val[1] << ") "
1749 << rt_gfc_dist << std::endl;
1750 }
1751 if (log > 0 && l2_gfc_dist > tol)
1752 {
1753 std::cout << e << ":" << j << " l2 gfc ("
1754 << f_val[0] << "," << f_val[1] << ") vs. ("
1755 << l2_gfc_val[0] << "," << l2_gfc_val[1] << ") "
1756 << l2_gfc_dist << std::endl;
1757 }
1758 if (log > 0 && dgv_gfc_dist > tol)
1759 {
1760 std::cout << e << ":" << j << " dgv gfc ("
1761 << f_val[0] << "," << f_val[1] << ") vs. ("
1762 << dgv_gfc_val[0] << ","
1763 << dgv_gfc_val[1] << ") "
1764 << dgv_gfc_dist << std::endl;
1765 }
1766 if (log > 0 && dgi_gfc_dist > tol)
1767 {
1768 std::cout << e << ":" << j << " dgi gfc ("
1769 << f_val[0] << "," << f_val[1] << ") vs. ("
1770 << dgi_gfc_val[0] << ","
1771 << dgi_gfc_val[1] << ") "
1772 << dgi_gfc_dist << std::endl;
1773 }
1774 if (log > 0 && h1_gvv_dist > tol)
1775 {
1776 std::cout << e << ":" << j << " h1 gvv ("
1777 << f_val[0] << "," << f_val[1] << ") vs. ("
1778 << h1_gvv_val[0] << "," << h1_gvv_val[1] << ") "
1779 << h1_gvv_dist << std::endl;
1780 }
1781 if (log > 0 && nd_gvv_dist > tol)
1782 {
1783 std::cout << e << ":" << j << " nd gvv ("
1784 << f_val[0] << "," << f_val[1] << ") vs. ("
1785 << nd_gvv_val[0] << "," << nd_gvv_val[1] << ") "
1786 << nd_gvv_dist << std::endl;
1787 }
1788 if (log > 0 && rt_gvv_dist > tol)
1789 {
1790 std::cout << e << ":" << j << " rt gvv ("
1791 << f_val[0] << "," << f_val[1] << ") vs. ("
1792 << rt_gvv_val[0] << "," << rt_gvv_val[1] << ") "
1793 << rt_gvv_dist << std::endl;
1794 }
1795 if (log > 0 && l2_gvv_dist > tol)
1796 {
1797 std::cout << e << ":" << j << " l2 gvv ("
1798 << f_val[0] << "," << f_val[1] << ") vs. ("
1799 << l2_gvv_val[0] << "," << l2_gvv_val[1] << ") "
1800 << l2_gvv_dist << std::endl;
1801 }
1802 if (log > 0 && dgv_gvv_dist > tol)
1803 {
1804 std::cout << e << ":" << j << " dgv gvv ("
1805 << f_val[0] << "," << f_val[1] << ") vs. ("
1806 << dgv_gvv_val[0] << ","
1807 << dgv_gvv_val[1] << ") "
1808 << dgv_gvv_dist << std::endl;
1809 }
1810 if (log > 0 && dgi_gvv_dist > tol)
1811 {
1812 std::cout << e << ":" << j << " dgi gvv ("
1813 << f_val[0] << "," << f_val[1] << ") vs. ("
1814 << dgi_gvv_val[0] << ","
1815 << dgi_gvv_val[1] << ") "
1816 << dgi_gvv_dist << std::endl;
1817 }
1818 }
1819
1820 h1_gfc_err /= ir.GetNPoints();
1821 nd_gfc_err /= ir.GetNPoints();
1822 rt_gfc_err /= ir.GetNPoints();
1823 l2_gfc_err /= ir.GetNPoints();
1824 dgv_gfc_err /= ir.GetNPoints();
1825 dgi_gfc_err /= ir.GetNPoints();
1826
1827 h1_gvv_err /= ir.GetNPoints();
1828 nd_gvv_err /= ir.GetNPoints();
1829 rt_gvv_err /= ir.GetNPoints();
1830 l2_gvv_err /= ir.GetNPoints();
1831 dgv_gvv_err /= ir.GetNPoints();
1832 dgi_gvv_err /= ir.GetNPoints();
1833
1834 REQUIRE( h1_gfc_err == MFEM_Approx(0.0));
1835 REQUIRE( nd_gfc_err == MFEM_Approx(0.0));
1836 REQUIRE( rt_gfc_err == MFEM_Approx(0.0));
1837 REQUIRE( l2_gfc_err == MFEM_Approx(0.0));
1838 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
1839 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
1840
1841 REQUIRE( h1_gvv_err == MFEM_Approx(0.0));
1842 REQUIRE( nd_gvv_err == MFEM_Approx(0.0));
1843 REQUIRE( rt_gvv_err == MFEM_Approx(0.0));
1844 REQUIRE( l2_gvv_err == MFEM_Approx(0.0));
1845 REQUIRE(dgv_gvv_err == MFEM_Approx(0.0));
1846 REQUIRE(dgi_gvv_err == MFEM_Approx(0.0));
1847 }
1848 }
1849
1850 SECTION("Boundary Evaluation 2D (H1 Context)")
1851 {
1852 std::cout << "Boundary Evaluation 2D (H1 Context)" << std::endl;
1853 for (int be = 0; be < mesh.GetNBE(); be++)
1854 {
1855 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
1856 const FiniteElement *fe = h1_fespace.GetBE(be);
1857 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
1858 2*order + 2);
1859
1860 double h1_err = 0.0;
1861 double nd_err = 0.0;
1862 double rt_err = 0.0;
1863 double l2_err = 0.0;
1864 double dgv_err = 0.0;
1865 double dgi_err = 0.0;
1866
1867 for (int j=0; j<ir.GetNPoints(); j++)
1868 {
1869 npts++;
1870 const IntegrationPoint &ip = ir.IntPoint(j);
1871 T->SetIntPoint(&ip);
1872
1873 funcCoef.Eval(f_val, *T, ip);
1874 h1_xCoef.Eval(h1_gfc_val, *T, ip);
1875 nd_xCoef.Eval(nd_gfc_val, *T, ip);
1876 rt_xCoef.Eval(rt_gfc_val, *T, ip);
1877 l2_xCoef.Eval(l2_gfc_val, *T, ip);
1878 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
1879 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
1880
1881 double h1_dist = Distance(f_val, h1_gfc_val, dim);
1882 double nd_dist = Distance(f_val, nd_gfc_val, dim);
1883 double rt_dist = Distance(f_val, rt_gfc_val, dim);
1884 double l2_dist = Distance(f_val, l2_gfc_val, dim);
1885 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
1886 double dgi_dist = Distance(f_val, dgi_gfc_val, dim);
1887
1888 h1_err += h1_dist;
1889 nd_err += nd_dist;
1890 rt_err += rt_dist;
1891 l2_err += l2_dist;
1892 dgv_err += dgv_dist;
1893 dgi_err += dgi_dist;
1894
1895 if (log > 0 && h1_dist > tol)
1896 {
1897 std::cout << be << ":" << j << " h1 ("
1898 << f_val[0] << "," << f_val[1] << ") vs. ("
1899 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ") "
1900 << h1_dist << std::endl;
1901 }
1902 if (log > 0 && nd_dist > tol)
1903 {
1904 std::cout << be << ":" << j << " nd ("
1905 << f_val[0] << "," << f_val[1] << ") vs. ("
1906 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ") "
1907 << nd_dist << std::endl;
1908 }
1909 if (log > 0 && rt_dist > tol)
1910 {
1911 std::cout << be << ":" << j << " rt ("
1912 << f_val[0] << "," << f_val[1] << ") vs. ("
1913 << rt_gfc_val[0] << "," << rt_gfc_val[1] << ") "
1914 << rt_dist << std::endl;
1915 }
1916 if (log > 0 && l2_dist > tol)
1917 {
1918 std::cout << be << ":" << j << " l2 ("
1919 << f_val[0] << "," << f_val[1] << ") vs. ("
1920 << l2_gfc_val[0] << "," << l2_gfc_val[1] << ") "
1921 << l2_dist << std::endl;
1922 }
1923 if (log > 0 && dgv_dist > tol)
1924 {
1925 std::cout << be << ":" << j << " dgv ("
1926 << f_val[0] << "," << f_val[1] << ") vs. ("
1927 << dgv_gfc_val[0] << ","
1928 << dgv_gfc_val[1] << ") "
1929 << dgv_dist << std::endl;
1930 }
1931 if (log > 0 && dgi_dist > tol)
1932 {
1933 std::cout << be << ":" << j << " dgi ("
1934 << f_val[0] << "," << f_val[1] << ") vs. ("
1935 << dgi_gfc_val[0] << ","
1936 << dgi_gfc_val[1] << ") "
1937 << dgi_dist << std::endl;
1938 }
1939 }
1940 h1_err /= ir.GetNPoints();
1941 nd_err /= ir.GetNPoints();
1942 rt_err /= ir.GetNPoints();
1943 l2_err /= ir.GetNPoints();
1944 dgv_err /= ir.GetNPoints();
1945 dgi_err /= ir.GetNPoints();
1946
1947 REQUIRE( h1_err == MFEM_Approx(0.0));
1948 REQUIRE( nd_err == MFEM_Approx(0.0));
1949 REQUIRE( rt_err == MFEM_Approx(0.0));
1950 REQUIRE( l2_err == MFEM_Approx(0.0));
1951 REQUIRE(dgv_err == MFEM_Approx(0.0));
1952 REQUIRE(dgi_err == MFEM_Approx(0.0));
1953 }
1954 }
1955
1956 SECTION("Boundary Evaluation 2D (DG Context)")
1957 {
1958 std::cout << "Boundary Evaluation 2D (DG Context)" << std::endl;
1959 for (int be = 0; be < mesh.GetNBE(); be++)
1960 {
1961 FaceElementTransformations *T =
1962 mesh.GetBdrFaceTransformations(be);
1963 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
1964 2*order + 2);
1965
1966 double h1_err = 0.0;
1967 double nd_err = 0.0;
1968 double rt_err = 0.0;
1969 double l2_err = 0.0;
1970 double dgv_err = 0.0;
1971 double dgi_err = 0.0;
1972
1973 for (int j=0; j<ir.GetNPoints(); j++)
1974 {
1975 npts++;
1976 const IntegrationPoint &ip = ir.IntPoint(j);
1977
1978 T->SetIntPoint(&ip);
1979
1980 funcCoef.Eval(f_val, *T, ip);
1981 h1_xCoef.Eval(h1_gfc_val, *T, ip);
1982 nd_xCoef.Eval(nd_gfc_val, *T, ip);
1983 rt_xCoef.Eval(rt_gfc_val, *T, ip);
1984 l2_xCoef.Eval(l2_gfc_val, *T, ip);
1985 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
1986 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
1987
1988 double h1_dist = Distance(f_val, h1_gfc_val, dim);
1989 double nd_dist = Distance(f_val, nd_gfc_val, dim);
1990 double rt_dist = Distance(f_val, rt_gfc_val, dim);
1991 double l2_dist = Distance(f_val, l2_gfc_val, dim);
1992 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
1993 double dgi_dist = Distance(f_val, dgi_gfc_val, dim);
1994
1995 h1_err += h1_dist;
1996 nd_err += nd_dist;
1997 rt_err += rt_dist;
1998 l2_err += l2_dist;
1999 dgv_err += dgv_dist;
2000 dgi_err += dgi_dist;
2001
2002 if (log > 0 && h1_dist > tol)
2003 {
2004 std::cout << be << ":" << j << " h1 ("
2005 << f_val[0] << "," << f_val[1] << ") vs. ("
2006 << h1_gfc_val[0] << ","
2007 << h1_gfc_val[1] << ") "
2008 << h1_dist << std::endl;
2009 }
2010 if (log > 0 && nd_dist > tol)
2011 {
2012 std::cout << be << ":" << j << " nd ("
2013 << f_val[0] << "," << f_val[1] << ") vs. ("
2014 << nd_gfc_val[0] << ","
2015 << nd_gfc_val[1] << ") "
2016 << nd_dist << std::endl;
2017 }
2018 if (log > 0 && rt_dist > tol)
2019 {
2020 std::cout << be << ":" << j << " rt ("
2021 << f_val[0] << "," << f_val[1] << ") vs. ("
2022 << rt_gfc_val[0] << ","
2023 << rt_gfc_val[1] << ") "
2024 << rt_dist << std::endl;
2025 }
2026 if (log > 0 && l2_dist > tol)
2027 {
2028 std::cout << be << ":" << j << " l2 ("
2029 << f_val[0] << "," << f_val[1] << ") vs. ("
2030 << l2_gfc_val[0] << ","
2031 << l2_gfc_val[1] << ") "
2032 << l2_dist << std::endl;
2033 }
2034 if (log > 0 && dgv_dist > tol)
2035 {
2036 std::cout << be << ":" << j << " dgv ("
2037 << f_val[0] << "," << f_val[1] << ") vs. ("
2038 << dgv_gfc_val[0] << ","
2039 << dgv_gfc_val[1] << ") "
2040 << dgv_dist << std::endl;
2041 }
2042 if (log > 0 && dgi_dist > tol)
2043 {
2044 std::cout << be << ":" << j << " dgi ("
2045 << f_val[0] << "," << f_val[1] << ") vs. ("
2046 << dgi_gfc_val[0] << ","
2047 << dgi_gfc_val[1] << ") "
2048 << dgi_dist << std::endl;
2049 }
2050 }
2051 h1_err /= ir.GetNPoints();
2052 nd_err /= ir.GetNPoints();
2053 rt_err /= ir.GetNPoints();
2054 l2_err /= ir.GetNPoints();
2055 dgv_err /= ir.GetNPoints();
2056 dgi_err /= ir.GetNPoints();
2057
2058 REQUIRE( h1_err == MFEM_Approx(0.0));
2059 REQUIRE( nd_err == MFEM_Approx(0.0));
2060 REQUIRE( rt_err == MFEM_Approx(0.0));
2061 REQUIRE( l2_err == MFEM_Approx(0.0));
2062 REQUIRE(dgv_err == MFEM_Approx(0.0));
2063 REQUIRE(dgi_err == MFEM_Approx(0.0));
2064 }
2065 }
2066
2067 SECTION("Edge Evaluation 2D")
2068 {
2069 std::cout << "Edge Evaluation 2D" << std::endl;
2070 for (int e = 0; e < mesh.GetNEdges(); e++)
2071 {
2072 ElementTransformation *T = mesh.GetEdgeTransformation(e);
2073 const FiniteElement *fe = h1_fespace.GetEdgeElement(e);
2074 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
2075 2*order + 2);
2076
2077 double h1_err = 0.0;
2078
2079 for (int j=0; j<ir.GetNPoints(); j++)
2080 {
2081 npts++;
2082 const IntegrationPoint &ip = ir.IntPoint(j);
2083 T->SetIntPoint(&ip);
2084
2085 funcCoef.Eval(f_val, *T, ip);
2086 h1_xCoef.Eval(h1_gfc_val, *T, ip);
2087
2088 double h1_dist = Distance(f_val, h1_gfc_val, 2);
2089
2090 h1_err += h1_dist;
2091
2092 if (log > 0 && h1_dist > tol)
2093 {
2094 std::cout << e << ":" << j << " h1 ("
2095 << f_val[0] << "," << f_val[1] << ") vs. ("
2096 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ") "
2097 << h1_dist << std::endl;
2098 }
2099 }
2100 h1_err /= ir.GetNPoints();
2101
2102 REQUIRE( h1_err == MFEM_Approx(0.0));
2103 }
2104 }
2105 }
2106 }
2107 std::cout << "Checked GridFunction::GetVectorValue at "
2108 << npts << " 2D points" << std::endl;
2109 }
2110
2111 #ifdef MFEM_USE_MPI
2112 #
2113 TEST_CASE("2D GetVectorValue in Parallel",
2114 "[ParGridFunction]"
2115 "[VectorGridFunctionCoefficient]"
2116 "[Parallel]")
2117 {
2118 int num_procs;
2119 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
2120
2121 int my_rank;
2122 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
2123
2124 int log = 1;
2125 int n = (int)ceil(sqrt(2*num_procs));
2126 int dim = 2;
2127 int order = 1;
2128 int npts = 0;
2129
2130 double tol = 1e-6;
2131
2132 for (int type = (int)Element::TRIANGLE;
2133 type <= (int)Element::QUADRILATERAL; type++)
2134 {
2135 Mesh mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
2136 ParMesh pmesh(MPI_COMM_WORLD, mesh);
2137 mesh.Clear();
2138
2139 VectorFunctionCoefficient funcCoef(dim, Func_2D_lin);
2140
2141 SECTION("2D GetVectorValue tests for element type " +
to_string(type)2142 std::to_string(type))
2143 {
2144 H1_FECollection h1_fec(order, dim);
2145 ND_FECollection nd_fec(order+1, dim);
2146 RT_FECollection rt_fec(order+1, dim);
2147 L2_FECollection l2_fec(order, dim);
2148 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
2149 FiniteElement::VALUE);
2150 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
2151 FiniteElement::INTEGRAL);
2152
2153 ParFiniteElementSpace h1_fespace(&pmesh, &h1_fec, dim);
2154 ParFiniteElementSpace nd_fespace(&pmesh, &nd_fec);
2155 ParFiniteElementSpace rt_fespace(&pmesh, &rt_fec);
2156 ParFiniteElementSpace l2_fespace(&pmesh, &l2_fec, dim);
2157 ParFiniteElementSpace dgv_fespace(&pmesh, &dgv_fec, dim);
2158 ParFiniteElementSpace dgi_fespace(&pmesh, &dgi_fec, dim);
2159
2160 ParGridFunction h1_x( &h1_fespace);
2161 ParGridFunction nd_x( &nd_fespace);
2162 ParGridFunction rt_x( &rt_fespace);
2163 ParGridFunction l2_x( &l2_fespace);
2164 ParGridFunction dgv_x(&dgv_fespace);
2165 ParGridFunction dgi_x(&dgi_fespace);
2166
2167 VectorGridFunctionCoefficient h1_xCoef( &h1_x);
2168 VectorGridFunctionCoefficient nd_xCoef( &nd_x);
2169 VectorGridFunctionCoefficient rt_xCoef( &rt_x);
2170 VectorGridFunctionCoefficient l2_xCoef( &l2_x);
2171 VectorGridFunctionCoefficient dgv_xCoef(&dgv_x);
2172 VectorGridFunctionCoefficient dgi_xCoef(&dgi_x);
2173
2174 h1_x.ProjectCoefficient(funcCoef);
2175 nd_x.ProjectCoefficient(funcCoef);
2176 rt_x.ProjectCoefficient(funcCoef);
2177 l2_x.ProjectCoefficient(funcCoef);
2178 dgv_x.ProjectCoefficient(funcCoef);
2179 dgi_x.ProjectCoefficient(funcCoef);
2180
2181 h1_x.ExchangeFaceNbrData();
2182 nd_x.ExchangeFaceNbrData();
2183 rt_x.ExchangeFaceNbrData();
2184 l2_x.ExchangeFaceNbrData();
2185 dgv_x.ExchangeFaceNbrData();
2186 dgi_x.ExchangeFaceNbrData();
2187
2188 Vector f_val(dim); f_val = 0.0;
2189
2190 Vector h1_gfc_val(dim); h1_gfc_val = 0.0;
2191 Vector nd_gfc_val(dim); nd_gfc_val = 0.0;
2192 Vector rt_gfc_val(dim); rt_gfc_val = 0.0;
2193 Vector l2_gfc_val(dim); l2_gfc_val = 0.0;
2194 Vector dgv_gfc_val(dim); dgv_gfc_val = 0.0;
2195 Vector dgi_gfc_val(dim); dgi_gfc_val = 0.0;
2196
2197 Vector h1_gvv_val(dim); h1_gvv_val = 0.0;
2198 Vector nd_gvv_val(dim); nd_gvv_val = 0.0;
2199 Vector rt_gvv_val(dim); rt_gvv_val = 0.0;
2200 Vector l2_gvv_val(dim); l2_gvv_val = 0.0;
2201 Vector dgv_gvv_val(dim); dgv_gvv_val = 0.0;
2202 Vector dgi_gvv_val(dim); dgi_gvv_val = 0.0;
2203
2204 SECTION("Shared Face Evaluation 2D")
2205 {
2206 if (my_rank == 0)
2207 {
2208 std::cout << "Shared Face Evaluation 2D" << std::endl;
2209 }
2210 for (int sf = 0; sf < pmesh.GetNSharedFaces(); sf++)
2211 {
2212 FaceElementTransformations *FET =
2213 pmesh.GetSharedFaceTransformations(sf);
2214 ElementTransformation *T = &FET->GetElement2Transformation();
2215 int e = FET->Elem2No;
2216 int e_nbr = e - pmesh.GetNE();
2217 const FiniteElement *fe = dgv_fespace.GetFaceNbrFE(e_nbr);
2218 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
2219 2*order + 2);
2220
2221 double h1_gfc_err = 0.0;
2222 double nd_gfc_err = 0.0;
2223 double rt_gfc_err = 0.0;
2224 double l2_gfc_err = 0.0;
2225 double dgv_gfc_err = 0.0;
2226 double dgi_gfc_err = 0.0;
2227
2228 double h1_gvv_err = 0.0;
2229 double nd_gvv_err = 0.0;
2230 double rt_gvv_err = 0.0;
2231 double l2_gvv_err = 0.0;
2232 double dgv_gvv_err = 0.0;
2233 double dgi_gvv_err = 0.0;
2234
2235 for (int j=0; j<ir.GetNPoints(); j++)
2236 {
2237 npts++;
2238 const IntegrationPoint &ip = ir.IntPoint(j);
2239 T->SetIntPoint(&ip);
2240
2241 funcCoef.Eval(f_val, *T, ip);
2242
2243 h1_xCoef.Eval(h1_gfc_val, *T, ip);
2244 nd_xCoef.Eval(nd_gfc_val, *T, ip);
2245 rt_xCoef.Eval(rt_gfc_val, *T, ip);
2246 l2_xCoef.Eval(l2_gfc_val, *T, ip);
2247 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
2248 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
2249
2250 h1_x.GetVectorValue(e, ip, h1_gvv_val);
2251 nd_x.GetVectorValue(e, ip, nd_gvv_val);
2252 rt_x.GetVectorValue(e, ip, rt_gvv_val);
2253 l2_x.GetVectorValue(e, ip, l2_gvv_val);
2254 dgv_x.GetVectorValue(e, ip, dgv_gvv_val);
2255 dgi_x.GetVectorValue(e, ip, dgi_gvv_val);
2256
2257 double h1_gfc_dist = Distance(f_val, h1_gfc_val, dim);
2258 double nd_gfc_dist = Distance(f_val, nd_gfc_val, dim);
2259 double rt_gfc_dist = Distance(f_val, rt_gfc_val, dim);
2260 double l2_gfc_dist = Distance(f_val, l2_gfc_val, dim);
2261 double dgv_gfc_dist = Distance(f_val, dgv_gfc_val, dim);
2262 double dgi_gfc_dist = Distance(f_val, dgi_gfc_val, dim);
2263
2264 double h1_gvv_dist = Distance(f_val, h1_gvv_val, dim);
2265 double nd_gvv_dist = Distance(f_val, nd_gvv_val, dim);
2266 double rt_gvv_dist = Distance(f_val, rt_gvv_val, dim);
2267 double l2_gvv_dist = Distance(f_val, l2_gvv_val, dim);
2268 double dgv_gvv_dist = Distance(f_val, dgv_gvv_val, dim);
2269 double dgi_gvv_dist = Distance(f_val, dgi_gvv_val, dim);
2270
2271 h1_gfc_err += h1_gfc_dist;
2272 nd_gfc_err += nd_gfc_dist;
2273 rt_gfc_err += rt_gfc_dist;
2274 l2_gfc_err += l2_gfc_dist;
2275 dgv_gfc_err += dgv_gfc_dist;
2276 dgi_gfc_err += dgi_gfc_dist;
2277
2278 h1_gvv_err += h1_gvv_dist;
2279 nd_gvv_err += nd_gvv_dist;
2280 rt_gvv_err += rt_gvv_dist;
2281 l2_gvv_err += l2_gvv_dist;
2282 dgv_gvv_err += dgv_gvv_dist;
2283 dgi_gvv_err += dgi_gvv_dist;
2284
2285 if (log > 0 && h1_gfc_dist > tol)
2286 {
2287 std::cout << e << ":" << j << " h1 ("
2288 << f_val[0] << "," << f_val[1] << ") vs. ("
2289 << h1_gfc_val[0] << ","
2290 << h1_gfc_val[1] << ") "
2291 << h1_gfc_dist << std::endl;
2292 }
2293 if (log > 0 && nd_gfc_dist > tol)
2294 {
2295 std::cout << e << ":" << j << " nd ("
2296 << f_val[0] << "," << f_val[1] << ") vs. ("
2297 << nd_gfc_val[0] << ","
2298 << nd_gfc_val[1] << ") "
2299 << nd_gfc_dist << std::endl;
2300 }
2301 if (log > 0 && rt_gfc_dist > tol)
2302 {
2303 std::cout << e << ":" << j << " rt ("
2304 << f_val[0] << "," << f_val[1] << ") vs. ("
2305 << rt_gfc_val[0] << ","
2306 << rt_gfc_val[1] << ") "
2307 << rt_gfc_dist << std::endl;
2308 }
2309 if (log > 0 && l2_gfc_dist > tol)
2310 {
2311 std::cout << e << ":" << j << " l2 ("
2312 << f_val[0] << "," << f_val[1] << ") vs. ("
2313 << l2_gfc_val[0] << ","
2314 << l2_gfc_val[1] << ") "
2315 << l2_gfc_dist << std::endl;
2316 }
2317 if (log > 0 && dgv_gfc_dist > tol)
2318 {
2319 std::cout << e << ":" << j << " dgv ("
2320 << f_val[0] << "," << f_val[1] << ") vs. ("
2321 << dgv_gfc_val[0] << ","
2322 << dgv_gfc_val[1] << ") "
2323 << dgv_gfc_dist << std::endl;
2324 }
2325 if (log > 0 && dgi_gfc_dist > tol)
2326 {
2327 std::cout << e << ":" << j << " dgi ("
2328 << f_val[0] << "," << f_val[1] << ") vs. ("
2329 << dgi_gfc_val[0] << ","
2330 << dgi_gfc_val[1] << ") "
2331 << dgi_gfc_dist << std::endl;
2332 }
2333 if (log > 0 && h1_gvv_dist > tol)
2334 {
2335 std::cout << e << ":" << j << " h1 gvv ("
2336 << f_val[0] << "," << f_val[1] << ","
2337 << f_val[2] << ") vs. ("
2338 << h1_gvv_val[0] << "," << h1_gvv_val[1] << ") "
2339 << h1_gvv_dist << std::endl;
2340 }
2341 if (log > 0 && nd_gvv_dist > tol)
2342 {
2343 std::cout << e << ":" << j << " nd gvv ("
2344 << f_val[0] << "," << f_val[1] << ","
2345 << f_val[2] << ") vs. ("
2346 << nd_gvv_val[0] << "," << nd_gvv_val[1] << ") "
2347 << nd_gvv_dist << std::endl;
2348 }
2349 if (log > 0 && rt_gvv_dist > tol)
2350 {
2351 std::cout << e << ":" << j << " rt gvv ("
2352 << f_val[0] << "," << f_val[1] << ","
2353 << f_val[2] << ") vs. ("
2354 << rt_gvv_val[0] << "," << rt_gvv_val[1] << ") "
2355 << rt_gvv_dist << std::endl;
2356 }
2357 if (log > 0 && l2_gvv_dist > tol)
2358 {
2359 std::cout << e << ":" << j << " l2 gvv ("
2360 << f_val[0] << "," << f_val[1] << ","
2361 << f_val[2] << ") vs. ("
2362 << l2_gvv_val[0] << "," << l2_gvv_val[1] << ") "
2363 << l2_gvv_dist << std::endl;
2364 }
2365 if (log > 0 && dgv_gvv_dist > tol)
2366 {
2367 std::cout << e << ":" << j << " dgv gvv ("
2368 << f_val[0] << "," << f_val[1] << ","
2369 << f_val[2] << ") vs. ("
2370 << dgv_gvv_val[0] << ","
2371 << dgv_gvv_val[1] << ") "
2372 << dgv_gvv_dist << std::endl;
2373 }
2374 if (log > 0 && dgi_gvv_dist > tol)
2375 {
2376 std::cout << e << ":" << j << " dgi gvv ("
2377 << f_val[0] << "," << f_val[1] << ","
2378 << f_val[2] << ") vs. ("
2379 << dgi_gvv_val[0] << ","
2380 << dgi_gvv_val[1] << ") "
2381 << dgi_gvv_dist << std::endl;
2382 }
2383 }
2384
2385 h1_gfc_err /= ir.GetNPoints();
2386 nd_gfc_err /= ir.GetNPoints();
2387 rt_gfc_err /= ir.GetNPoints();
2388 l2_gfc_err /= ir.GetNPoints();
2389 dgv_gfc_err /= ir.GetNPoints();
2390 dgi_gfc_err /= ir.GetNPoints();
2391
2392 h1_gvv_err /= ir.GetNPoints();
2393 nd_gvv_err /= ir.GetNPoints();
2394 rt_gvv_err /= ir.GetNPoints();
2395 l2_gvv_err /= ir.GetNPoints();
2396 dgv_gvv_err /= ir.GetNPoints();
2397 dgi_gvv_err /= ir.GetNPoints();
2398
2399 REQUIRE( h1_gfc_err == MFEM_Approx(0.0));
2400 REQUIRE( nd_gfc_err == MFEM_Approx(0.0));
2401 REQUIRE( rt_gfc_err == MFEM_Approx(0.0));
2402 REQUIRE( l2_gfc_err == MFEM_Approx(0.0));
2403 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
2404 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
2405
2406 REQUIRE( h1_gvv_err == MFEM_Approx(0.0));
2407 REQUIRE( nd_gvv_err == MFEM_Approx(0.0));
2408 REQUIRE( rt_gvv_err == MFEM_Approx(0.0));
2409 REQUIRE( l2_gvv_err == MFEM_Approx(0.0));
2410 REQUIRE(dgv_gvv_err == MFEM_Approx(0.0));
2411 REQUIRE(dgi_gvv_err == MFEM_Approx(0.0));
2412 }
2413 }
2414 }
2415 }
2416 std::cout << my_rank << ": Checked GridFunction::GetVectorValue at "
2417 << npts << " 2D points" << std::endl;
2418 }
2419
2420 #endif // MFEM_USE_MPI
2421
2422 TEST_CASE("3D GetVectorValue",
2423 "[GridFunction]"
2424 "[VectorGridFunctionCoefficient]")
2425 {
2426 int log = 1;
2427 int n = 1;
2428 int dim = 3;
2429 int order = 1;
2430 int npts = 0;
2431
2432 double tol = 1e-6;
2433
2434 for (int type = (int)Element::TETRAHEDRON;
2435 type <= (int)Element::HEXAHEDRON; type++)
2436 {
2437 Mesh mesh = Mesh::MakeCartesian3D(
2438 n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
2439
2440 VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
2441
2442 SECTION("3D GetVectorValue tests for element type " +
to_string(type)2443 std::to_string(type))
2444 {
2445 H1_FECollection h1_fec(order, dim);
2446 ND_FECollection nd_fec(order+1, dim);
2447 RT_FECollection rt_fec(order+1, dim);
2448 L2_FECollection l2_fec(order, dim);
2449 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
2450 FiniteElement::VALUE);
2451 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
2452 FiniteElement::INTEGRAL);
2453
2454 FiniteElementSpace h1_fespace(&mesh, &h1_fec, dim);
2455 FiniteElementSpace nd_fespace(&mesh, &nd_fec);
2456 FiniteElementSpace rt_fespace(&mesh, &rt_fec);
2457 FiniteElementSpace l2_fespace(&mesh, &l2_fec, dim);
2458 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec, dim);
2459 FiniteElementSpace dgi_fespace(&mesh, &dgi_fec, dim);
2460
2461 GridFunction h1_x( &h1_fespace);
2462 GridFunction nd_x( &nd_fespace);
2463 GridFunction rt_x( &rt_fespace);
2464 GridFunction l2_x( &l2_fespace);
2465 GridFunction dgv_x(&dgv_fespace);
2466 GridFunction dgi_x(&dgi_fespace);
2467
2468 VectorGridFunctionCoefficient h1_xCoef( &h1_x);
2469 VectorGridFunctionCoefficient nd_xCoef( &nd_x);
2470 VectorGridFunctionCoefficient rt_xCoef( &rt_x);
2471 VectorGridFunctionCoefficient l2_xCoef( &l2_x);
2472 VectorGridFunctionCoefficient dgv_xCoef(&dgv_x);
2473 VectorGridFunctionCoefficient dgi_xCoef(&dgi_x);
2474
2475 h1_x.ProjectCoefficient(funcCoef);
2476 nd_x.ProjectCoefficient(funcCoef);
2477 rt_x.ProjectCoefficient(funcCoef);
2478 l2_x.ProjectCoefficient(funcCoef);
2479 dgv_x.ProjectCoefficient(funcCoef);
2480 dgi_x.ProjectCoefficient(funcCoef);
2481
2482 Vector f_val(dim); f_val = 0.0;
2483
2484 Vector h1_gfc_val(dim); h1_gfc_val = 0.0;
2485 Vector nd_gfc_val(dim); nd_gfc_val = 0.0;
2486 Vector rt_gfc_val(dim); rt_gfc_val = 0.0;
2487 Vector l2_gfc_val(dim); l2_gfc_val = 0.0;
2488 Vector dgv_gfc_val(dim); dgv_gfc_val = 0.0;
2489 Vector dgi_gfc_val(dim); dgi_gfc_val = 0.0;
2490
2491 Vector h1_gvv_val(dim); h1_gvv_val = 0.0;
2492 Vector nd_gvv_val(dim); nd_gvv_val = 0.0;
2493 Vector rt_gvv_val(dim); rt_gvv_val = 0.0;
2494 Vector l2_gvv_val(dim); l2_gvv_val = 0.0;
2495 Vector dgv_gvv_val(dim); dgv_gvv_val = 0.0;
2496 Vector dgi_gvv_val(dim); dgi_gvv_val = 0.0;
2497
2498 SECTION("Domain Evaluation 3D")
2499 {
2500 std::cout << "Domain Evaluation 3D" << std::endl;
2501 for (int e = 0; e < mesh.GetNE(); e++)
2502 {
2503 ElementTransformation *T = mesh.GetElementTransformation(e);
2504 const FiniteElement *fe = h1_fespace.GetFE(e);
2505 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
2506 2*order + 2);
2507
2508 double h1_gfc_err = 0.0;
2509 double nd_gfc_err = 0.0;
2510 double rt_gfc_err = 0.0;
2511 double l2_gfc_err = 0.0;
2512 double dgv_gfc_err = 0.0;
2513 double dgi_gfc_err = 0.0;
2514
2515 double h1_gvv_err = 0.0;
2516 double nd_gvv_err = 0.0;
2517 double rt_gvv_err = 0.0;
2518 double l2_gvv_err = 0.0;
2519 double dgv_gvv_err = 0.0;
2520 double dgi_gvv_err = 0.0;
2521
2522 for (int j=0; j<ir.GetNPoints(); j++)
2523 {
2524 npts++;
2525 const IntegrationPoint &ip = ir.IntPoint(j);
2526 T->SetIntPoint(&ip);
2527
2528 funcCoef.Eval(f_val, *T, ip);
2529
2530 h1_xCoef.Eval(h1_gfc_val, *T, ip);
2531 nd_xCoef.Eval(nd_gfc_val, *T, ip);
2532 rt_xCoef.Eval(rt_gfc_val, *T, ip);
2533 l2_xCoef.Eval(l2_gfc_val, *T, ip);
2534 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
2535 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
2536
2537 h1_x.GetVectorValue(e, ip, h1_gvv_val);
2538 nd_x.GetVectorValue(e, ip, nd_gvv_val);
2539 rt_x.GetVectorValue(e, ip, rt_gvv_val);
2540 l2_x.GetVectorValue(e, ip, l2_gvv_val);
2541 dgv_x.GetVectorValue(e, ip, dgv_gvv_val);
2542 dgi_x.GetVectorValue(e, ip, dgi_gvv_val);
2543
2544 double h1_gfc_dist = Distance(f_val, h1_gfc_val, dim);
2545 double nd_gfc_dist = Distance(f_val, nd_gfc_val, dim);
2546 double rt_gfc_dist = Distance(f_val, rt_gfc_val, dim);
2547 double l2_gfc_dist = Distance(f_val, l2_gfc_val, dim);
2548 double dgv_gfc_dist = Distance(f_val, dgv_gfc_val, dim);
2549 double dgi_gfc_dist = Distance(f_val, dgi_gfc_val, dim);
2550
2551 double h1_gvv_dist = Distance(f_val, h1_gvv_val, dim);
2552 double nd_gvv_dist = Distance(f_val, nd_gvv_val, dim);
2553 double rt_gvv_dist = Distance(f_val, rt_gvv_val, dim);
2554 double l2_gvv_dist = Distance(f_val, l2_gvv_val, dim);
2555 double dgv_gvv_dist = Distance(f_val, dgv_gvv_val, dim);
2556 double dgi_gvv_dist = Distance(f_val, dgi_gvv_val, dim);
2557
2558 h1_gfc_err += h1_gfc_dist;
2559 nd_gfc_err += nd_gfc_dist;
2560 rt_gfc_err += rt_gfc_dist;
2561 l2_gfc_err += l2_gfc_dist;
2562 dgv_gfc_err += dgv_gfc_dist;
2563 dgi_gfc_err += dgi_gfc_dist;
2564
2565 h1_gvv_err += h1_gvv_dist;
2566 nd_gvv_err += nd_gvv_dist;
2567 rt_gvv_err += rt_gvv_dist;
2568 l2_gvv_err += l2_gvv_dist;
2569 dgv_gvv_err += dgv_gvv_dist;
2570 dgi_gvv_err += dgi_gvv_dist;
2571
2572 if (log > 0 && h1_gfc_dist > tol)
2573 {
2574 std::cout << e << ":" << j << " h1 gfc ("
2575 << f_val[0] << "," << f_val[1] << ","
2576 << f_val[2] << ") vs. ("
2577 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
2578 << h1_gfc_val[2] << ") "
2579 << h1_gfc_dist << std::endl;
2580 }
2581 if (log > 0 && nd_gfc_dist > tol)
2582 {
2583 std::cout << e << ":" << j << " nd gfc ("
2584 << f_val[0] << "," << f_val[1] << ","
2585 << f_val[2] << ") vs. ("
2586 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ","
2587 << nd_gfc_val[2] << ") "
2588 << nd_gfc_dist << std::endl;
2589 }
2590 if (log > 0 && rt_gfc_dist > tol)
2591 {
2592 std::cout << e << ":" << j << " rt gfc ("
2593 << f_val[0] << "," << f_val[1] << ","
2594 << f_val[2] << ") vs. ("
2595 << rt_gfc_val[0] << "," << rt_gfc_val[1] << ","
2596 << rt_gfc_val[2] << ") "
2597 << rt_gfc_dist << std::endl;
2598 }
2599 if (log > 0 && l2_gfc_dist > tol)
2600 {
2601 std::cout << e << ":" << j << " l2 gfc ("
2602 << f_val[0] << "," << f_val[1] << ","
2603 << f_val[2] << ") vs. ("
2604 << l2_gfc_val[0] << "," << l2_gfc_val[1] << ","
2605 << l2_gfc_val[2] << ") "
2606 << l2_gfc_dist << std::endl;
2607 }
2608 if (log > 0 && dgv_gfc_dist > tol)
2609 {
2610 std::cout << e << ":" << j << " dgv gfc ("
2611 << f_val[0] << "," << f_val[1] << ","
2612 << f_val[2] << ") vs. ("
2613 << dgv_gfc_val[0] << ","
2614 << dgv_gfc_val[1] << ","
2615 << dgv_gfc_val[2] << ") "
2616 << dgv_gfc_dist << std::endl;
2617 }
2618 if (log > 0 && dgi_gfc_dist > tol)
2619 {
2620 std::cout << e << ":" << j << " dgi gfc ("
2621 << f_val[0] << "," << f_val[1] << ","
2622 << f_val[2] << ") vs. ("
2623 << dgi_gfc_val[0] << ","
2624 << dgi_gfc_val[1] << ","
2625 << dgi_gfc_val[2] << ") "
2626 << dgi_gfc_dist << std::endl;
2627 }
2628 if (log > 0 && h1_gvv_dist > tol)
2629 {
2630 std::cout << e << ":" << j << " h1 gvv ("
2631 << f_val[0] << "," << f_val[1] << ","
2632 << f_val[2] << ") vs. ("
2633 << h1_gvv_val[0] << "," << h1_gvv_val[1] << ","
2634 << h1_gvv_val[2] << ") "
2635 << h1_gvv_dist << std::endl;
2636 }
2637 if (log > 0 && nd_gvv_dist > tol)
2638 {
2639 std::cout << e << ":" << j << " nd gvv ("
2640 << f_val[0] << "," << f_val[1] << ","
2641 << f_val[2] << ") vs. ("
2642 << nd_gvv_val[0] << "," << nd_gvv_val[1] << ","
2643 << nd_gvv_val[2] << ") "
2644 << nd_gvv_dist << std::endl;
2645 }
2646 if (log > 0 && rt_gvv_dist > tol)
2647 {
2648 std::cout << e << ":" << j << " rt gvv ("
2649 << f_val[0] << "," << f_val[1] << ","
2650 << f_val[2] << ") vs. ("
2651 << rt_gvv_val[0] << "," << rt_gvv_val[1] << ","
2652 << rt_gvv_val[2] << ") "
2653 << rt_gvv_dist << std::endl;
2654 }
2655 if (log > 0 && l2_gvv_dist > tol)
2656 {
2657 std::cout << e << ":" << j << " l2 gvv ("
2658 << f_val[0] << "," << f_val[1] << ","
2659 << f_val[2] << ") vs. ("
2660 << l2_gvv_val[0] << "," << l2_gvv_val[1] << ","
2661 << l2_gvv_val[2] << ") "
2662 << l2_gvv_dist << std::endl;
2663 }
2664 if (log > 0 && dgv_gvv_dist > tol)
2665 {
2666 std::cout << e << ":" << j << " dgv gvv ("
2667 << f_val[0] << "," << f_val[1] << ","
2668 << f_val[2] << ") vs. ("
2669 << dgv_gvv_val[0] << ","
2670 << dgv_gvv_val[1] << ","
2671 << dgv_gvv_val[2] << ") "
2672 << dgv_gvv_dist << std::endl;
2673 }
2674 if (log > 0 && dgi_gvv_dist > tol)
2675 {
2676 std::cout << e << ":" << j << " dgi gvv ("
2677 << f_val[0] << "," << f_val[1] << ","
2678 << f_val[2] << ") vs. ("
2679 << dgi_gvv_val[0] << ","
2680 << dgi_gvv_val[1] << ","
2681 << dgi_gvv_val[2] << ") "
2682 << dgi_gvv_dist << std::endl;
2683 }
2684 }
2685
2686 h1_gfc_err /= ir.GetNPoints();
2687 nd_gfc_err /= ir.GetNPoints();
2688 rt_gfc_err /= ir.GetNPoints();
2689 l2_gfc_err /= ir.GetNPoints();
2690 dgv_gfc_err /= ir.GetNPoints();
2691 dgi_gfc_err /= ir.GetNPoints();
2692
2693 h1_gvv_err /= ir.GetNPoints();
2694 nd_gvv_err /= ir.GetNPoints();
2695 rt_gvv_err /= ir.GetNPoints();
2696 l2_gvv_err /= ir.GetNPoints();
2697 dgv_gvv_err /= ir.GetNPoints();
2698 dgi_gvv_err /= ir.GetNPoints();
2699
2700 REQUIRE( h1_gfc_err == MFEM_Approx(0.0));
2701 REQUIRE( nd_gfc_err == MFEM_Approx(0.0));
2702 REQUIRE( rt_gfc_err == MFEM_Approx(0.0));
2703 REQUIRE( l2_gfc_err == MFEM_Approx(0.0));
2704 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
2705 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
2706
2707 REQUIRE( h1_gvv_err == MFEM_Approx(0.0));
2708 REQUIRE( nd_gvv_err == MFEM_Approx(0.0));
2709 REQUIRE( rt_gvv_err == MFEM_Approx(0.0));
2710 REQUIRE( l2_gvv_err == MFEM_Approx(0.0));
2711 REQUIRE(dgv_gvv_err == MFEM_Approx(0.0));
2712 REQUIRE(dgi_gvv_err == MFEM_Approx(0.0));
2713 }
2714 }
2715
2716 SECTION("Boundary Evaluation 3D (H1 Context)")
2717 {
2718 std::cout << "Boundary Evaluation 3D (H1 Context)" << std::endl;
2719 for (int be = 0; be < mesh.GetNBE(); be++)
2720 {
2721 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
2722 const FiniteElement *fe = h1_fespace.GetBE(be);
2723 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
2724 2*order + 2);
2725
2726 double h1_err = 0.0;
2727 double nd_err = 0.0;
2728 double rt_err = 0.0;
2729 double l2_err = 0.0;
2730 double dgv_err = 0.0;
2731 double dgi_err = 0.0;
2732
2733 for (int j=0; j<ir.GetNPoints(); j++)
2734 {
2735 npts++;
2736 const IntegrationPoint &ip = ir.IntPoint(j);
2737 T->SetIntPoint(&ip);
2738
2739 funcCoef.Eval(f_val, *T, ip);
2740 h1_xCoef.Eval(h1_gfc_val, *T, ip);
2741 nd_xCoef.Eval(nd_gfc_val, *T, ip);
2742 rt_xCoef.Eval(rt_gfc_val, *T, ip);
2743 l2_xCoef.Eval(l2_gfc_val, *T, ip);
2744 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
2745 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
2746
2747 double h1_dist = Distance(f_val, h1_gfc_val, dim);
2748 double nd_dist = Distance(f_val, nd_gfc_val, dim);
2749 double rt_dist = Distance(f_val, rt_gfc_val, dim);
2750 double l2_dist = Distance(f_val, l2_gfc_val, dim);
2751 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
2752 double dgi_dist = Distance(f_val, dgi_gfc_val, dim);
2753
2754 h1_err += h1_dist;
2755 nd_err += nd_dist;
2756 rt_err += rt_dist;
2757 l2_err += l2_dist;
2758 dgv_err += dgv_dist;
2759 dgi_err += dgi_dist;
2760
2761 if (log > 0 && h1_dist > tol)
2762 {
2763 std::cout << be << ":" << j << " h1 ("
2764 << f_val[0] << "," << f_val[1] << ","
2765 << f_val[2] << ") vs. ("
2766 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
2767 << h1_gfc_val[2] << ") " << h1_dist
2768 << std::endl;
2769 }
2770 if (log > 0 && nd_dist > tol)
2771 {
2772 std::cout << be << ":" << j << " nd ("
2773 << f_val[0] << "," << f_val[1] << ","
2774 << f_val[2] << ") vs. ("
2775 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ","
2776 << nd_gfc_val[2] << ") " << nd_dist
2777 << std::endl;
2778 }
2779 if (log > 0 && rt_dist > tol)
2780 {
2781 std::cout << be << ":" << j << " rt ("
2782 << f_val[0] << "," << f_val[1] << ","
2783 << f_val[2] << ") vs. ("
2784 << rt_gfc_val[0] << "," << rt_gfc_val[1] << ","
2785 << rt_gfc_val[2] << ") " << rt_dist
2786 << std::endl;
2787 }
2788 if (log > 0 && l2_dist > tol)
2789 {
2790 std::cout << be << ":" << j << " l2 ("
2791 << f_val[0] << "," << f_val[1] << ","
2792 << f_val[2] << ") vs. ("
2793 << l2_gfc_val[0] << "," << l2_gfc_val[1] << ","
2794 << l2_gfc_val[2] << ") " << l2_dist
2795 << std::endl;
2796 }
2797 if (log > 0 && dgv_dist > tol)
2798 {
2799 std::cout << be << ":" << j << " dgv ("
2800 << f_val[0] << "," << f_val[1] << ","
2801 << f_val[2] << ") vs. ("
2802 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
2803 << dgv_gfc_val[2] << ") " << dgv_dist
2804 << std::endl;
2805 }
2806 if (log > 0 && dgi_dist > tol)
2807 {
2808 std::cout << be << ":" << j << " dgi ("
2809 << f_val[0] << "," << f_val[1] << ","
2810 << f_val[2] << ") vs. ("
2811 << dgi_gfc_val[0] << "," << dgi_gfc_val[1] << ","
2812 << dgi_gfc_val[2] << ") " << dgi_dist
2813 << std::endl;
2814 }
2815 }
2816 h1_err /= ir.GetNPoints();
2817 nd_err /= ir.GetNPoints();
2818 rt_err /= ir.GetNPoints();
2819 l2_err /= ir.GetNPoints();
2820 dgv_err /= ir.GetNPoints();
2821 dgi_err /= ir.GetNPoints();
2822
2823 REQUIRE( h1_err == MFEM_Approx(0.0));
2824 REQUIRE( nd_err == MFEM_Approx(0.0));
2825 REQUIRE( rt_err == MFEM_Approx(0.0));
2826 REQUIRE( l2_err == MFEM_Approx(0.0));
2827 REQUIRE(dgv_err == MFEM_Approx(0.0));
2828 REQUIRE(dgi_err == MFEM_Approx(0.0));
2829 }
2830 }
2831
2832 SECTION("Boundary Evaluation 3D (DG Context)")
2833 {
2834 std::cout << "Boundary Evaluation 3D (DG Context)" << std::endl;
2835 for (int be = 0; be < mesh.GetNBE(); be++)
2836 {
2837 FaceElementTransformations *T =
2838 mesh.GetBdrFaceTransformations(be);
2839 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
2840 2*order + 2);
2841
2842 double h1_err = 0.0;
2843 double nd_err = 0.0;
2844 double rt_err = 0.0;
2845 double l2_err = 0.0;
2846 double dgv_err = 0.0;
2847 double dgi_err = 0.0;
2848
2849 for (int j=0; j<ir.GetNPoints(); j++)
2850 {
2851 npts++;
2852 const IntegrationPoint &ip = ir.IntPoint(j);
2853
2854 T->SetIntPoint(&ip);
2855
2856 funcCoef.Eval(f_val, *T, ip);
2857 h1_xCoef.Eval(h1_gfc_val, *T, ip);
2858 nd_xCoef.Eval(nd_gfc_val, *T, ip);
2859 rt_xCoef.Eval(rt_gfc_val, *T, ip);
2860 l2_xCoef.Eval(l2_gfc_val, *T, ip);
2861 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
2862 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
2863
2864 double h1_dist = Distance(f_val, h1_gfc_val, dim);
2865 double nd_dist = Distance(f_val, nd_gfc_val, dim);
2866 double rt_dist = Distance(f_val, rt_gfc_val, dim);
2867 double l2_dist = Distance(f_val, l2_gfc_val, dim);
2868 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
2869 double dgi_dist = Distance(f_val, dgi_gfc_val, dim);
2870
2871 h1_err += h1_dist;
2872 nd_err += nd_dist;
2873 rt_err += rt_dist;
2874 l2_err += l2_dist;
2875 dgv_err += dgv_dist;
2876 dgi_err += dgi_dist;
2877
2878 if (log > 0 && h1_dist > tol)
2879 {
2880 std::cout << be << ":" << j << " h1 ("
2881 << f_val[0] << "," << f_val[1] << ","
2882 << f_val[2] << ") vs. ("
2883 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
2884 << h1_gfc_val[2] << ") " << h1_dist
2885 << std::endl;
2886 }
2887 if (log > 0 && nd_dist > tol)
2888 {
2889 std::cout << be << ":" << j << " nd ("
2890 << f_val[0] << "," << f_val[1] << ","
2891 << f_val[2] << ") vs. ("
2892 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ","
2893 << nd_gfc_val[2] << ") " << nd_dist
2894 << std::endl;
2895 }
2896 if (log > 0 && rt_dist > tol)
2897 {
2898 std::cout << be << ":" << j << " rt ("
2899 << f_val[0] << "," << f_val[1] << ","
2900 << f_val[2] << ") vs. ("
2901 << rt_gfc_val[0] << "," << rt_gfc_val[1] << ","
2902 << rt_gfc_val[2] << ") " << rt_dist
2903 << std::endl;
2904 }
2905 if (log > 0 && l2_dist > tol)
2906 {
2907 std::cout << be << ":" << j << " l2 ("
2908 << f_val[0] << "," << f_val[1] << ","
2909 << f_val[2] << ") vs. ("
2910 << l2_gfc_val[0] << "," << l2_gfc_val[1] << ","
2911 << l2_gfc_val[2] << ") " << l2_dist
2912 << std::endl;
2913 }
2914 if (log > 0 && dgv_dist > tol)
2915 {
2916 std::cout << be << ":" << j << " dgv ("
2917 << f_val[0] << "," << f_val[1] << ","
2918 << f_val[2] << ") vs. ("
2919 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
2920 << dgv_gfc_val[2] << ") " << dgv_dist
2921 << std::endl;
2922 }
2923 if (log > 0 && dgi_dist > tol)
2924 {
2925 std::cout << be << ":" << j << " dgi ("
2926 << f_val[0] << "," << f_val[1] << ","
2927 << f_val[2] << ") vs. ("
2928 << dgi_gfc_val[0] << "," << dgi_gfc_val[1] << ","
2929 << dgi_gfc_val[2] << ") " << dgi_dist
2930 << std::endl;
2931 }
2932 }
2933 h1_err /= ir.GetNPoints();
2934 nd_err /= ir.GetNPoints();
2935 rt_err /= ir.GetNPoints();
2936 l2_err /= ir.GetNPoints();
2937 dgv_err /= ir.GetNPoints();
2938 dgi_err /= ir.GetNPoints();
2939
2940 REQUIRE( h1_err == MFEM_Approx(0.0));
2941 REQUIRE( nd_err == MFEM_Approx(0.0));
2942 REQUIRE( rt_err == MFEM_Approx(0.0));
2943 REQUIRE( l2_err == MFEM_Approx(0.0));
2944 REQUIRE(dgv_err == MFEM_Approx(0.0));
2945 REQUIRE(dgi_err == MFEM_Approx(0.0));
2946 }
2947 }
2948
2949 SECTION("Edge Evaluation 3D")
2950 {
2951 std::cout << "Edge Evaluation 3D" << std::endl;
2952 for (int e = 0; e < mesh.GetNEdges(); e++)
2953 {
2954 ElementTransformation *T = mesh.GetEdgeTransformation(e);
2955 const FiniteElement *fe = h1_fespace.GetEdgeElement(e);
2956 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
2957 2*order + 2);
2958
2959 double h1_err = 0.0;
2960
2961 for (int j=0; j<ir.GetNPoints(); j++)
2962 {
2963 npts++;
2964 const IntegrationPoint &ip = ir.IntPoint(j);
2965 T->SetIntPoint(&ip);
2966
2967 funcCoef.Eval(f_val, *T, ip);
2968 h1_xCoef.Eval(h1_gfc_val, *T, ip);
2969
2970 double h1_dist = Distance(f_val, h1_gfc_val, dim);
2971
2972 h1_err += h1_dist;
2973
2974 if (log > 0 && h1_dist > tol)
2975 {
2976 std::cout << e << ":" << j << " h1 ("
2977 << f_val[0] << "," << f_val[1] << ","
2978 << f_val[2] << ") vs. ("
2979 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
2980 << h1_gfc_val[2] << ") " << h1_dist
2981 << std::endl;
2982 }
2983 }
2984 h1_err /= ir.GetNPoints();
2985
2986 REQUIRE( h1_err == MFEM_Approx(0.0));
2987 }
2988 }
2989
2990 SECTION("Face Evaluation 3D")
2991 {
2992 std::cout << "Face Evaluation 3D" << std::endl;
2993 for (int f = 0; f < mesh.GetNFaces(); f++)
2994 {
2995 ElementTransformation *T = mesh.GetFaceTransformation(f);
2996 const FiniteElement *fe = h1_fespace.GetFaceElement(f);
2997 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
2998 2*order + 2);
2999
3000 double h1_err = 0.0;
3001
3002 for (int j=0; j<ir.GetNPoints(); j++)
3003 {
3004 npts++;
3005 const IntegrationPoint &ip = ir.IntPoint(j);
3006 T->SetIntPoint(&ip);
3007
3008 funcCoef.Eval(f_val, *T, ip);
3009 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3010
3011 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3012
3013 h1_err += h1_dist;
3014
3015 if (log > 0 && h1_dist > tol)
3016 {
3017 std::cout << f << ":" << j << " h1 ("
3018 << f_val[0] << "," << f_val[1] << ","
3019 << f_val[2] << ") vs. ("
3020 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
3021 << h1_gfc_val[2] << ") " << h1_dist
3022 << std::endl;
3023 }
3024 }
3025 h1_err /= ir.GetNPoints();
3026
3027 REQUIRE( h1_err == MFEM_Approx(0.0));
3028 }
3029 }
3030 }
3031 }
3032 std::cout << "Checked GridFunction::GetVectorValue at "
3033 << npts << " 3D points" << std::endl;
3034 }
3035
3036 #ifdef MFEM_USE_MPI
3037 #
3038 TEST_CASE("3D GetVectorValue in Parallel",
3039 "[ParGridFunction]"
3040 "[VectorGridFunctionCoefficient]"
3041 "[Parallel]")
3042 {
3043 int num_procs;
3044 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
3045
3046 int my_rank;
3047 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
3048
3049 int log = 1;
3050 int n = (int)ceil(pow(2*num_procs, 1.0 / 3.0));
3051 int dim = 3;
3052 int order = 1;
3053 int npts = 0;
3054
3055 double tol = 1e-6;
3056
3057 for (int type = (int)Element::TETRAHEDRON;
3058 type <= (int)Element::HEXAHEDRON; type++)
3059 {
3060 Mesh mesh = Mesh::MakeCartesian3D(
3061 n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
3062 ParMesh pmesh(MPI_COMM_WORLD, mesh);
3063 if (type == Element::TETRAHEDRON)
3064 {
3065 pmesh.ReorientTetMesh();
3066 }
3067 mesh.Clear();
3068
3069 VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
3070
3071 SECTION("3D GetVectorValue tests for element type " +
to_string(type)3072 std::to_string(type))
3073 {
3074 H1_FECollection h1_fec(order, dim);
3075 ND_FECollection nd_fec(order+1, dim);
3076 RT_FECollection rt_fec(order+1, dim);
3077 L2_FECollection l2_fec(order, dim);
3078 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
3079 FiniteElement::VALUE);
3080 DG_FECollection dgi_fec(order, dim, BasisType::GaussLegendre,
3081 FiniteElement::INTEGRAL);
3082
3083 ParFiniteElementSpace h1_fespace(&pmesh, &h1_fec, dim);
3084 ParFiniteElementSpace nd_fespace(&pmesh, &nd_fec);
3085 ParFiniteElementSpace rt_fespace(&pmesh, &rt_fec);
3086 ParFiniteElementSpace l2_fespace(&pmesh, &l2_fec, dim);
3087 ParFiniteElementSpace dgv_fespace(&pmesh, &dgv_fec, dim);
3088 ParFiniteElementSpace dgi_fespace(&pmesh, &dgi_fec, dim);
3089
3090 ParGridFunction h1_x( &h1_fespace);
3091 ParGridFunction nd_x( &nd_fespace);
3092 ParGridFunction rt_x( &rt_fespace);
3093 ParGridFunction l2_x( &l2_fespace);
3094 ParGridFunction dgv_x(&dgv_fespace);
3095 ParGridFunction dgi_x(&dgi_fespace);
3096
3097 VectorGridFunctionCoefficient h1_xCoef( &h1_x);
3098 VectorGridFunctionCoefficient nd_xCoef( &nd_x);
3099 VectorGridFunctionCoefficient rt_xCoef( &rt_x);
3100 VectorGridFunctionCoefficient l2_xCoef( &l2_x);
3101 VectorGridFunctionCoefficient dgv_xCoef(&dgv_x);
3102 VectorGridFunctionCoefficient dgi_xCoef(&dgi_x);
3103
3104 h1_x.ProjectCoefficient(funcCoef);
3105 nd_x.ProjectCoefficient(funcCoef);
3106 rt_x.ProjectCoefficient(funcCoef);
3107 l2_x.ProjectCoefficient(funcCoef);
3108 dgv_x.ProjectCoefficient(funcCoef);
3109 dgi_x.ProjectCoefficient(funcCoef);
3110
3111 h1_x.ExchangeFaceNbrData();
3112 nd_x.ExchangeFaceNbrData();
3113 rt_x.ExchangeFaceNbrData();
3114 l2_x.ExchangeFaceNbrData();
3115 dgv_x.ExchangeFaceNbrData();
3116 dgi_x.ExchangeFaceNbrData();
3117
3118 Vector f_val(dim); f_val = 0.0;
3119
3120 Vector h1_gfc_val(dim); h1_gfc_val = 0.0;
3121 Vector nd_gfc_val(dim); nd_gfc_val = 0.0;
3122 Vector rt_gfc_val(dim); rt_gfc_val = 0.0;
3123 Vector l2_gfc_val(dim); l2_gfc_val = 0.0;
3124 Vector dgv_gfc_val(dim); dgv_gfc_val = 0.0;
3125 Vector dgi_gfc_val(dim); dgi_gfc_val = 0.0;
3126
3127 Vector h1_gvv_val(dim); h1_gvv_val = 0.0;
3128 Vector nd_gvv_val(dim); nd_gvv_val = 0.0;
3129 Vector rt_gvv_val(dim); rt_gvv_val = 0.0;
3130 Vector l2_gvv_val(dim); l2_gvv_val = 0.0;
3131 Vector dgv_gvv_val(dim); dgv_gvv_val = 0.0;
3132 Vector dgi_gvv_val(dim); dgi_gvv_val = 0.0;
3133
3134 SECTION("Shared Face Evaluation 3D")
3135 {
3136 if (my_rank == 0)
3137 {
3138 std::cout << "Shared Face Evaluation 3D" << std::endl;
3139 }
3140 for (int sf = 0; sf < pmesh.GetNSharedFaces(); sf++)
3141 {
3142 FaceElementTransformations *FET =
3143 pmesh.GetSharedFaceTransformations(sf);
3144 ElementTransformation *T = &FET->GetElement2Transformation();
3145 int e = FET->Elem2No;
3146 int e_nbr = e - pmesh.GetNE();
3147 const FiniteElement *fe = dgv_fespace.GetFaceNbrFE(e_nbr);
3148 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
3149 2*order + 2);
3150
3151 double h1_gfc_err = 0.0;
3152 double nd_gfc_err = 0.0;
3153 double rt_gfc_err = 0.0;
3154 double l2_gfc_err = 0.0;
3155 double dgv_gfc_err = 0.0;
3156 double dgi_gfc_err = 0.0;
3157
3158 double h1_gvv_err = 0.0;
3159 double nd_gvv_err = 0.0;
3160 double rt_gvv_err = 0.0;
3161 double l2_gvv_err = 0.0;
3162 double dgv_gvv_err = 0.0;
3163 double dgi_gvv_err = 0.0;
3164
3165 for (int j=0; j<ir.GetNPoints(); j++)
3166 {
3167 npts++;
3168 const IntegrationPoint &ip = ir.IntPoint(j);
3169 T->SetIntPoint(&ip);
3170
3171 funcCoef.Eval(f_val, *T, ip);
3172
3173 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3174 nd_xCoef.Eval(nd_gfc_val, *T, ip);
3175 rt_xCoef.Eval(rt_gfc_val, *T, ip);
3176 l2_xCoef.Eval(l2_gfc_val, *T, ip);
3177 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3178 dgi_xCoef.Eval(dgi_gfc_val, *T, ip);
3179
3180 h1_x.GetVectorValue(e, ip, h1_gvv_val);
3181 nd_x.GetVectorValue(e, ip, nd_gvv_val);
3182 rt_x.GetVectorValue(e, ip, rt_gvv_val);
3183 l2_x.GetVectorValue(e, ip, l2_gvv_val);
3184 dgv_x.GetVectorValue(e, ip, dgv_gvv_val);
3185 dgi_x.GetVectorValue(e, ip, dgi_gvv_val);
3186
3187 double h1_gfc_dist = Distance(f_val, h1_gfc_val, dim);
3188 double nd_gfc_dist = Distance(f_val, nd_gfc_val, dim);
3189 double rt_gfc_dist = Distance(f_val, rt_gfc_val, dim);
3190 double l2_gfc_dist = Distance(f_val, l2_gfc_val, dim);
3191 double dgv_gfc_dist = Distance(f_val, dgv_gfc_val, dim);
3192 double dgi_gfc_dist = Distance(f_val, dgi_gfc_val, dim);
3193
3194 double h1_gvv_dist = Distance(f_val, h1_gvv_val, dim);
3195 double nd_gvv_dist = Distance(f_val, nd_gvv_val, dim);
3196 double rt_gvv_dist = Distance(f_val, rt_gvv_val, dim);
3197 double l2_gvv_dist = Distance(f_val, l2_gvv_val, dim);
3198 double dgv_gvv_dist = Distance(f_val, dgv_gvv_val, dim);
3199 double dgi_gvv_dist = Distance(f_val, dgi_gvv_val, dim);
3200
3201 h1_gfc_err += h1_gfc_dist;
3202 nd_gfc_err += nd_gfc_dist;
3203 rt_gfc_err += rt_gfc_dist;
3204 l2_gfc_err += l2_gfc_dist;
3205 dgv_gfc_err += dgv_gfc_dist;
3206 dgi_gfc_err += dgi_gfc_dist;
3207
3208 h1_gvv_err += h1_gvv_dist;
3209 nd_gvv_err += nd_gvv_dist;
3210 rt_gvv_err += rt_gvv_dist;
3211 l2_gvv_err += l2_gvv_dist;
3212 dgv_gvv_err += dgv_gvv_dist;
3213 dgi_gvv_err += dgi_gvv_dist;
3214
3215 if (log > 0 && h1_gfc_dist > tol)
3216 {
3217 std::cout << e << ":" << j << " h1 gfc ("
3218 << f_val[0] << "," << f_val[1] << ","
3219 << f_val[2] << ") vs. ("
3220 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
3221 << h1_gfc_val[2] << ") "
3222 << h1_gfc_dist << std::endl;
3223 }
3224 if (log > 0 && nd_gfc_dist > tol)
3225 {
3226 std::cout << e << ":" << j << " nd gfc ("
3227 << f_val[0] << "," << f_val[1] << ","
3228 << f_val[2] << ") vs. ("
3229 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ","
3230 << nd_gfc_val[2] << ") "
3231 << nd_gfc_dist << std::endl;
3232 }
3233 if (log > 0 && rt_gfc_dist > tol)
3234 {
3235 std::cout << e << ":" << j << " rt gfc ("
3236 << f_val[0] << "," << f_val[1] << ","
3237 << f_val[2] << ") vs. ("
3238 << rt_gfc_val[0] << "," << rt_gfc_val[1] << ","
3239 << rt_gfc_val[2] << ") "
3240 << rt_gfc_dist << std::endl;
3241 }
3242 if (log > 0 && l2_gfc_dist > tol)
3243 {
3244 std::cout << e << ":" << j << " l2 gfc ("
3245 << f_val[0] << "," << f_val[1] << ","
3246 << f_val[2] << ") vs. ("
3247 << l2_gfc_val[0] << "," << l2_gfc_val[1] << ","
3248 << l2_gfc_val[2] << ") "
3249 << l2_gfc_dist << std::endl;
3250 }
3251 if (log > 0 && dgv_gfc_dist > tol)
3252 {
3253 std::cout << e << ":" << j << " dgv gfc ("
3254 << f_val[0] << "," << f_val[1] << ","
3255 << f_val[2] << ") vs. ("
3256 << dgv_gfc_val[0] << ","
3257 << dgv_gfc_val[1] << ","
3258 << dgv_gfc_val[2] << ") "
3259 << dgv_gfc_dist << std::endl;
3260 }
3261 if (log > 0 && dgi_gfc_dist > tol)
3262 {
3263 std::cout << e << ":" << j << " dgi gfc ("
3264 << f_val[0] << "," << f_val[1] << ","
3265 << f_val[2] << ") vs. ("
3266 << dgi_gfc_val[0] << ","
3267 << dgi_gfc_val[1] << ","
3268 << dgi_gfc_val[2] << ") "
3269 << dgi_gfc_dist << std::endl;
3270 }
3271 if (log > 0 && h1_gvv_dist > tol)
3272 {
3273 std::cout << e << ":" << j << " h1 gvv ("
3274 << f_val[0] << "," << f_val[1] << ","
3275 << f_val[2] << ") vs. ("
3276 << h1_gvv_val[0] << "," << h1_gvv_val[1] << ","
3277 << h1_gvv_val[2] << ") "
3278 << h1_gvv_dist << std::endl;
3279 }
3280 if (log > 0 && nd_gvv_dist > tol)
3281 {
3282 std::cout << e << ":" << j << " nd gvv ("
3283 << f_val[0] << "," << f_val[1] << ","
3284 << f_val[2] << ") vs. ("
3285 << nd_gvv_val[0] << "," << nd_gvv_val[1] << ","
3286 << nd_gvv_val[2] << ") "
3287 << nd_gvv_dist << std::endl;
3288 }
3289 if (log > 0 && rt_gvv_dist > tol)
3290 {
3291 std::cout << e << ":" << j << " rt gvv ("
3292 << f_val[0] << "," << f_val[1] << ","
3293 << f_val[2] << ") vs. ("
3294 << rt_gvv_val[0] << "," << rt_gvv_val[1] << ","
3295 << rt_gvv_val[2] << ") "
3296 << rt_gvv_dist << std::endl;
3297 }
3298 if (log > 0 && l2_gvv_dist > tol)
3299 {
3300 std::cout << e << ":" << j << " l2 gvv ("
3301 << f_val[0] << "," << f_val[1] << ","
3302 << f_val[2] << ") vs. ("
3303 << l2_gvv_val[0] << "," << l2_gvv_val[1] << ","
3304 << l2_gvv_val[2] << ") "
3305 << l2_gvv_dist << std::endl;
3306 }
3307 if (log > 0 && dgv_gvv_dist > tol)
3308 {
3309 std::cout << e << ":" << j << " dgv gvv ("
3310 << f_val[0] << "," << f_val[1] << ","
3311 << f_val[2] << ") vs. ("
3312 << dgv_gvv_val[0] << ","
3313 << dgv_gvv_val[1] << ","
3314 << dgv_gvv_val[2] << ") "
3315 << dgv_gvv_dist << std::endl;
3316 }
3317 if (log > 0 && dgi_gvv_dist > tol)
3318 {
3319 std::cout << e << ":" << j << " dgi gvv ("
3320 << f_val[0] << "," << f_val[1] << ","
3321 << f_val[2] << ") vs. ("
3322 << dgi_gvv_val[0] << ","
3323 << dgi_gvv_val[1] << ","
3324 << dgi_gvv_val[2] << ") "
3325 << dgi_gvv_dist << std::endl;
3326 }
3327 }
3328
3329 h1_gfc_err /= ir.GetNPoints();
3330 nd_gfc_err /= ir.GetNPoints();
3331 rt_gfc_err /= ir.GetNPoints();
3332 l2_gfc_err /= ir.GetNPoints();
3333 dgv_gfc_err /= ir.GetNPoints();
3334 dgi_gfc_err /= ir.GetNPoints();
3335
3336 h1_gvv_err /= ir.GetNPoints();
3337 nd_gvv_err /= ir.GetNPoints();
3338 rt_gvv_err /= ir.GetNPoints();
3339 l2_gvv_err /= ir.GetNPoints();
3340 dgv_gvv_err /= ir.GetNPoints();
3341 dgi_gvv_err /= ir.GetNPoints();
3342
3343 REQUIRE( h1_gfc_err == MFEM_Approx(0.0));
3344 REQUIRE( nd_gfc_err == MFEM_Approx(0.0));
3345 REQUIRE( rt_gfc_err == MFEM_Approx(0.0));
3346 REQUIRE( l2_gfc_err == MFEM_Approx(0.0));
3347 REQUIRE(dgv_gfc_err == MFEM_Approx(0.0));
3348 REQUIRE(dgi_gfc_err == MFEM_Approx(0.0));
3349
3350 REQUIRE( h1_gvv_err == MFEM_Approx(0.0));
3351 REQUIRE( nd_gvv_err == MFEM_Approx(0.0));
3352 REQUIRE( rt_gvv_err == MFEM_Approx(0.0));
3353 REQUIRE( l2_gvv_err == MFEM_Approx(0.0));
3354 REQUIRE(dgv_gvv_err == MFEM_Approx(0.0));
3355 REQUIRE(dgi_gvv_err == MFEM_Approx(0.0));
3356 }
3357 }
3358 }
3359 }
3360 std::cout << my_rank << ": Checked GridFunction::GetVectorValue at "
3361 << npts << " 3D points" << std::endl;
3362 }
3363
3364 #endif // MFEM_USE_MPI
3365
3366 TEST_CASE("1D GetGradient",
3367 "[GridFunction]"
3368 "[GradientGridFunctionCoefficient]")
3369 {
3370 int log = 1;
3371 int n = 1;
3372 int dim = 1;
3373 int order = 2;
3374 int npts = 0;
3375
3376 double tol = 1e-6;
3377
3378 for (int type = (int)Element::SEGMENT;
3379 type <= (int)Element::SEGMENT; type++)
3380 {
3381 Mesh mesh = Mesh::MakeCartesian1D(n, 2.0);
3382
3383 FunctionCoefficient funcCoef(func_1D_quad);
3384 VectorFunctionCoefficient dFuncCoef(dim, dfunc_1D_quad);
3385
to_string(type)3386 SECTION("1D GetGradient tests for element type " + std::to_string(type))
3387 {
3388 H1_FECollection h1_fec(order, dim);
3389 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
3390 FiniteElement::VALUE);
3391
3392 FiniteElementSpace h1_fespace(&mesh, &h1_fec);
3393 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
3394
3395 GridFunction h1_x(&h1_fespace);
3396 GridFunction dgv_x(&dgv_fespace);
3397
3398 GradientGridFunctionCoefficient h1_xCoef(&h1_x);
3399 GradientGridFunctionCoefficient dgv_xCoef(&dgv_x);
3400
3401 h1_x.ProjectCoefficient(funcCoef);
3402 dgv_x.ProjectCoefficient(funcCoef);
3403
3404 Vector f_val(dim); f_val = 0.0;
3405 Vector h1_gfc_val(dim); h1_gfc_val = 0.0;
3406 Vector dgv_gfc_val(dim); dgv_gfc_val = 0.0;
3407
3408 SECTION("Domain Evaluation 1D")
3409 {
3410 std::cout << "Domain Evaluation 1D" << std::endl;
3411 for (int e = 0; e < mesh.GetNE(); e++)
3412 {
3413 ElementTransformation *T = mesh.GetElementTransformation(e);
3414 const FiniteElement *fe = h1_fespace.GetFE(e);
3415 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
3416 2*order + 2);
3417
3418 double h1_err = 0.0;
3419 double dgv_err = 0.0;
3420
3421 for (int j=0; j<ir.GetNPoints(); j++)
3422 {
3423 npts++;
3424 const IntegrationPoint &ip = ir.IntPoint(j);
3425 T->SetIntPoint(&ip);
3426
3427 dFuncCoef.Eval(f_val, *T, ip);
3428 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3429 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3430
3431 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3432 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3433
3434 h1_err += h1_dist;
3435 dgv_err += dgv_dist;
3436
3437 if (log > 0 && h1_dist > tol)
3438 {
3439 std::cout << e << ":" << j << " h1 ("
3440 << f_val[0] << ") vs. ("
3441 << h1_gfc_val[0] << ") "
3442 << h1_dist << std::endl;
3443 }
3444 if (log > 0 && dgv_dist > tol)
3445 {
3446 std::cout << e << ":" << j << " dgv ("
3447 << f_val[0] << ") vs. ("
3448 << dgv_gfc_val[0] << ") "
3449 << dgv_dist << std::endl;
3450 }
3451 }
3452 h1_err /= ir.GetNPoints();
3453 dgv_err /= ir.GetNPoints();
3454
3455 REQUIRE(h1_err == MFEM_Approx(0.0));
3456 REQUIRE(dgv_err == MFEM_Approx(0.0));
3457 }
3458 }
3459
3460 SECTION("Boundary Evaluation 1D (H1 Context)")
3461 {
3462 std::cout << "Boundary Evaluation 1D (H1 Context)" << std::endl;
3463 for (int be = 0; be < mesh.GetNBE(); be++)
3464 {
3465 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
3466 const FiniteElement *fe = h1_fespace.GetBE(be);
3467 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
3468 2*order + 2);
3469
3470 double h1_err = 0.0;
3471 double dgv_err = 0.0;
3472
3473 for (int j=0; j<ir.GetNPoints(); j++)
3474 {
3475 npts++;
3476 const IntegrationPoint &ip = ir.IntPoint(j);
3477 T->SetIntPoint(&ip);
3478
3479 dFuncCoef.Eval(f_val, *T, ip);
3480 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3481 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3482
3483 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3484 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3485
3486 h1_err += h1_dist;
3487 dgv_err += dgv_dist;
3488
3489 if (log > 0 && h1_dist > tol)
3490 {
3491 std::cout << be << ":" << j << " h1 ("
3492 << f_val[0] << ") vs. ("
3493 << h1_gfc_val[0] << ") "
3494 << h1_dist << std::endl;
3495 }
3496 if (log > 0 && dgv_dist > tol)
3497 {
3498 std::cout << be << ":" << j << " dgv ("
3499 << f_val[0] << ") vs. ("
3500 << dgv_gfc_val[0] << ") "
3501 << dgv_dist << std::endl;
3502 }
3503 }
3504 h1_err /= ir.GetNPoints();
3505 dgv_err /= ir.GetNPoints();
3506
3507 REQUIRE(h1_err == MFEM_Approx(0.0));
3508 REQUIRE(dgv_err == MFEM_Approx(0.0));
3509 }
3510 }
3511
3512 SECTION("Boundary Evaluation 1D (DG Context)")
3513 {
3514 std::cout << "Boundary Evaluation 1D (DG Context)" << std::endl;
3515 for (int be = 0; be < mesh.GetNBE(); be++)
3516 {
3517 FaceElementTransformations *T =
3518 mesh.GetBdrFaceTransformations(be);
3519 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
3520 2*order + 2);
3521
3522 double h1_err = 0.0;
3523 double dgv_err = 0.0;
3524
3525 for (int j=0; j<ir.GetNPoints(); j++)
3526 {
3527 npts++;
3528 const IntegrationPoint &ip = ir.IntPoint(j);
3529
3530 T->SetIntPoint(&ip);
3531
3532 dFuncCoef.Eval(f_val, *T, ip);
3533 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3534 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3535
3536 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3537 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3538
3539 h1_err += h1_dist;
3540 dgv_err += dgv_dist;
3541
3542 if (log > 0 && h1_dist > tol)
3543 {
3544 std::cout << be << ":" << j << " h1 ("
3545 << f_val[0] << ") vs. ("
3546 << h1_gfc_val[0] << ") "
3547 << h1_dist << std::endl;
3548 }
3549 if (log > 0 && dgv_dist > tol)
3550 {
3551 std::cout << be << ":" << j << " dgv ("
3552 << f_val[0] << ") vs. ("
3553 << dgv_gfc_val[0] << ") "
3554 << dgv_dist << std::endl;
3555 }
3556 }
3557 h1_err /= ir.GetNPoints();
3558 dgv_err /= ir.GetNPoints();
3559
3560 REQUIRE(h1_err == MFEM_Approx(0.0));
3561 REQUIRE(dgv_err == MFEM_Approx(0.0));
3562 }
3563 }
3564 }
3565 }
3566 std::cout << "Checked GridFunction::GetGradient at "
3567 << npts << " 1D points" << std::endl;
3568 }
3569
3570 TEST_CASE("2D GetGradient",
3571 "[GridFunction]"
3572 "[GradientGridFunctionCoefficient]")
3573 {
3574 int log = 1;
3575 int n = 1;
3576 int dim = 2;
3577 int order = 2;
3578 int npts = 0;
3579
3580 double tol = 1e-6;
3581
3582 for (int type = (int)Element::TRIANGLE;
3583 type <= (int)Element::QUADRILATERAL; type++)
3584 {
3585 Mesh mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
3586
3587 FunctionCoefficient funcCoef(func_2D_quad);
3588 VectorFunctionCoefficient dFuncCoef(dim, dfunc_2D_quad);
3589
to_string(type)3590 SECTION("2D GetGradient tests for element type " + std::to_string(type))
3591 {
3592 H1_FECollection h1_fec(order, dim);
3593 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
3594 FiniteElement::VALUE);
3595
3596 FiniteElementSpace h1_fespace(&mesh, &h1_fec);
3597 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
3598
3599 GridFunction h1_x(&h1_fespace);
3600 GridFunction dgv_x(&dgv_fespace);
3601
3602 GradientGridFunctionCoefficient h1_xCoef(&h1_x);
3603 GradientGridFunctionCoefficient dgv_xCoef(&dgv_x);
3604
3605 h1_x.ProjectCoefficient(funcCoef);
3606 dgv_x.ProjectCoefficient(funcCoef);
3607
3608 Vector f_val(dim); f_val = 0.0;
3609 Vector h1_gfc_val(dim); h1_gfc_val = 0.0;
3610 Vector dgv_gfc_val(dim); dgv_gfc_val = 0.0;
3611
3612 SECTION("Domain Evaluation 2D")
3613 {
3614 std::cout << "Domain Evaluation 2D" << std::endl;
3615 for (int e = 0; e < mesh.GetNE(); e++)
3616 {
3617 ElementTransformation *T = mesh.GetElementTransformation(e);
3618 const FiniteElement *fe = h1_fespace.GetFE(e);
3619 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
3620 2*order + 2);
3621
3622 double h1_err = 0.0;
3623 double dgv_err = 0.0;
3624
3625 for (int j=0; j<ir.GetNPoints(); j++)
3626 {
3627 npts++;
3628 const IntegrationPoint &ip = ir.IntPoint(j);
3629 T->SetIntPoint(&ip);
3630
3631 dFuncCoef.Eval(f_val, *T, ip);
3632 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3633 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3634
3635 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3636 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3637
3638 h1_err += h1_dist;
3639 dgv_err += dgv_dist;
3640
3641 if (log > 0 && h1_dist > tol)
3642 {
3643 std::cout << e << ":" << j << " h1 ("
3644 << f_val[0] << "," << f_val[1] << ") vs. ("
3645 << h1_gfc_val[0] << ","
3646 << h1_gfc_val[1] << ") "
3647 << h1_dist << std::endl;
3648 }
3649 if (log > 0 && dgv_dist > tol)
3650 {
3651 std::cout << e << ":" << j << " dgv ("
3652 << f_val[0] << "," << f_val[1] << ") vs. ("
3653 << dgv_gfc_val[0] << ","
3654 << dgv_gfc_val[1] << ") "
3655 << dgv_dist << std::endl;
3656 }
3657 }
3658 h1_err /= ir.GetNPoints();
3659 dgv_err /= ir.GetNPoints();
3660
3661 REQUIRE(h1_err == MFEM_Approx(0.0));
3662 REQUIRE(dgv_err == MFEM_Approx(0.0));
3663 }
3664 }
3665
3666 SECTION("Boundary Evaluation 2D (H1 Context)")
3667 {
3668 std::cout << "Boundary Evaluation 2D (H1 Context)" << std::endl;
3669 for (int be = 0; be < mesh.GetNBE(); be++)
3670 {
3671 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
3672 const FiniteElement *fe = h1_fespace.GetBE(be);
3673 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
3674 2*order + 2);
3675
3676 double h1_err = 0.0;
3677 double dgv_err = 0.0;
3678
3679 for (int j=0; j<ir.GetNPoints(); j++)
3680 {
3681 npts++;
3682 const IntegrationPoint &ip = ir.IntPoint(j);
3683 T->SetIntPoint(&ip);
3684
3685 dFuncCoef.Eval(f_val, *T, ip);
3686 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3687 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3688
3689 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3690 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3691
3692 h1_err += h1_dist;
3693 dgv_err += dgv_dist;
3694
3695 if (log > 0 && h1_dist > tol)
3696 {
3697 std::cout << be << ":" << j << " h1 ("
3698 << f_val[0] << "," << f_val[1] << ") vs. ("
3699 << h1_gfc_val[0] << ","
3700 << h1_gfc_val[1] << ") "
3701 << h1_dist << std::endl;
3702 }
3703 if (log > 0 && dgv_dist > tol)
3704 {
3705 std::cout << be << ":" << j << " dgv ("
3706 << f_val[0] << "," << f_val[1] << ") vs. ("
3707 << dgv_gfc_val[0] << ","
3708 << dgv_gfc_val[1] << ") "
3709 << dgv_dist << std::endl;
3710 }
3711 }
3712 h1_err /= ir.GetNPoints();
3713 dgv_err /= ir.GetNPoints();
3714
3715 REQUIRE(h1_err == MFEM_Approx(0.0));
3716 REQUIRE(dgv_err == MFEM_Approx(0.0));
3717 }
3718 }
3719
3720 SECTION("Boundary Evaluation 2D (DG Context)")
3721 {
3722 std::cout << "Boundary Evaluation 2D (DG Context)" << std::endl;
3723 for (int be = 0; be < mesh.GetNBE(); be++)
3724 {
3725 FaceElementTransformations *T =
3726 mesh.GetBdrFaceTransformations(be);
3727 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
3728 2*order + 2);
3729
3730 double h1_err = 0.0;
3731 double dgv_err = 0.0;
3732
3733 for (int j=0; j<ir.GetNPoints(); j++)
3734 {
3735 npts++;
3736 const IntegrationPoint &ip = ir.IntPoint(j);
3737
3738 T->SetIntPoint(&ip);
3739
3740 dFuncCoef.Eval(f_val, *T, ip);
3741 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3742 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3743
3744 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3745 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3746
3747 h1_err += h1_dist;
3748 dgv_err += dgv_dist;
3749
3750 if (log > 0 && h1_dist > tol)
3751 {
3752 std::cout << be << ":" << j << " h1 ("
3753 << f_val[0] << "," << f_val[1] << ") vs. ("
3754 << h1_gfc_val[0] << ","
3755 << h1_gfc_val[1] << ") "
3756 << h1_dist << std::endl;
3757 }
3758 if (log > 0 && dgv_dist > tol)
3759 {
3760 std::cout << be << ":" << j << " dgv ("
3761 << f_val[0] << "," << f_val[1] << ") vs. ("
3762 << dgv_gfc_val[0] << ","
3763 << dgv_gfc_val[1] << ") "
3764 << dgv_dist << std::endl;
3765 }
3766 }
3767 h1_err /= ir.GetNPoints();
3768 dgv_err /= ir.GetNPoints();
3769
3770 REQUIRE(h1_err == MFEM_Approx(0.0));
3771 REQUIRE(dgv_err == MFEM_Approx(0.0));
3772 }
3773 }
3774 }
3775 }
3776 std::cout << "Checked GridFunction::GetGradient at "
3777 << npts << " 2D points" << std::endl;
3778 }
3779
3780 TEST_CASE("3D GetGradient",
3781 "[GridFunction]"
3782 "[GradientGridFunctionCoefficient]")
3783 {
3784 int log = 1;
3785 int n = 1;
3786 int dim = 3;
3787 int order = 2;
3788 int npts = 0;
3789
3790 double tol = 1e-6;
3791
3792 for (int type = (int)Element::TETRAHEDRON;
3793 type <= (int)Element::WEDGE; type++)
3794 {
3795 Mesh mesh = Mesh::MakeCartesian3D(
3796 n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
3797
3798 FunctionCoefficient funcCoef(func_3D_quad);
3799 VectorFunctionCoefficient dFuncCoef(dim, dfunc_3D_quad);
3800
to_string(type)3801 SECTION("3D GetGradient tests for element type " + std::to_string(type))
3802 {
3803 H1_FECollection h1_fec(order, dim);
3804 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
3805 FiniteElement::VALUE);
3806
3807 FiniteElementSpace h1_fespace(&mesh, &h1_fec);
3808 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec);
3809
3810 GridFunction h1_x(&h1_fespace);
3811 GridFunction dgv_x(&dgv_fespace);
3812
3813 GradientGridFunctionCoefficient h1_xCoef(&h1_x);
3814 GradientGridFunctionCoefficient dgv_xCoef(&dgv_x);
3815
3816 h1_x.ProjectCoefficient(funcCoef);
3817 dgv_x.ProjectCoefficient(funcCoef);
3818
3819 Vector f_val(dim); f_val = 0.0;
3820 Vector h1_gfc_val(dim); h1_gfc_val = 0.0;
3821 Vector dgv_gfc_val(dim); dgv_gfc_val = 0.0;
3822
3823 SECTION("Domain Evaluation 3D")
3824 {
3825 std::cout << "Domain Evaluation 3D" << std::endl;
3826 for (int e = 0; e < mesh.GetNE(); e++)
3827 {
3828 ElementTransformation *T = mesh.GetElementTransformation(e);
3829 const FiniteElement *fe = h1_fespace.GetFE(e);
3830 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
3831 2*order + 2);
3832
3833 double h1_err = 0.0;
3834 double dgv_err = 0.0;
3835
3836 for (int j=0; j<ir.GetNPoints(); j++)
3837 {
3838 npts++;
3839 const IntegrationPoint &ip = ir.IntPoint(j);
3840 T->SetIntPoint(&ip);
3841
3842 dFuncCoef.Eval(f_val, *T, ip);
3843 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3844 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3845
3846 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3847 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3848
3849 h1_err += h1_dist;
3850 dgv_err += dgv_dist;
3851
3852 if (log > 0 && h1_dist > tol)
3853 {
3854 std::cout << e << ":" << j << " h1 ("
3855 << f_val[0] << "," << f_val[1] << ","
3856 << f_val[2] << ") vs. ("
3857 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
3858 << h1_gfc_val[2] << ") "
3859 << h1_dist << std::endl;
3860 }
3861 if (log > 0 && dgv_dist > tol)
3862 {
3863 std::cout << e << ":" << j << " dgv ("
3864 << f_val[0] << "," << f_val[1] << ","
3865 << f_val[2] << ") vs. ("
3866 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
3867 << dgv_gfc_val[2] << ") " << dgv_dist
3868 << std::endl;
3869 }
3870 }
3871 h1_err /= ir.GetNPoints();
3872 dgv_err /= ir.GetNPoints();
3873
3874 REQUIRE(h1_err == MFEM_Approx(0.0));
3875 REQUIRE(dgv_err == MFEM_Approx(0.0));
3876 }
3877 }
3878
3879 SECTION("Boundary Evaluation 3D (H1 Context)")
3880 {
3881 std::cout << "Boundary Evaluation 3D (H1 Context)" << std::endl;
3882 for (int be = 0; be < mesh.GetNBE(); be++)
3883 {
3884 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
3885 const FiniteElement *fe = h1_fespace.GetBE(be);
3886 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
3887 2*order + 2);
3888
3889 double h1_err = 0.0;
3890 double dgv_err = 0.0;
3891
3892 for (int j=0; j<ir.GetNPoints(); j++)
3893 {
3894 npts++;
3895 const IntegrationPoint &ip = ir.IntPoint(j);
3896 T->SetIntPoint(&ip);
3897
3898 dFuncCoef.Eval(f_val, *T, ip);
3899 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3900 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3901
3902 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3903 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3904
3905 h1_err += h1_dist;
3906 dgv_err += dgv_dist;
3907
3908 if (log > 0 && h1_dist > tol)
3909 {
3910 std::cout << be << ":" << j << " h1 ("
3911 << f_val[0] << "," << f_val[1] << ","
3912 << f_val[2] << ") vs. ("
3913 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
3914 << h1_gfc_val[2] << ") "
3915 << h1_dist << std::endl;
3916 }
3917 if (log > 0 && dgv_dist > tol)
3918 {
3919 std::cout << be << ":" << j << " dgv ("
3920 << f_val[0] << "," << f_val[1] << ","
3921 << f_val[2] << ") vs. ("
3922 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
3923 << dgv_gfc_val[2] << ") " << dgv_dist
3924 << std::endl;
3925 }
3926 }
3927 h1_err /= ir.GetNPoints();
3928 dgv_err /= ir.GetNPoints();
3929
3930 REQUIRE(h1_err == MFEM_Approx(0.0));
3931 REQUIRE(dgv_err == MFEM_Approx(0.0));
3932 }
3933 }
3934
3935 SECTION("Boundary Evaluation 3D (DG Context)")
3936 {
3937 std::cout << "Boundary Evaluation 3D (DG Context)" << std::endl;
3938 for (int be = 0; be < mesh.GetNBE(); be++)
3939 {
3940 FaceElementTransformations *T =
3941 mesh.GetBdrFaceTransformations(be);
3942 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
3943 2*order + 2);
3944
3945 double h1_err = 0.0;
3946 double dgv_err = 0.0;
3947
3948 for (int j=0; j<ir.GetNPoints(); j++)
3949 {
3950 npts++;
3951 const IntegrationPoint &ip = ir.IntPoint(j);
3952
3953 T->SetIntPoint(&ip);
3954
3955 dFuncCoef.Eval(f_val, *T, ip);
3956 h1_xCoef.Eval(h1_gfc_val, *T, ip);
3957 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
3958
3959 double h1_dist = Distance(f_val, h1_gfc_val, dim);
3960 double dgv_dist = Distance(f_val, dgv_gfc_val, dim);
3961
3962 h1_err += h1_dist;
3963 dgv_err += dgv_dist;
3964
3965 if (log > 0 && h1_dist > tol)
3966 {
3967 std::cout << be << ":" << j << " h1 ("
3968 << f_val[0] << "," << f_val[1] << ","
3969 << f_val[2] << ") vs. ("
3970 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
3971 << h1_gfc_val[2] << ") "
3972 << h1_dist << std::endl;
3973 }
3974 if (log > 0 && dgv_dist > tol)
3975 {
3976 std::cout << be << ":" << j << " dgv ("
3977 << f_val[0] << "," << f_val[1] << ","
3978 << f_val[2] << ") vs. ("
3979 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
3980 << dgv_gfc_val[2] << ") " << dgv_dist
3981 << std::endl;
3982 }
3983 }
3984 h1_err /= ir.GetNPoints();
3985 dgv_err /= ir.GetNPoints();
3986
3987 REQUIRE(h1_err == MFEM_Approx(0.0));
3988 REQUIRE(dgv_err == MFEM_Approx(0.0));
3989 }
3990 }
3991 }
3992 }
3993 std::cout << "Checked GridFunction::GetGradient at "
3994 << npts << " 3D points" << std::endl;
3995 }
3996
3997 TEST_CASE("2D GetCurl",
3998 "[GridFunction]"
3999 "[CurlGridFunctionCoefficient]")
4000 {
4001 int log = 1;
4002 int n = 1;
4003 int dim = 2;
4004 int order = 2;
4005 int npts = 0;
4006
4007 double tol = 1e-6;
4008
4009 for (int type = (int)Element::TRIANGLE;
4010 type <= (int)Element::QUADRILATERAL; type++)
4011 {
4012 Mesh mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
4013
4014 VectorFunctionCoefficient funcCoef(dim, Func_2D_quad);
4015 VectorFunctionCoefficient dFuncCoef(1, RotFunc_2D_quad);
4016
4017 SECTION("2D GetCurl tests for element type " +
to_string(type)4018 std::to_string(type))
4019 {
4020 H1_FECollection h1_fec(order, dim);
4021 ND_FECollection nd_fec(order+1, dim);
4022 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
4023 FiniteElement::VALUE);
4024
4025 FiniteElementSpace h1_fespace(&mesh, &h1_fec, dim);
4026 FiniteElementSpace nd_fespace(&mesh, &nd_fec);
4027 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec, dim);
4028
4029 GridFunction h1_x( &h1_fespace);
4030 GridFunction nd_x( &nd_fespace);
4031 GridFunction dgv_x(&dgv_fespace);
4032
4033 CurlGridFunctionCoefficient h1_xCoef( &h1_x);
4034 CurlGridFunctionCoefficient nd_xCoef( &nd_x);
4035 CurlGridFunctionCoefficient dgv_xCoef(&dgv_x);
4036
4037 h1_x.ProjectCoefficient(funcCoef);
4038 nd_x.ProjectCoefficient(funcCoef);
4039 dgv_x.ProjectCoefficient(funcCoef);
4040
4041 Vector f_val(2*dim-3); f_val = 0.0;
4042 Vector h1_gfc_val(2*dim-3); h1_gfc_val = 0.0;
4043 Vector nd_gfc_val(2*dim-3); nd_gfc_val = 0.0;
4044 Vector dgv_gfc_val(2*dim-3); dgv_gfc_val = 0.0;
4045
4046 SECTION("Domain Evaluation 2D")
4047 {
4048 std::cout << "Domain Evaluation 2D" << std::endl;
4049 for (int e = 0; e < mesh.GetNE(); e++)
4050 {
4051 ElementTransformation *T = mesh.GetElementTransformation(e);
4052 const FiniteElement *fe = h1_fespace.GetFE(e);
4053 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4054 2*order + 2);
4055
4056 double h1_err = 0.0;
4057 double nd_err = 0.0;
4058 double dgv_err = 0.0;
4059
4060 for (int j=0; j<ir.GetNPoints(); j++)
4061 {
4062 npts++;
4063 const IntegrationPoint &ip = ir.IntPoint(j);
4064 T->SetIntPoint(&ip);
4065
4066 dFuncCoef.Eval(f_val, *T, ip);
4067 h1_xCoef.Eval(h1_gfc_val, *T, ip);
4068 nd_xCoef.Eval(nd_gfc_val, *T, ip);
4069 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
4070
4071 double h1_dist = Distance(f_val, h1_gfc_val, 2*dim-3);
4072 double nd_dist = Distance(f_val, nd_gfc_val, 2*dim-3);
4073 double dgv_dist = Distance(f_val, dgv_gfc_val, 2*dim-3);
4074
4075 h1_err += h1_dist;
4076 nd_err += nd_dist;
4077 dgv_err += dgv_dist;
4078
4079 if (log > 0 && h1_dist > tol)
4080 {
4081 std::cout << e << ":" << j << " h1 ("
4082 << f_val[0] << ") vs. ("
4083 << h1_gfc_val[0] << ") " << h1_dist
4084 << std::endl;
4085 }
4086 if (log > 0 && nd_dist > tol)
4087 {
4088 std::cout << e << ":" << j << " nd ("
4089 << f_val[0] << ") vs. ("
4090 << nd_gfc_val[0] << ") " << nd_dist
4091 << std::endl;
4092 }
4093 if (log > 0 && dgv_dist > tol)
4094 {
4095 std::cout << e << ":" << j << " dgv ("
4096 << f_val[0] << ") vs. ("
4097 << dgv_gfc_val[0] << ") " << dgv_dist
4098 << std::endl;
4099 }
4100 }
4101 h1_err /= ir.GetNPoints();
4102 nd_err /= ir.GetNPoints();
4103 dgv_err /= ir.GetNPoints();
4104
4105 REQUIRE( h1_err == MFEM_Approx(0.0));
4106 REQUIRE( nd_err == MFEM_Approx(0.0));
4107 REQUIRE(dgv_err == MFEM_Approx(0.0));
4108 }
4109 }
4110
4111 SECTION("Boundary Evaluation 2D (H1 Context)")
4112 {
4113 std::cout << "Boundary Evaluation 2D (H1 Context)" << std::endl;
4114 for (int be = 0; be < mesh.GetNBE(); be++)
4115 {
4116 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
4117 const FiniteElement *fe = h1_fespace.GetBE(be);
4118 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4119 2*order + 2);
4120
4121 double h1_err = 0.0;
4122 double nd_err = 0.0;
4123 double dgv_err = 0.0;
4124
4125 for (int j=0; j<ir.GetNPoints(); j++)
4126 {
4127 npts++;
4128 const IntegrationPoint &ip = ir.IntPoint(j);
4129 T->SetIntPoint(&ip);
4130
4131 dFuncCoef.Eval(f_val, *T, ip);
4132 h1_xCoef.Eval(h1_gfc_val, *T, ip);
4133 nd_xCoef.Eval(nd_gfc_val, *T, ip);
4134 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
4135
4136 double h1_dist = Distance(f_val, h1_gfc_val, 2*dim-3);
4137 double nd_dist = Distance(f_val, nd_gfc_val, 2*dim-3);
4138 double dgv_dist = Distance(f_val, dgv_gfc_val, 2*dim-3);
4139
4140 h1_err += h1_dist;
4141 nd_err += nd_dist;
4142 dgv_err += dgv_dist;
4143
4144 if (log > 0 && h1_dist > tol)
4145 {
4146 std::cout << be << ":" << j << " h1 ("
4147 << f_val[0] << ") vs. ("
4148 << h1_gfc_val[0] << ") " << h1_dist
4149 << std::endl;
4150 }
4151 if (log > 0 && nd_dist > tol)
4152 {
4153 std::cout << be << ":" << j << " nd ("
4154 << f_val[0] << ") vs. ("
4155 << nd_gfc_val[0] << ") " << nd_dist
4156 << std::endl;
4157 }
4158 if (log > 0 && dgv_dist > tol)
4159 {
4160 std::cout << be << ":" << j << " dgv ("
4161 << f_val[0] << ") vs. ("
4162 << dgv_gfc_val[0] << ") " << dgv_dist
4163 << std::endl;
4164 }
4165 }
4166 h1_err /= ir.GetNPoints();
4167 nd_err /= ir.GetNPoints();
4168 dgv_err /= ir.GetNPoints();
4169
4170 REQUIRE( h1_err == MFEM_Approx(0.0));
4171 REQUIRE( nd_err == MFEM_Approx(0.0));
4172 REQUIRE(dgv_err == MFEM_Approx(0.0));
4173 }
4174 }
4175
4176 SECTION("Boundary Evaluation 2D (DG Context)")
4177 {
4178 std::cout << "Boundary Evaluation 2D (DG Context)" << std::endl;
4179 for (int be = 0; be < mesh.GetNBE(); be++)
4180 {
4181 FaceElementTransformations *T =
4182 mesh.GetBdrFaceTransformations(be);
4183 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
4184 2*order + 2);
4185
4186 double h1_err = 0.0;
4187 double nd_err = 0.0;
4188 double dgv_err = 0.0;
4189
4190 for (int j=0; j<ir.GetNPoints(); j++)
4191 {
4192 npts++;
4193 const IntegrationPoint &ip = ir.IntPoint(j);
4194
4195 T->SetIntPoint(&ip);
4196
4197 dFuncCoef.Eval(f_val, *T, ip);
4198 h1_xCoef.Eval(h1_gfc_val, *T, ip);
4199 nd_xCoef.Eval(nd_gfc_val, *T, ip);
4200 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
4201
4202 double h1_dist = Distance(f_val, h1_gfc_val, 2*dim-3);
4203 double nd_dist = Distance(f_val, nd_gfc_val, 2*dim-3);
4204 double dgv_dist = Distance(f_val, dgv_gfc_val, 2*dim-3);
4205
4206 h1_err += h1_dist;
4207 nd_err += nd_dist;
4208 dgv_err += dgv_dist;
4209
4210 if (log > 0 && h1_dist > tol)
4211 {
4212 std::cout << be << ":" << j << " h1 ("
4213 << f_val[0] << ") vs. ("
4214 << h1_gfc_val[0] << ") " << h1_dist
4215 << std::endl;
4216 }
4217 if (log > 0 && nd_dist > tol)
4218 {
4219 std::cout << be << ":" << j << " nd ("
4220 << f_val[0] << ") vs. ("
4221 << nd_gfc_val[0] << ") " << nd_dist
4222 << std::endl;
4223 }
4224 if (log > 0 && dgv_dist > tol)
4225 {
4226 std::cout << be << ":" << j << " dgv ("
4227 << f_val[0] << ") vs. ("
4228 << dgv_gfc_val[0] << ") " << dgv_dist
4229 << std::endl;
4230 }
4231 }
4232 h1_err /= ir.GetNPoints();
4233 nd_err /= ir.GetNPoints();
4234 dgv_err /= ir.GetNPoints();
4235
4236 REQUIRE( h1_err == MFEM_Approx(0.0));
4237 REQUIRE( nd_err == MFEM_Approx(0.0));
4238 REQUIRE(dgv_err == MFEM_Approx(0.0));
4239 }
4240 }
4241 }
4242 }
4243 std::cout << "Checked GridFunction::GetCurl at "
4244 << npts << " 2D points" << std::endl;
4245 }
4246
4247 TEST_CASE("3D GetCurl",
4248 "[GridFunction]"
4249 "[CurlGridFunctionCoefficient]")
4250 {
4251 int log = 1;
4252 int n = 1;
4253 int dim = 3;
4254 int order = 2;
4255 int npts = 0;
4256
4257 double tol = 1e-6;
4258
4259 for (int type = (int)Element::TETRAHEDRON;
4260 type <= (int)Element::HEXAHEDRON; type++)
4261 {
4262 Mesh mesh = Mesh::MakeCartesian3D(
4263 n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
4264
4265 VectorFunctionCoefficient funcCoef(dim, Func_3D_quad);
4266 VectorFunctionCoefficient dFuncCoef(dim, CurlFunc_3D_quad);
4267
4268 SECTION("3D GetCurl tests for element type " +
to_string(type)4269 std::to_string(type))
4270 {
4271 H1_FECollection h1_fec(order, dim);
4272 ND_FECollection nd_fec(order+1, dim);
4273 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
4274 FiniteElement::VALUE);
4275
4276 FiniteElementSpace h1_fespace(&mesh, &h1_fec, dim);
4277 FiniteElementSpace nd_fespace(&mesh, &nd_fec);
4278 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec, dim);
4279
4280 GridFunction h1_x( &h1_fespace);
4281 GridFunction nd_x( &nd_fespace);
4282 GridFunction dgv_x(&dgv_fespace);
4283
4284 CurlGridFunctionCoefficient h1_xCoef( &h1_x);
4285 CurlGridFunctionCoefficient nd_xCoef( &nd_x);
4286 CurlGridFunctionCoefficient dgv_xCoef(&dgv_x);
4287
4288 h1_x.ProjectCoefficient(funcCoef);
4289 nd_x.ProjectCoefficient(funcCoef);
4290 dgv_x.ProjectCoefficient(funcCoef);
4291
4292 Vector f_val(2*dim-3); f_val = 0.0;
4293 Vector h1_gfc_val(2*dim-3); h1_gfc_val = 0.0;
4294 Vector nd_gfc_val(2*dim-3); nd_gfc_val = 0.0;
4295 Vector dgv_gfc_val(2*dim-3); dgv_gfc_val = 0.0;
4296
4297 SECTION("Domain Evaluation 3D")
4298 {
4299 std::cout << "Domain Evaluation 3D" << std::endl;
4300 for (int e = 0; e < mesh.GetNE(); e++)
4301 {
4302 ElementTransformation *T = mesh.GetElementTransformation(e);
4303 const FiniteElement *fe = h1_fespace.GetFE(e);
4304 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4305 2*order + 2);
4306
4307 double h1_err = 0.0;
4308 double nd_err = 0.0;
4309 double dgv_err = 0.0;
4310
4311 for (int j=0; j<ir.GetNPoints(); j++)
4312 {
4313 npts++;
4314 const IntegrationPoint &ip = ir.IntPoint(j);
4315 T->SetIntPoint(&ip);
4316
4317 dFuncCoef.Eval(f_val, *T, ip);
4318 h1_xCoef.Eval(h1_gfc_val, *T, ip);
4319 nd_xCoef.Eval(nd_gfc_val, *T, ip);
4320 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
4321
4322 double h1_dist = Distance(f_val, h1_gfc_val, 2*dim-3);
4323 double nd_dist = Distance(f_val, nd_gfc_val, 2*dim-3);
4324 double dgv_dist = Distance(f_val, dgv_gfc_val, 2*dim-3);
4325
4326 h1_err += h1_dist;
4327 nd_err += nd_dist;
4328 dgv_err += dgv_dist;
4329
4330 if (log > 0 && h1_dist > tol)
4331 {
4332 std::cout << e << ":" << j << " h1 ("
4333 << f_val[0] << "," << f_val[1] << ","
4334 << f_val[2] << ") vs. ("
4335 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
4336 << h1_gfc_val[2] << ") " << h1_dist
4337 << std::endl;
4338 }
4339 if (log > 0 && nd_dist > tol)
4340 {
4341 std::cout << e << ":" << j << " nd ("
4342 << f_val[0] << "," << f_val[1] << ","
4343 << f_val[2] << ") vs. ("
4344 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ","
4345 << nd_gfc_val[2] << ") " << nd_dist
4346 << std::endl;
4347 }
4348 if (log > 0 && dgv_dist > tol)
4349 {
4350 std::cout << e << ":" << j << " dgv ("
4351 << f_val[0] << "," << f_val[1] << ","
4352 << f_val[2] << ") vs. ("
4353 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
4354 << dgv_gfc_val[2] << ") " << dgv_dist
4355 << std::endl;
4356 }
4357 }
4358 h1_err /= ir.GetNPoints();
4359 nd_err /= ir.GetNPoints();
4360 dgv_err /= ir.GetNPoints();
4361
4362 REQUIRE( h1_err == MFEM_Approx(0.0));
4363 REQUIRE( nd_err == MFEM_Approx(0.0));
4364 REQUIRE(dgv_err == MFEM_Approx(0.0));
4365 }
4366 }
4367
4368 SECTION("Boundary Evaluation 3D (H1 Context)")
4369 {
4370 std::cout << "Boundary Evaluation 3D (H1 Context)" << std::endl;
4371 for (int be = 0; be < mesh.GetNBE(); be++)
4372 {
4373 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
4374 const FiniteElement *fe = h1_fespace.GetBE(be);
4375 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4376 2*order + 2);
4377
4378 double h1_err = 0.0;
4379 double nd_err = 0.0;
4380 double dgv_err = 0.0;
4381
4382 for (int j=0; j<ir.GetNPoints(); j++)
4383 {
4384 npts++;
4385 const IntegrationPoint &ip = ir.IntPoint(j);
4386 T->SetIntPoint(&ip);
4387
4388 dFuncCoef.Eval(f_val, *T, ip);
4389 h1_xCoef.Eval(h1_gfc_val, *T, ip);
4390 nd_xCoef.Eval(nd_gfc_val, *T, ip);
4391 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
4392
4393 double h1_dist = Distance(f_val, h1_gfc_val, 2*dim-3);
4394 double nd_dist = Distance(f_val, nd_gfc_val, 2*dim-3);
4395 double dgv_dist = Distance(f_val, dgv_gfc_val, 2*dim-3);
4396
4397 h1_err += h1_dist;
4398 nd_err += nd_dist;
4399 dgv_err += dgv_dist;
4400
4401 if (log > 0 && h1_dist > tol)
4402 {
4403 std::cout << be << ":" << j << " h1 ("
4404 << f_val[0] << "," << f_val[1] << ","
4405 << f_val[2] << ") vs. ("
4406 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
4407 << h1_gfc_val[2] << ") " << h1_dist
4408 << std::endl;
4409 }
4410 if (log > 0 && nd_dist > tol)
4411 {
4412 std::cout << be << ":" << j << " nd ("
4413 << f_val[0] << "," << f_val[1] << ","
4414 << f_val[2] << ") vs. ("
4415 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ","
4416 << nd_gfc_val[2] << ") " << nd_dist
4417 << std::endl;
4418 }
4419 if (log > 0 && dgv_dist > tol)
4420 {
4421 std::cout << be << ":" << j << " dgv ("
4422 << f_val[0] << "," << f_val[1] << ","
4423 << f_val[2] << ") vs. ("
4424 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
4425 << dgv_gfc_val[2] << ") " << dgv_dist
4426 << std::endl;
4427 }
4428 }
4429 h1_err /= ir.GetNPoints();
4430 nd_err /= ir.GetNPoints();
4431 dgv_err /= ir.GetNPoints();
4432
4433 REQUIRE( h1_err == MFEM_Approx(0.0));
4434 REQUIRE( nd_err == MFEM_Approx(0.0));
4435 REQUIRE(dgv_err == MFEM_Approx(0.0));
4436 }
4437 }
4438
4439 SECTION("Boundary Evaluation 3D (DG Context)")
4440 {
4441 std::cout << "Boundary Evaluation 3D (DG Context)" << std::endl;
4442 for (int be = 0; be < mesh.GetNBE(); be++)
4443 {
4444 FaceElementTransformations *T =
4445 mesh.GetBdrFaceTransformations(be);
4446 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
4447 2*order + 2);
4448
4449 double h1_err = 0.0;
4450 double nd_err = 0.0;
4451 double dgv_err = 0.0;
4452
4453 for (int j=0; j<ir.GetNPoints(); j++)
4454 {
4455 npts++;
4456 const IntegrationPoint &ip = ir.IntPoint(j);
4457
4458 T->SetIntPoint(&ip);
4459
4460 dFuncCoef.Eval(f_val, *T, ip);
4461 h1_xCoef.Eval(h1_gfc_val, *T, ip);
4462 nd_xCoef.Eval(nd_gfc_val, *T, ip);
4463 dgv_xCoef.Eval(dgv_gfc_val, *T, ip);
4464
4465 double h1_dist = Distance(f_val, h1_gfc_val, 2*dim-3);
4466 double nd_dist = Distance(f_val, nd_gfc_val, 2*dim-3);
4467 double dgv_dist = Distance(f_val, dgv_gfc_val, 2*dim-3);
4468
4469 h1_err += h1_dist;
4470 nd_err += nd_dist;
4471 dgv_err += dgv_dist;
4472
4473 if (log > 0 && h1_dist > tol)
4474 {
4475 std::cout << be << ":" << j << " h1 ("
4476 << f_val[0] << "," << f_val[1] << ","
4477 << f_val[2] << ") vs. ("
4478 << h1_gfc_val[0] << "," << h1_gfc_val[1] << ","
4479 << h1_gfc_val[2] << ") " << h1_dist
4480 << std::endl;
4481 }
4482 if (log > 0 && nd_dist > tol)
4483 {
4484 std::cout << be << ":" << j << " nd ("
4485 << f_val[0] << "," << f_val[1] << ","
4486 << f_val[2] << ") vs. ("
4487 << nd_gfc_val[0] << "," << nd_gfc_val[1] << ","
4488 << nd_gfc_val[2] << ") " << nd_dist
4489 << std::endl;
4490 }
4491 if (log > 0 && dgv_dist > tol)
4492 {
4493 std::cout << be << ":" << j << " dgv ("
4494 << f_val[0] << "," << f_val[1] << ","
4495 << f_val[2] << ") vs. ("
4496 << dgv_gfc_val[0] << "," << dgv_gfc_val[1] << ","
4497 << dgv_gfc_val[2] << ") " << dgv_dist
4498 << std::endl;
4499 }
4500 }
4501 h1_err /= ir.GetNPoints();
4502 nd_err /= ir.GetNPoints();
4503 dgv_err /= ir.GetNPoints();
4504
4505 REQUIRE( h1_err == MFEM_Approx(0.0));
4506 REQUIRE( nd_err == MFEM_Approx(0.0));
4507 REQUIRE(dgv_err == MFEM_Approx(0.0));
4508 }
4509 }
4510 }
4511 }
4512 std::cout << "Checked GridFunction::GetCurl at "
4513 << npts << " 3D points" << std::endl;
4514 }
4515
4516 TEST_CASE("2D GetDivergence",
4517 "[GridFunction]"
4518 "[DivergenceGridFunctionCoefficient]")
4519 {
4520 int log = 1;
4521 int n = 1;
4522 int dim = 2;
4523 int order = 2;
4524 int npts = 0;
4525
4526 double tol = 1e-6;
4527
4528 for (int type = (int)Element::TRIANGLE;
4529 type <= (int)Element::QUADRILATERAL; type++)
4530 {
4531 Mesh mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
4532
4533 VectorFunctionCoefficient funcCoef(dim, Func_2D_quad);
4534 FunctionCoefficient dFuncCoef(DivFunc_2D_quad);
4535
to_string(type)4536 SECTION("2D GetValue tests for element type " + std::to_string(type))
4537 {
4538 H1_FECollection h1_fec(order, dim);
4539 RT_FECollection rt_fec(order+1, dim);
4540 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
4541 FiniteElement::VALUE);
4542
4543 FiniteElementSpace h1_fespace(&mesh, &h1_fec, dim);
4544 FiniteElementSpace rt_fespace(&mesh, &rt_fec);
4545 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec, dim);
4546
4547 GridFunction h1_x(&h1_fespace);
4548 GridFunction rt_x(&rt_fespace);
4549 GridFunction dgv_x(&dgv_fespace);
4550
4551 DivergenceGridFunctionCoefficient h1_xCoef(&h1_x);
4552 DivergenceGridFunctionCoefficient rt_xCoef(&rt_x);
4553 DivergenceGridFunctionCoefficient dgv_xCoef(&dgv_x);
4554
4555 h1_x.ProjectCoefficient(funcCoef);
4556 rt_x.ProjectCoefficient(funcCoef);
4557 dgv_x.ProjectCoefficient(funcCoef);
4558
4559 SECTION("Domain Evaluation 2D")
4560 {
4561 std::cout << "Domain Evaluation 2D" << std::endl;
4562 for (int e = 0; e < mesh.GetNE(); e++)
4563 {
4564 ElementTransformation *T = mesh.GetElementTransformation(e);
4565 const FiniteElement *fe = h1_fespace.GetFE(e);
4566 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4567 2*order + 2);
4568
4569 double h1_err = 0.0;
4570 double rt_err = 0.0;
4571 double dgv_err = 0.0;
4572
4573 for (int j=0; j<ir.GetNPoints(); j++)
4574 {
4575 npts++;
4576 const IntegrationPoint &ip = ir.IntPoint(j);
4577 T->SetIntPoint(&ip);
4578
4579 double f_val = dFuncCoef.Eval(*T, ip);
4580 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
4581 double rt_gfc_val = rt_xCoef.Eval(*T, ip);
4582 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
4583
4584 h1_err += fabs(f_val - h1_gfc_val);
4585 rt_err += fabs(f_val - rt_gfc_val);
4586 dgv_err += fabs(f_val - dgv_gfc_val);
4587
4588 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
4589 {
4590 std::cout << e << ":" << j << " h1 " << f_val << " "
4591 << h1_gfc_val << " "
4592 << fabs(f_val - h1_gfc_val)
4593 << std::endl;
4594 }
4595 if (log > 0 && fabs(f_val - rt_gfc_val) > tol)
4596 {
4597 std::cout << e << ":" << j << " rt " << f_val << " "
4598 << rt_gfc_val << " "
4599 << fabs(f_val - rt_gfc_val)
4600 << std::endl;
4601 }
4602 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
4603 {
4604 std::cout << e << ":" << j << " dgv " << f_val << " "
4605 << dgv_gfc_val << " "
4606 << fabs(f_val - dgv_gfc_val)
4607 << std::endl;
4608 }
4609 }
4610 h1_err /= ir.GetNPoints();
4611 rt_err /= ir.GetNPoints();
4612 dgv_err /= ir.GetNPoints();
4613
4614 REQUIRE(h1_err == MFEM_Approx(0.0));
4615 REQUIRE(rt_err == MFEM_Approx(0.0));
4616 REQUIRE(dgv_err == MFEM_Approx(0.0));
4617 }
4618 }
4619
4620 SECTION("Boundary Evaluation 2D (H1 Context)")
4621 {
4622 std::cout << "Boundary Evaluation 2D (H1 Context)" << std::endl;
4623 for (int be = 0; be < mesh.GetNBE(); be++)
4624 {
4625 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
4626 const FiniteElement *fe = h1_fespace.GetBE(be);
4627 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4628 2*order + 2);
4629
4630 double h1_err = 0.0;
4631 double rt_err = 0.0;
4632 double dgv_err = 0.0;
4633
4634 for (int j=0; j<ir.GetNPoints(); j++)
4635 {
4636 npts++;
4637 const IntegrationPoint &ip = ir.IntPoint(j);
4638 T->SetIntPoint(&ip);
4639
4640 double f_val = dFuncCoef.Eval(*T, ip);
4641 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
4642 double rt_gfc_val = rt_xCoef.Eval(*T, ip);
4643 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
4644
4645 h1_err += fabs(f_val - h1_gfc_val);
4646 rt_err += fabs(f_val - rt_gfc_val);
4647 dgv_err += fabs(f_val - dgv_gfc_val);
4648
4649 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
4650 {
4651 std::cout << be << ":" << j << " h1 " << f_val << " "
4652 << h1_gfc_val << " "
4653 << fabs(f_val - h1_gfc_val)
4654 << std::endl;
4655 }
4656 if (log > 0 && fabs(f_val - rt_gfc_val) > tol)
4657 {
4658 std::cout << be << ":" << j << " rt " << f_val << " "
4659 << rt_gfc_val << " "
4660 << fabs(f_val - rt_gfc_val)
4661 << std::endl;
4662 }
4663 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
4664 {
4665 std::cout << be << ":" << j << " dgv " << f_val << " "
4666 << dgv_gfc_val << " "
4667 << fabs(f_val - dgv_gfc_val)
4668 << std::endl;
4669 }
4670 }
4671 h1_err /= ir.GetNPoints();
4672 rt_err /= ir.GetNPoints();
4673 dgv_err /= ir.GetNPoints();
4674
4675 REQUIRE(h1_err == MFEM_Approx(0.0));
4676 REQUIRE(rt_err == MFEM_Approx(0.0));
4677 REQUIRE(dgv_err == MFEM_Approx(0.0));
4678 }
4679 }
4680
4681 SECTION("Boundary Evaluation 2D (DG Context)")
4682 {
4683 std::cout << "Boundary Evaluation 2D (DG Context)" << std::endl;
4684 for (int be = 0; be < mesh.GetNBE(); be++)
4685 {
4686 FaceElementTransformations *T =
4687 mesh.GetBdrFaceTransformations(be);
4688 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
4689 2*order + 2);
4690
4691 double h1_err = 0.0;
4692 double rt_err = 0.0;
4693 double dgv_err = 0.0;
4694
4695 for (int j=0; j<ir.GetNPoints(); j++)
4696 {
4697 npts++;
4698 const IntegrationPoint &ip = ir.IntPoint(j);
4699
4700 T->SetIntPoint(&ip);
4701
4702 double f_val = dFuncCoef.Eval(*T, ip);
4703 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
4704 double rt_gfc_val = rt_xCoef.Eval(*T, ip);
4705 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
4706
4707 h1_err += fabs(f_val - h1_gfc_val);
4708 rt_err += fabs(f_val - rt_gfc_val);
4709 dgv_err += fabs(f_val - dgv_gfc_val);
4710
4711 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
4712 {
4713 std::cout << be << ":" << j << " h1 " << f_val << " "
4714 << h1_gfc_val << " "
4715 << fabs(f_val - h1_gfc_val)
4716 << std::endl;
4717 }
4718 if (log > 0 && fabs(f_val - rt_gfc_val) > tol)
4719 {
4720 std::cout << be << ":" << j << " rt " << f_val << " "
4721 << rt_gfc_val << " "
4722 << fabs(f_val - rt_gfc_val)
4723 << std::endl;
4724 }
4725 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
4726 {
4727 std::cout << be << ":" << j << " dgv " << f_val << " "
4728 << dgv_gfc_val << " "
4729 << fabs(f_val - dgv_gfc_val)
4730 << std::endl;
4731 }
4732 }
4733 h1_err /= ir.GetNPoints();
4734 rt_err /= ir.GetNPoints();
4735 dgv_err /= ir.GetNPoints();
4736
4737 REQUIRE(h1_err == MFEM_Approx(0.0));
4738 REQUIRE(rt_err == MFEM_Approx(0.0));
4739 REQUIRE(dgv_err == MFEM_Approx(0.0));
4740 }
4741 }
4742 }
4743 }
4744 std::cout << "Checked GridFunction::GetDivergence at "
4745 << npts << " 2D points" << std::endl;
4746 }
4747
4748 TEST_CASE("3D GetDivergence",
4749 "[GridFunction]"
4750 "[DivergenceGridFunctionCoefficient]")
4751 {
4752 int log = 1;
4753 int n = 1;
4754 int dim = 3;
4755 int order = 2;
4756 int npts = 0;
4757
4758 double tol = 1e-6;
4759
4760 for (int type = (int)Element::TETRAHEDRON;
4761 type <= (int)Element::HEXAHEDRON; type++)
4762 {
4763 Mesh mesh = Mesh::MakeCartesian3D(
4764 n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
4765
4766 VectorFunctionCoefficient funcCoef(dim, Func_3D_quad);
4767 FunctionCoefficient dFuncCoef(DivFunc_3D_quad);
4768
to_string(type)4769 SECTION("3D GetValue tests for element type " + std::to_string(type))
4770 {
4771 H1_FECollection h1_fec(order, dim);
4772 RT_FECollection rt_fec(order+1, dim);
4773 DG_FECollection dgv_fec(order, dim, BasisType::GaussLegendre,
4774 FiniteElement::VALUE);
4775
4776 FiniteElementSpace h1_fespace(&mesh, &h1_fec, dim);
4777 FiniteElementSpace rt_fespace(&mesh, &rt_fec);
4778 FiniteElementSpace dgv_fespace(&mesh, &dgv_fec, dim);
4779
4780 GridFunction h1_x(&h1_fespace);
4781 GridFunction rt_x(&rt_fespace);
4782 GridFunction dgv_x(&dgv_fespace);
4783
4784 DivergenceGridFunctionCoefficient h1_xCoef(&h1_x);
4785 DivergenceGridFunctionCoefficient rt_xCoef(&rt_x);
4786 DivergenceGridFunctionCoefficient dgv_xCoef(&dgv_x);
4787
4788 h1_x.ProjectCoefficient(funcCoef);
4789 rt_x.ProjectCoefficient(funcCoef);
4790 dgv_x.ProjectCoefficient(funcCoef);
4791
4792 SECTION("Domain Evaluation 3D")
4793 {
4794 std::cout << "Domain Evaluation 3D" << std::endl;
4795 for (int e = 0; e < mesh.GetNE(); e++)
4796 {
4797 ElementTransformation *T = mesh.GetElementTransformation(e);
4798 const FiniteElement *fe = h1_fespace.GetFE(e);
4799 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4800 2*order + 2);
4801
4802 double h1_err = 0.0;
4803 double rt_err = 0.0;
4804 double dgv_err = 0.0;
4805
4806 for (int j=0; j<ir.GetNPoints(); j++)
4807 {
4808 npts++;
4809 const IntegrationPoint &ip = ir.IntPoint(j);
4810 T->SetIntPoint(&ip);
4811
4812 double f_val = dFuncCoef.Eval(*T, ip);
4813 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
4814 double rt_gfc_val = rt_xCoef.Eval(*T, ip);
4815 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
4816
4817 h1_err += fabs(f_val - h1_gfc_val);
4818 rt_err += fabs(f_val - rt_gfc_val);
4819 dgv_err += fabs(f_val - dgv_gfc_val);
4820
4821 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
4822 {
4823 std::cout << e << ":" << j << " h1 " << f_val << " "
4824 << h1_gfc_val << " "
4825 << fabs(f_val - h1_gfc_val)
4826 << std::endl;
4827 }
4828 if (log > 0 && fabs(f_val - rt_gfc_val) > tol)
4829 {
4830 std::cout << e << ":" << j << " rt " << f_val << " "
4831 << rt_gfc_val << " "
4832 << fabs(f_val - rt_gfc_val)
4833 << std::endl;
4834 }
4835 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
4836 {
4837 std::cout << e << ":" << j << " dgv " << f_val << " "
4838 << dgv_gfc_val << " "
4839 << fabs(f_val - dgv_gfc_val)
4840 << std::endl;
4841 }
4842 }
4843 h1_err /= ir.GetNPoints();
4844 rt_err /= ir.GetNPoints();
4845 dgv_err /= ir.GetNPoints();
4846
4847 REQUIRE(h1_err == MFEM_Approx(0.0));
4848 REQUIRE(rt_err == MFEM_Approx(0.0));
4849 REQUIRE(dgv_err == MFEM_Approx(0.0));
4850 }
4851 }
4852
4853 SECTION("Boundary Evaluation 3D (H1 Context)")
4854 {
4855 std::cout << "Boundary Evaluation 3D (H1 Context)" << std::endl;
4856 for (int be = 0; be < mesh.GetNBE(); be++)
4857 {
4858 ElementTransformation *T = mesh.GetBdrElementTransformation(be);
4859 const FiniteElement *fe = h1_fespace.GetBE(be);
4860 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
4861 2*order + 2);
4862
4863 double h1_err = 0.0;
4864 double rt_err = 0.0;
4865 double dgv_err = 0.0;
4866
4867 for (int j=0; j<ir.GetNPoints(); j++)
4868 {
4869 npts++;
4870 const IntegrationPoint &ip = ir.IntPoint(j);
4871 T->SetIntPoint(&ip);
4872
4873 double f_val = dFuncCoef.Eval(*T, ip);
4874 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
4875 double rt_gfc_val = rt_xCoef.Eval(*T, ip);
4876 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
4877
4878 h1_err += fabs(f_val - h1_gfc_val);
4879 rt_err += fabs(f_val - rt_gfc_val);
4880 dgv_err += fabs(f_val - dgv_gfc_val);
4881
4882 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
4883 {
4884 std::cout << be << ":" << j << " h1 " << f_val << " "
4885 << h1_gfc_val << " "
4886 << fabs(f_val - h1_gfc_val)
4887 << std::endl;
4888 }
4889 if (log > 0 && fabs(f_val - rt_gfc_val) > tol)
4890 {
4891 std::cout << be << ":" << j << " rt " << f_val << " "
4892 << rt_gfc_val << " "
4893 << fabs(f_val - rt_gfc_val)
4894 << std::endl;
4895 }
4896 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
4897 {
4898 std::cout << be << ":" << j << " dgv " << f_val << " "
4899 << dgv_gfc_val << " "
4900 << fabs(f_val - dgv_gfc_val)
4901 << std::endl;
4902 }
4903 }
4904 h1_err /= ir.GetNPoints();
4905 rt_err /= ir.GetNPoints();
4906 dgv_err /= ir.GetNPoints();
4907
4908 REQUIRE(h1_err == MFEM_Approx(0.0));
4909 REQUIRE(rt_err == MFEM_Approx(0.0));
4910 REQUIRE(dgv_err == MFEM_Approx(0.0));
4911 }
4912 }
4913
4914 SECTION("Boundary Evaluation 3D (DG Context)")
4915 {
4916 std::cout << "Boundary Evaluation 3D (DG Context)" << std::endl;
4917 for (int be = 0; be < mesh.GetNBE(); be++)
4918 {
4919 FaceElementTransformations *T =
4920 mesh.GetBdrFaceTransformations(be);
4921 const IntegrationRule &ir = IntRules.Get(T->GetGeometryType(),
4922 2*order + 2);
4923
4924 double h1_err = 0.0;
4925 double rt_err = 0.0;
4926 double dgv_err = 0.0;
4927
4928 for (int j=0; j<ir.GetNPoints(); j++)
4929 {
4930 npts++;
4931 const IntegrationPoint &ip = ir.IntPoint(j);
4932
4933 T->SetIntPoint(&ip);
4934
4935 double f_val = dFuncCoef.Eval(*T, ip);
4936 double h1_gfc_val = h1_xCoef.Eval(*T, ip);
4937 double rt_gfc_val = rt_xCoef.Eval(*T, ip);
4938 double dgv_gfc_val = dgv_xCoef.Eval(*T, ip);
4939
4940 h1_err += fabs(f_val - h1_gfc_val);
4941 rt_err += fabs(f_val - rt_gfc_val);
4942 dgv_err += fabs(f_val - dgv_gfc_val);
4943
4944 if (log > 0 && fabs(f_val - h1_gfc_val) > tol)
4945 {
4946 std::cout << be << ":" << j << " h1 " << f_val << " "
4947 << h1_gfc_val << " "
4948 << fabs(f_val - h1_gfc_val)
4949 << std::endl;
4950 }
4951 if (log > 0 && fabs(f_val - rt_gfc_val) > tol)
4952 {
4953 std::cout << be << ":" << j << " rt " << f_val << " "
4954 << rt_gfc_val << " "
4955 << fabs(f_val - rt_gfc_val)
4956 << std::endl;
4957 }
4958 if (log > 0 && fabs(f_val - dgv_gfc_val) > tol)
4959 {
4960 std::cout << be << ":" << j << " dgv " << f_val << " "
4961 << dgv_gfc_val << " "
4962 << fabs(f_val - dgv_gfc_val)
4963 << std::endl;
4964 }
4965 }
4966 h1_err /= ir.GetNPoints();
4967 rt_err /= ir.GetNPoints();
4968 dgv_err /= ir.GetNPoints();
4969
4970 REQUIRE(h1_err == MFEM_Approx(0.0));
4971 REQUIRE(rt_err == MFEM_Approx(0.0));
4972 REQUIRE(dgv_err == MFEM_Approx(0.0));
4973 }
4974 }
4975 }
4976 }
4977 std::cout << "Checked GridFunction::GetDivergence at "
4978 << npts << " 3D points" << std::endl;
4979 }
4980
4981 } // namespace get_value
4982