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 lin_interp
18 {
19
f1(const Vector & x)20 double f1(const Vector & x) { return 2.345 * x[0]; }
grad_f1(const Vector & x)21 double grad_f1(const Vector & x) { return 2.345; }
Grad_f1(const Vector & x,Vector & df)22 void Grad_f1(const Vector & x, Vector & df)
23 {
24 df.SetSize(1);
25 df[0] = grad_f1(x);
26 }
27
f2(const Vector & x)28 double f2(const Vector & x) { return 2.345 * x[0] + 3.579 * x[1]; }
F2(const Vector & x,Vector & v)29 void F2(const Vector & x, Vector & v)
30 {
31 v.SetSize(2);
32 v[0] = 1.234 * x[0] - 2.357 * x[1];
33 v[1] = 3.572 * x[0] + 4.321 * x[1];
34 }
35
Grad_f2(const Vector & x,Vector & df)36 void Grad_f2(const Vector & x, Vector & df)
37 {
38 df.SetSize(2);
39 df[0] = 2.345;
40 df[1] = 3.579;
41 }
curlF2(const Vector & x)42 double curlF2(const Vector & x) { return 3.572 + 2.357; }
CurlF2(const Vector & x,Vector & v)43 void CurlF2(const Vector & x, Vector & v)
44 {
45 v.SetSize(1);
46 v[0] = curlF2(x);
47 }
DivF2(const Vector & x)48 double DivF2(const Vector & x) { return 1.234 + 4.321; }
49
f3(const Vector & x)50 double f3(const Vector & x)
51 { return 2.345 * x[0] + 3.579 * x[1] + 4.680 * x[2]; }
F3(const Vector & x,Vector & v)52 void F3(const Vector & x, Vector & v)
53 {
54 v.SetSize(3);
55 v[0] = 1.234 * x[0] - 2.357 * x[1] + 3.572 * x[2];
56 v[1] = 2.537 * x[0] + 4.321 * x[1] - 1.234 * x[2];
57 v[2] = -2.572 * x[0] + 1.321 * x[1] + 3.234 * x[2];
58 }
59
Grad_f3(const Vector & x,Vector & df)60 void Grad_f3(const Vector & x, Vector & df)
61 {
62 df.SetSize(3);
63 df[0] = 2.345;
64 df[1] = 3.579;
65 df[2] = 4.680;
66 }
CurlF3(const Vector & x,Vector & df)67 void CurlF3(const Vector & x, Vector & df)
68 {
69 df.SetSize(3);
70 df[0] = 1.321 + 1.234;
71 df[1] = 3.572 + 2.572;
72 df[2] = 2.537 + 2.357;
73 }
DivF3(const Vector & x)74 double DivF3(const Vector & x)
75 { return 1.234 + 4.321 + 3.234; }
76
g1(const Vector & x)77 double g1(const Vector & x) { return 4.234 * x[0]; }
g2(const Vector & x)78 double g2(const Vector & x) { return 4.234 * x[0] + 3.357 * x[1]; }
g3(const Vector & x)79 double g3(const Vector & x)
80 { return 4.234 * x[0] + 3.357 * x[1] + 1.572 * x[2]; }
81
G2(const Vector & x,Vector & v)82 void G2(const Vector & x, Vector & v)
83 {
84 v.SetSize(2);
85 v[0] = 4.234 * x[0] + 3.357 * x[1];
86 v[1] = 4.537 * x[0] + 1.321 * x[1];
87 }
G3(const Vector & x,Vector & v)88 void G3(const Vector & x, Vector & v)
89 {
90 v.SetSize(3);
91 v[0] = 4.234 * x[0] + 3.357 * x[1] + 1.572 * x[2];
92 v[1] = 4.537 * x[0] + 1.321 * x[1] + 2.234 * x[2];
93 v[2] = 1.572 * x[0] + 2.321 * x[1] + 3.234 * x[2];
94 }
95
fg1(const Vector & x)96 double fg1(const Vector & x) { return f1(x) * g1(x); }
97
fg2(const Vector & x)98 double fg2(const Vector & x) { return f2(x) * g2(x); }
fG2(const Vector & x,Vector & v)99 void fG2(const Vector & x, Vector & v) { G2(x, v); v *= f2(x); }
Fg2(const Vector & x,Vector & v)100 void Fg2(const Vector & x, Vector & v) { F2(x, v); v *= g2(x); }
101
FcrossG2(const Vector & x,Vector & FxG)102 void FcrossG2(const Vector & x, Vector & FxG)
103 {
104 Vector F; F2(x, F);
105 Vector G; G2(x, G);
106 FxG.SetSize(1);
107 FxG(0) = F(0) * G(1) - F(1) * G(0);
108 }
109
FdotG2(const Vector & x)110 double FdotG2(const Vector & x)
111 {
112 Vector F; F2(x, F);
113 Vector G; G2(x, G);
114 return F * G;
115 }
116
fg3(const Vector & x)117 double fg3(const Vector & x) { return f3(x) * g3(x); }
fG3(const Vector & x,Vector & v)118 void fG3(const Vector & x, Vector & v) { G3(x, v); v *= f3(x); }
Fg3(const Vector & x,Vector & v)119 void Fg3(const Vector & x, Vector & v) { F3(x, v); v *= g3(x); }
120
FcrossG3(const Vector & x,Vector & FxG)121 void FcrossG3(const Vector & x, Vector & FxG)
122 {
123 Vector F; F3(x, F);
124 Vector G; G3(x, G);
125 FxG.SetSize(3);
126 FxG(0) = F(1) * G(2) - F(2) * G(1);
127 FxG(1) = F(2) * G(0) - F(0) * G(2);
128 FxG(2) = F(0) * G(1) - F(1) * G(0);
129 }
130
FdotG3(const Vector & x)131 double FdotG3(const Vector & x)
132 {
133 Vector F; F3(x, F);
134 Vector G; G3(x, G);
135 return F * G;
136 }
137
138 TEST_CASE("Identity Linear Interpolators",
139 "[IdentityInterpolator]")
140 {
141 int order_h1 = 1, order_nd = 2, order_rt = 1, order_l2 = 1, n = 3, dim = -1;
142 double tol = 1e-9;
143
144 for (int type = (int)Element::SEGMENT;
145 type <= (int)Element::HEXAHEDRON; type++)
146 {
147 Mesh mesh;
148
149 if (type < (int)Element::TRIANGLE)
150 {
151 dim = 1;
152 mesh = Mesh::MakeCartesian1D(n, 2.0);
153
154 }
155 else if (type < (int)Element::TETRAHEDRON)
156 {
157 dim = 2;
158 mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
159 }
160 else
161 {
162 dim = 3;
163 mesh = Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
164
165 if (type == Element::TETRAHEDRON)
166 {
167 mesh.ReorientTetMesh();
168 }
169 }
170
171 FunctionCoefficient fCoef((dim==1) ? f1 :
172 ((dim==2) ? f2 : f3));
173 VectorFunctionCoefficient dfCoef(dim,
174 (dim==1) ? Grad_f1 :
175 ((dim==2)? Grad_f2 : Grad_f3));
176
to_string(type)177 SECTION("Operators on H1 for element type " + std::to_string(type))
178 {
179 H1_FECollection fec_h1(order_h1, dim);
180 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
181
182 GridFunction f0(&fespace_h1);
183 f0.ProjectCoefficient(fCoef);
184
185 SECTION("Mapping to H1")
186 {
187 H1_FECollection fec_h1p(order_h1+1, dim);
188 FiniteElementSpace fespace_h1p(&mesh, &fec_h1p);
189
190 GridFunction f0p(&fespace_h1p);
191
192 DiscreteLinearOperator Op(&fespace_h1,&fespace_h1p);
193 Op.AddDomainInterpolator(new IdentityInterpolator());
194 Op.Assemble();
195
196 Op.Mult(f0,f0p);
197
198 REQUIRE( f0p.ComputeH1Error(&fCoef, &dfCoef) < tol );
199 }
200 SECTION("Mapping to L2")
201 {
202 L2_FECollection fec_l2(order_l2, dim);
203 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
204
205 GridFunction f1(&fespace_l2);
206
207 DiscreteLinearOperator Op(&fespace_h1,&fespace_l2);
208 Op.AddDomainInterpolator(new IdentityInterpolator());
209 Op.Assemble();
210
211 Op.Mult(f0,f1);
212
213 REQUIRE( f1.ComputeL2Error(fCoef) < tol );
214 }
215 SECTION("Mapping to L2 (INTEGRAL)")
216 {
217 L2_FECollection fec_l2(order_l2, dim,
218 BasisType::GaussLegendre,
219 FiniteElement::INTEGRAL);
220 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
221
222 GridFunction f1(&fespace_l2);
223
224 DiscreteLinearOperator Op(&fespace_h1,&fespace_l2);
225 Op.AddDomainInterpolator(new IdentityInterpolator());
226 Op.Assemble();
227
228 Op.Mult(f0,f1);
229
230 REQUIRE( f1.ComputeL2Error(fCoef) < tol );
231 }
232 }
to_string(type)233 SECTION("Operators on L2 for element type " + std::to_string(type))
234 {
235 L2_FECollection fec_l2(order_l2, dim);
236 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
237
238 GridFunction f0(&fespace_l2);
239 f0.ProjectCoefficient(fCoef);
240
241 SECTION("Mapping to L2")
242 {
243 L2_FECollection fec_l2p(order_l2+1, dim);
244 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
245
246 GridFunction f1(&fespace_l2p);
247
248 DiscreteLinearOperator Op(&fespace_l2,&fespace_l2p);
249 Op.AddDomainInterpolator(new IdentityInterpolator());
250 Op.Assemble();
251
252 Op.Mult(f0,f1);
253
254 REQUIRE( f1.ComputeL2Error(fCoef) < tol );
255 }
256 SECTION("Mapping to L2 (INTEGRAL)")
257 {
258 L2_FECollection fec_l2p(order_l2+1, dim,
259 BasisType::GaussLegendre,
260 FiniteElement::INTEGRAL);
261 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
262
263 GridFunction f1(&fespace_l2p);
264
265 DiscreteLinearOperator Op(&fespace_l2,&fespace_l2p);
266 Op.AddDomainInterpolator(new IdentityInterpolator());
267 Op.Assemble();
268
269 Op.Mult(f0,f1);
270
271 REQUIRE( f1.ComputeL2Error(fCoef) < tol );
272 }
273 }
274 SECTION("Operators on L2 (INTEGRAL) for element type " +
to_string(type)275 std::to_string(type))
276 {
277 L2_FECollection fec_l2(order_l2, dim,
278 BasisType::GaussLegendre,
279 FiniteElement::INTEGRAL);
280 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
281
282 GridFunction f0(&fespace_l2);
283 f0.ProjectCoefficient(fCoef);
284
285 SECTION("Mapping to L2")
286 {
287 L2_FECollection fec_l2p(order_l2+1, dim);
288 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
289
290 GridFunction f1(&fespace_l2p);
291
292 DiscreteLinearOperator Op(&fespace_l2,&fespace_l2p);
293 Op.AddDomainInterpolator(new IdentityInterpolator());
294 Op.Assemble();
295
296 Op.Mult(f0,f1);
297
298 REQUIRE( f1.ComputeL2Error(fCoef) < tol );
299 }
300 SECTION("Mapping to L2 (INTEGRAL)")
301 {
302 L2_FECollection fec_l2p(order_l2+1, dim,
303 BasisType::GaussLegendre,
304 FiniteElement::INTEGRAL);
305 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
306
307 GridFunction f1(&fespace_l2p);
308
309 DiscreteLinearOperator Op(&fespace_l2,&fespace_l2p);
310 Op.AddDomainInterpolator(new IdentityInterpolator());
311 Op.Assemble();
312
313 Op.Mult(f0,f1);
314
315 REQUIRE( f1.ComputeL2Error(fCoef) < tol );
316 }
317 }
318 if (dim > 1)
319 {
320 VectorFunctionCoefficient FCoef(dim,
321 (dim==2) ? F2 : F3);
322 VectorFunctionCoefficient curlFCoef(dim,
323 (dim==2) ? CurlF2 : CurlF3);
324 FunctionCoefficient divFCoef((dim==2) ? DivF2 : DivF3);
325
to_string(type)326 SECTION("Operators on HCurl for element type " + std::to_string(type))
327 {
328 ND_FECollection fec_nd(order_nd, dim);
329 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
330
331 GridFunction f1(&fespace_nd);
332 f1.ProjectCoefficient(FCoef);
333
334 SECTION("Mapping to HCurl")
335 {
336 ND_FECollection fec_ndp(order_nd+1, dim);
337 FiniteElementSpace fespace_ndp(&mesh, &fec_ndp);
338
339 GridFunction f1p(&fespace_ndp);
340
341 DiscreteLinearOperator Op(&fespace_nd,&fespace_ndp);
342 Op.AddDomainInterpolator(new IdentityInterpolator());
343 Op.Assemble();
344
345 Op.Mult(f1,f1p);
346
347 REQUIRE( f1p.ComputeHCurlError(&FCoef, &curlFCoef) < tol );
348 }
349 SECTION("Mapping to L2^d")
350 {
351 L2_FECollection fec_l2(order_l2, dim);
352 FiniteElementSpace fespace_l2(&mesh, &fec_l2, dim);
353
354 GridFunction f2d(&fespace_l2);
355
356 DiscreteLinearOperator Op(&fespace_nd,&fespace_l2);
357 Op.AddDomainInterpolator(new IdentityInterpolator());
358 Op.Assemble();
359
360 Op.Mult(f1,f2d);
361
362 REQUIRE( f2d.ComputeL2Error(FCoef) < tol );
363 }
364 SECTION("Mapping to L2^d (INTEGRAL)")
365 {
366 L2_FECollection fec_l2(order_l2, dim,
367 BasisType::GaussLegendre,
368 FiniteElement::INTEGRAL);
369 FiniteElementSpace fespace_l2(&mesh, &fec_l2, dim);
370
371 GridFunction f2d(&fespace_l2);
372
373 DiscreteLinearOperator Op(&fespace_nd,&fespace_l2);
374 Op.AddDomainInterpolator(new IdentityInterpolator());
375 Op.Assemble();
376
377 Op.Mult(f1,f2d);
378
379 REQUIRE( f2d.ComputeL2Error(FCoef) < tol );
380 }
381 }
to_string(type)382 SECTION("Operators on HDiv for element type " + std::to_string(type))
383 {
384 RT_FECollection fec_rt(order_rt, dim);
385 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
386
387 GridFunction f2(&fespace_rt);
388 f2.ProjectCoefficient(FCoef);
389
390 SECTION("Mapping to HDiv")
391 {
392 RT_FECollection fec_rtp(order_rt+1, dim);
393 FiniteElementSpace fespace_rtp(&mesh, &fec_rtp);
394
395 GridFunction f2p(&fespace_rtp);
396
397 DiscreteLinearOperator Op(&fespace_rt,&fespace_rtp);
398 Op.AddDomainInterpolator(new IdentityInterpolator());
399 Op.Assemble();
400
401 Op.Mult(f2,f2p);
402
403 REQUIRE( f2p.ComputeHDivError(&FCoef, &divFCoef) < tol );
404 }
405 SECTION("Mapping to L2^d")
406 {
407 L2_FECollection fec_l2(order_l2, dim);
408 FiniteElementSpace fespace_l2(&mesh, &fec_l2, dim);
409
410 GridFunction f2d(&fespace_l2);
411
412 DiscreteLinearOperator Op(&fespace_rt,&fespace_l2);
413 Op.AddDomainInterpolator(new IdentityInterpolator());
414 Op.Assemble();
415
416 Op.Mult(f2,f2d);
417
418 REQUIRE( f2d.ComputeL2Error(FCoef) < tol );
419 }
420 SECTION("Mapping to L2^d (INTEGRAL)")
421 {
422 L2_FECollection fec_l2(order_l2, dim,
423 BasisType::GaussLegendre,
424 FiniteElement::INTEGRAL);
425 FiniteElementSpace fespace_l2(&mesh, &fec_l2, dim);
426
427 GridFunction f2d(&fespace_l2);
428
429 DiscreteLinearOperator Op(&fespace_rt,&fespace_l2);
430 Op.AddDomainInterpolator(new IdentityInterpolator());
431 Op.Assemble();
432
433 Op.Mult(f2,f2d);
434
435 REQUIRE( f2d.ComputeL2Error(FCoef) < tol );
436 }
437 }
to_string(type)438 SECTION("Operators on H1^d for element type " + std::to_string(type))
439 {
440 H1_FECollection fec_h1(order_h1, dim);
441 FiniteElementSpace fespace_h1(&mesh, &fec_h1, dim);
442
443 GridFunction f0(&fespace_h1);
444 f0.ProjectCoefficient(FCoef);
445
446 SECTION("Mapping to HCurl")
447 {
448 ND_FECollection fec_ndp(order_nd, dim);
449 FiniteElementSpace fespace_ndp(&mesh, &fec_ndp);
450
451 GridFunction f1(&fespace_ndp);
452
453 DiscreteLinearOperator Op(&fespace_h1,&fespace_ndp);
454 Op.AddDomainInterpolator(new IdentityInterpolator());
455 Op.Assemble();
456
457 Op.Mult(f0,f1);
458
459 REQUIRE( f1.ComputeHCurlError(&FCoef, &curlFCoef) < tol );
460 }
461 SECTION("Mapping to HDiv")
462 {
463 RT_FECollection fec_rtp(order_rt, dim);
464 FiniteElementSpace fespace_rtp(&mesh, &fec_rtp);
465
466 GridFunction f2(&fespace_rtp);
467
468 DiscreteLinearOperator Op(&fespace_h1,&fespace_rtp);
469 Op.AddDomainInterpolator(new IdentityInterpolator());
470 Op.Assemble();
471
472 Op.Mult(f0,f2);
473
474 REQUIRE( f2.ComputeHDivError(&FCoef, &divFCoef) < tol );
475 }
476 }
477 /// The following tests would fail. The reason for the failure would
478 /// not be obvious from the user's point of view. I recommend keeping
479 /// these tests here as a reminder that we should consider supporting
480 /// this, or a very similar, usage.
481 /*
482 SECTION("Mapping to L2^d")
483 {
484 L2_FECollection fec_l2(order_l2, dim);
485 FiniteElementSpace fespace_l2(&mesh, &fec_l2, dim);
486
487 GridFunction f2d(&fespace_l2);
488
489 DiscreteLinearOperator Op(&fespace_h1,&fespace_l2);
490 Op.AddDomainInterpolator(new IdentityInterpolator());
491 Op.Assemble();
492
493 Op.Mult(f0,f2d);
494
495 REQUIRE( f2d.ComputeL2Error(FCoef) < tol );
496 }
497 SECTION("Mapping to L2^d (INTEGRAL)")
498 {
499 L2_FECollection fec_l2(order_l2, dim,
500 BasisType::GaussLegendre,
501 FiniteElement::INTEGRAL);
502 FiniteElementSpace fespace_l2(&mesh, &fec_l2, dim);
503
504 GridFunction f2d(&fespace_l2);
505
506 DiscreteLinearOperator Op(&fespace_h1,&fespace_l2);
507 Op.AddDomainInterpolator(new IdentityInterpolator());
508 Op.Assemble();
509
510 Op.Mult(f0,f2d);
511
512 REQUIRE( f2d.ComputeL2Error(FCoef) < tol );
513 }
514 */
515 }
516 }
517 }
518
519 TEST_CASE("Derivative Linear Interpolators",
520 "[GradientInterpolator]"
521 "[CurlInterpolator]"
522 "[DivergenceInterpolator]")
523 {
524 int order_h1 = 1, order_nd = 1, order_rt = 0, order_l2 = 0, n = 3, dim = -1;
525 double tol = 1e-9;
526
527 for (int type = (int)Element::SEGMENT;
528 type <= (int)Element::HEXAHEDRON; type++)
529 {
530 Mesh mesh;
531
532 if (type < (int)Element::TRIANGLE)
533 {
534 dim = 1;
535 mesh = Mesh::MakeCartesian1D(n, 2.0);
536
537 }
538 else if (type < (int)Element::TETRAHEDRON)
539 {
540 dim = 2;
541 mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
542 }
543 else
544 {
545 dim = 3;
546 mesh = Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
547
548 if (type == Element::TETRAHEDRON)
549 {
550 mesh.ReorientTetMesh();
551 }
552 }
553
554 FunctionCoefficient fCoef((dim==1) ? f1 :
555 ((dim==2) ? f2 : f3));
556 FunctionCoefficient gradfCoef(grad_f1);
557 VectorFunctionCoefficient GradfCoef(dim,
558 (dim==1) ? Grad_f1 :
559 ((dim==2)? Grad_f2 : Grad_f3));
560
to_string(type)561 SECTION("Operators on H1 for element type " + std::to_string(type))
562 {
563 H1_FECollection fec_h1(order_h1, dim);
564 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
565
566 GridFunction f0(&fespace_h1);
567 f0.ProjectCoefficient(fCoef);
568
569 if (dim ==1)
570 {
571 SECTION("Mapping to L2")
572 {
573 L2_FECollection fec_l2(order_l2, dim);
574 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
575
576 GridFunction df0(&fespace_l2);
577
578 DiscreteLinearOperator Op(&fespace_h1,&fespace_l2);
579 Op.AddDomainInterpolator(new GradientInterpolator());
580 Op.Assemble();
581
582 Op.Mult(f0,df0);
583
584 REQUIRE( df0.ComputeL2Error(gradfCoef) < tol );
585 }
586 SECTION("Mapping to L2 (INTEGRAL)")
587 {
588 L2_FECollection fec_l2(order_l2, dim,
589 BasisType::GaussLegendre,
590 FiniteElement::INTEGRAL);
591 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
592
593 GridFunction df0(&fespace_l2);
594
595 DiscreteLinearOperator Op(&fespace_h1,&fespace_l2);
596 Op.AddDomainInterpolator(new GradientInterpolator());
597 Op.Assemble();
598
599 Op.Mult(f0,df0);
600
601 REQUIRE( df0.ComputeL2Error(gradfCoef) < tol );
602 }
603
604 }
605 else
606 {
607 SECTION("Mapping to HCurl")
608 {
609 ND_FECollection fec_nd(order_nd, dim);
610 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
611
612 GridFunction df0(&fespace_nd);
613
614 DiscreteLinearOperator Op(&fespace_h1,&fespace_nd);
615 Op.AddDomainInterpolator(new GradientInterpolator());
616 Op.Assemble();
617
618 Op.Mult(f0,df0);
619
620 REQUIRE( df0.ComputeL2Error(GradfCoef) < tol );
621 }
622 }
623 }
624 if (dim > 1)
625 {
626 VectorFunctionCoefficient FCoef(dim,
627 (dim==2) ? F2 : F3);
628 FunctionCoefficient curlFCoef(curlF2);
629 VectorFunctionCoefficient CurlFCoef(dim,
630 (dim==2) ? CurlF2 : CurlF3);
631 FunctionCoefficient DivFCoef((dim==2) ? DivF2 : DivF3);
632
to_string(type)633 SECTION("Operators on HCurl for element type " + std::to_string(type))
634 {
635 ND_FECollection fec_nd(order_nd, dim);
636 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
637
638 GridFunction F1(&fespace_nd);
639 F1.ProjectCoefficient(FCoef);
640
641 if (dim == 2)
642 {
643 SECTION("Mapping to L2")
644 {
645 L2_FECollection fec_l2(order_l2, dim);
646 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
647
648 GridFunction dF1(&fespace_l2);
649
650 DiscreteLinearOperator Op(&fespace_nd,&fespace_l2);
651 Op.AddDomainInterpolator(new CurlInterpolator());
652 Op.Assemble();
653
654 Op.Mult(F1,dF1);
655
656 REQUIRE( dF1.ComputeL2Error(curlFCoef) < tol );
657 }
658 SECTION("Mapping to L2 (INTEGRAL)")
659 {
660 L2_FECollection fec_l2(order_l2, dim,
661 BasisType::GaussLegendre,
662 FiniteElement::INTEGRAL);
663 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
664
665 GridFunction dF1(&fespace_l2);
666
667 DiscreteLinearOperator Op(&fespace_nd,&fespace_l2);
668 Op.AddDomainInterpolator(new CurlInterpolator());
669 Op.Assemble();
670
671 Op.Mult(F1,dF1);
672
673 REQUIRE( dF1.ComputeL2Error(curlFCoef) < tol );
674 }
675 }
676 else
677 {
678 SECTION("Mapping to HDiv")
679 {
680 RT_FECollection fec_rt(order_rt, dim);
681 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
682
683 GridFunction dF1(&fespace_rt);
684
685 DiscreteLinearOperator Op(&fespace_nd,&fespace_rt);
686 Op.AddDomainInterpolator(new CurlInterpolator());
687 Op.Assemble();
688
689 Op.Mult(F1,dF1);
690
691 REQUIRE( dF1.ComputeL2Error(CurlFCoef) < tol );
692 }
693 }
694 }
to_string(type)695 SECTION("Operators on HDiv for element type " + std::to_string(type))
696 {
697 RT_FECollection fec_rt(order_rt, dim);
698 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
699
700 GridFunction F2(&fespace_rt);
701 F2.ProjectCoefficient(FCoef);
702
703 SECTION("Mapping to L2")
704 {
705 L2_FECollection fec_l2(order_l2, dim);
706 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
707
708 GridFunction dF2(&fespace_l2);
709
710 DiscreteLinearOperator Op(&fespace_rt,&fespace_l2);
711 Op.AddDomainInterpolator(new DivergenceInterpolator());
712 Op.Assemble();
713
714 Op.Mult(F2,dF2);
715
716 REQUIRE( dF2.ComputeL2Error(DivFCoef) < tol );
717 }
718 SECTION("Mapping to L2 (INTEGRAL)")
719 {
720 L2_FECollection fec_l2(order_l2, dim,
721 BasisType::GaussLegendre,
722 FiniteElement::INTEGRAL);
723 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
724
725 GridFunction dF2(&fespace_l2);
726
727 DiscreteLinearOperator Op(&fespace_rt,&fespace_l2);
728 Op.AddDomainInterpolator(new DivergenceInterpolator());
729 Op.Assemble();
730
731 Op.Mult(F2,dF2);
732
733 REQUIRE( dF2.ComputeL2Error(DivFCoef) < tol );
734 }
735 }
736 }
737 }
738 }
739
740 TEST_CASE("Product Linear Interpolators",
741 "[ScalarProductInterpolator]"
742 "[VectorScalarProductInterpolator]"
743 "[ScalarVectorProductInterpolator]"
744 "[ScalarCrossProductInterpolator]"
745 "[VectorCrossProductInterpolator]"
746 "[VectorInnerProductInterpolator]")
747 {
748 int order_h1 = 1, order_nd = 2, order_rt = 2, n = 3, dim = -1;
749 double tol = 1e-9;
750
751 for (int type = (int)Element::SEGMENT;
752 type <= (int)Element::HEXAHEDRON; type++)
753 {
754 Mesh mesh;
755
756 if (type < (int)Element::TRIANGLE)
757 {
758 dim = 1;
759 mesh = Mesh::MakeCartesian1D(n, 2.0);
760
761 }
762 else if (type < (int)Element::TETRAHEDRON)
763 {
764 dim = 2;
765 mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
766 }
767 else
768 {
769 dim = 3;
770 mesh = Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
771
772 if (type == Element::TETRAHEDRON)
773 {
774 mesh.ReorientTetMesh();
775 }
776 }
777
778 FunctionCoefficient fCoef((dim==1) ? f1 :
779 ((dim==2) ? f2 : f3));
780 FunctionCoefficient gCoef((dim==1) ? g1 :
781 ((dim==2) ? g2 : g3));
782 FunctionCoefficient fgCoef((dim==1) ? fg1 :
783 ((dim==2) ? fg2 : fg3));
784
785 VectorFunctionCoefficient FCoef(dim,
786 (dim==2) ? F2 : F3);
787 VectorFunctionCoefficient GCoef(dim,
788 (dim==2) ? G2 : G3);
789
790 FunctionCoefficient FGCoef((dim==2) ? FdotG2 : FdotG3);
791 VectorFunctionCoefficient fGCoef(dim, (dim==2) ? fG2 : fG3);
792 VectorFunctionCoefficient FgCoef(dim, (dim==2) ? Fg2 : Fg3);
793 VectorFunctionCoefficient FxGCoef(dim, (dim==2) ? FcrossG2 : FcrossG3);
794
to_string(type)795 SECTION("Operators on H1 for element type " + std::to_string(type))
796 {
797 H1_FECollection fec_h1(order_h1, dim);
798 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
799
800 GridFunction g0(&fespace_h1);
801 g0.ProjectCoefficient(gCoef);
802
803 SECTION("Mapping H1 to H1")
804 {
805 H1_FECollection fec_h1p(2*order_h1, dim);
806 FiniteElementSpace fespace_h1p(&mesh, &fec_h1p);
807
808 DiscreteLinearOperator Opf0(&fespace_h1,&fespace_h1p);
809 Opf0.AddDomainInterpolator(
810 new ScalarProductInterpolator(fCoef));
811 Opf0.Assemble();
812
813 GridFunction fg0(&fespace_h1p);
814 Opf0.Mult(g0,fg0);
815
816 REQUIRE( fg0.ComputeL2Error(fgCoef) < tol );
817 }
818 if (dim > 1)
819 {
820 SECTION("Mapping to HCurl")
821 {
822 ND_FECollection fec_nd(order_nd, dim);
823 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
824
825 ND_FECollection fec_ndp(order_h1+order_nd, dim);
826 FiniteElementSpace fespace_ndp(&mesh, &fec_ndp);
827
828 DiscreteLinearOperator OpF1(&fespace_h1,&fespace_ndp);
829 OpF1.AddDomainInterpolator(
830 new VectorScalarProductInterpolator(FCoef));
831 OpF1.Assemble();
832
833 GridFunction Fg1(&fespace_ndp);
834 OpF1.Mult(g0,Fg1);
835
836 REQUIRE( Fg1.ComputeL2Error(FgCoef) < tol );
837 }
838 }
839 }
840 if (dim > 1)
841 {
to_string(type)842 SECTION("Operators on HCurl for element type " + std::to_string(type))
843 {
844 ND_FECollection fec_nd(order_nd, dim);
845 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
846
847 GridFunction G1(&fespace_nd);
848 G1.ProjectCoefficient(GCoef);
849
850 SECTION("Mapping HCurl to HCurl")
851 {
852 H1_FECollection fec_h1(order_h1, dim);
853 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
854
855 ND_FECollection fec_ndp(order_nd+order_h1, dim);
856 FiniteElementSpace fespace_ndp(&mesh, &fec_ndp);
857
858 DiscreteLinearOperator Opf0(&fespace_nd,&fespace_ndp);
859 Opf0.AddDomainInterpolator(
860 new ScalarVectorProductInterpolator(fCoef));
861 Opf0.Assemble();
862
863 GridFunction fG1(&fespace_ndp);
864 Opf0.Mult(G1,fG1);
865
866 REQUIRE( fG1.ComputeL2Error(fGCoef) < tol );
867 }
868 if (dim == 2)
869 {
870 SECTION("Mapping to L2")
871 {
872 L2_FECollection fec_l2p(2*order_nd-1, dim);
873 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
874
875 DiscreteLinearOperator OpF1(&fespace_nd,&fespace_l2p);
876 OpF1.AddDomainInterpolator(
877 new ScalarCrossProductInterpolator(FCoef));
878 OpF1.Assemble();
879
880 GridFunction FxG2(&fespace_l2p);
881 OpF1.Mult(G1,FxG2);
882
883 REQUIRE( FxG2.ComputeL2Error(FxGCoef) < tol );
884 }
885 }
886 else
887 {
888 SECTION("Mapping to HDiv")
889 {
890 RT_FECollection fec_rtp(2*order_nd-1, dim);
891 FiniteElementSpace fespace_rtp(&mesh, &fec_rtp);
892
893 DiscreteLinearOperator OpF1(&fespace_nd,&fespace_rtp);
894 OpF1.AddDomainInterpolator(
895 new VectorCrossProductInterpolator(FCoef));
896 OpF1.Assemble();
897
898 GridFunction FxG2(&fespace_rtp);
899 OpF1.Mult(G1,FxG2);
900
901 REQUIRE( FxG2.ComputeL2Error(FxGCoef) < tol );
902 }
903 }
904 SECTION("Mapping to L2")
905 {
906 RT_FECollection fec_rt(order_rt, dim);
907 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
908
909 L2_FECollection fec_l2p(order_nd+order_rt, dim);
910 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
911
912 DiscreteLinearOperator OpF2(&fespace_nd,&fespace_l2p);
913 OpF2.AddDomainInterpolator(
914 new VectorInnerProductInterpolator(FCoef));
915 OpF2.Assemble();
916
917 GridFunction FG3(&fespace_l2p);
918 OpF2.Mult(G1,FG3);
919
920 REQUIRE( FG3.ComputeL2Error(FGCoef) < tol );
921 }
922 }
to_string(type)923 SECTION("Operators on HDiv for element type " + std::to_string(type))
924 {
925 RT_FECollection fec_rt(order_rt, dim);
926 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
927
928 GridFunction G2(&fespace_rt);
929 G2.ProjectCoefficient(GCoef);
930
931 SECTION("Mapping to L2")
932 {
933 ND_FECollection fec_nd(order_nd, dim);
934 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
935
936 L2_FECollection fec_l2p(order_nd+order_rt, dim);
937 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
938
939 DiscreteLinearOperator OpF1(&fespace_rt,&fespace_l2p);
940 OpF1.AddDomainInterpolator(
941 new VectorInnerProductInterpolator(FCoef));
942 OpF1.Assemble();
943
944 GridFunction FG3(&fespace_l2p);
945 OpF1.Mult(G2,FG3);
946
947 REQUIRE( FG3.ComputeL2Error(FGCoef) < tol );
948 }
949 }
950 }
951 }
952 }
953
954 } // namespace lin_interp
955