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 bilininteg_2d
18 {
19
f2(const Vector & x)20 double f2(const Vector & x) { return 2.345 * x[0] + 3.579 * x[1]; }
F2(const Vector & x,Vector & v)21 void F2(const Vector & x, Vector & v)
22 {
23 v.SetSize(2);
24 v[0] = 1.234 * x[0] - 2.357 * x[1];
25 v[1] = 3.572 * x[0] + 4.321 * x[1];
26 }
27
q2(const Vector & x)28 double q2(const Vector & x) { return 4.234 * x[0] + 3.357 * x[1]; }
V2(const Vector & x,Vector & v)29 void V2(const Vector & x, Vector & v)
30 {
31 v.SetSize(2);
32 v[0] = 2.234 * x[0] + 1.357 * x[1];
33 v[1] = 4.572 * x[0] + 3.321 * x[1];
34 }
M2(const Vector & x,DenseMatrix & m)35 void M2(const Vector & x, DenseMatrix & m)
36 {
37 m.SetSize(2);
38
39 m(0,0) = 4.234 * x[0] + 3.357 * x[1];
40 m(0,1) = 0.234 * x[0] + 0.357 * x[1];
41
42 m(1,0) = -0.572 * x[0] - 0.321 * x[1];
43 m(1,1) = 4.537 * x[0] + 1.321 * x[1];
44 }
MT2(const Vector & x,DenseMatrix & m)45 void MT2(const Vector & x, DenseMatrix & m)
46 {
47 M2(x, m); m.Transpose();
48 }
49
qf2(const Vector & x)50 double qf2(const Vector & x) { return q2(x) * f2(x); }
qF2(const Vector & x,Vector & v)51 void qF2(const Vector & x, Vector & v) { F2(x, v); v *= q2(x); }
MF2(const Vector & x,Vector & v)52 void MF2(const Vector & x, Vector & v)
53 {
54 DenseMatrix M(2); M2(x, M);
55 Vector F(2); F2(x, F);
56 v.SetSize(2); M.Mult(F, v);
57 }
DF2(const Vector & x,Vector & v)58 void DF2(const Vector & x, Vector & v)
59 {
60 Vector D(2); V2(x, D);
61 Vector F(2); F2(x, v);
62 v[0] *= D[0]; v[1] *= D[1];
63 }
64
Grad_f2(const Vector & x,Vector & df)65 void Grad_f2(const Vector & x, Vector & df)
66 {
67 df.SetSize(2);
68 df[0] = 2.345;
69 df[1] = 3.579;
70 }
71
Curl_zf2(const Vector & x,Vector & df)72 void Curl_zf2(const Vector & x, Vector & df)
73 {
74 Grad_f2(x, df);
75 double df0 = df[0];
76 df[0] = df[1];
77 df[1] = -df0;
78 }
79
CurlF2(const Vector & x)80 double CurlF2(const Vector & x) { return 3.572 + 2.357; }
DivF2(const Vector & x)81 double DivF2(const Vector & x) { return 1.234 + 4.321; }
82
qGrad_f2(const Vector & x,Vector & df)83 void qGrad_f2(const Vector & x, Vector & df)
84 {
85 Grad_f2(x, df);
86 df *= q2(x);
87 }
DGrad_f2(const Vector & x,Vector & df)88 void DGrad_f2(const Vector & x, Vector & df)
89 {
90 Vector D(2); V2(x, D);
91 Grad_f2(x, df); df[0] *= D[0]; df[1] *= D[1];
92 }
MGrad_f2(const Vector & x,Vector & df)93 void MGrad_f2(const Vector & x, Vector & df)
94 {
95 DenseMatrix M(2); M2(x, M);
96 Vector gradf(2); Grad_f2(x, gradf);
97 M.Mult(gradf, df);
98 }
99
qCurlF2(const Vector & x)100 double qCurlF2(const Vector & x)
101 {
102 return CurlF2(x) * q2(x);
103 }
104
qDivF2(const Vector & x)105 double qDivF2(const Vector & x)
106 {
107 return q2(x) * DivF2(x);
108 }
109
Vf2(const Vector & x,Vector & vf)110 void Vf2(const Vector & x, Vector & vf) { V2(x, vf); vf *= f2(x); }
111
VdotF2(const Vector & x)112 double VdotF2(const Vector & x)
113 {
114 Vector v;
115 Vector f;
116 V2(x, v);
117 F2(x, f);
118 return v * f;
119 }
120
VdotGrad_f2(const Vector & x)121 double VdotGrad_f2(const Vector & x)
122 {
123 Vector v; V2(x, v);
124 Vector gradf; Grad_f2(x, gradf);
125 return v * gradf;
126 }
127
Vcross_f2(const Vector & x,Vector & Vf)128 void Vcross_f2(const Vector & x, Vector & Vf)
129 {
130 Vector v; V2(x, v);
131 Vf.SetSize(2);
132 Vf[0] = v[1]; Vf[1] = -v[0]; Vf *= f2(x);
133 }
134
VcrossF2(const Vector & x)135 double VcrossF2(const Vector & x)
136 {
137 Vector v; V2(x, v);
138 Vector f; F2(x, f);
139 return v(0) * f(1) - v(1) * f(0);
140 }
141
VcrossGrad_f2(const Vector & x)142 double VcrossGrad_f2(const Vector & x)
143 {
144 Vector V; V2(x, V);
145 Vector dF; Grad_f2(x, dF);
146 return V(0) * dF(1) - V(1) * dF(0);
147 }
148
VcrossCurlF2(const Vector & x,Vector & VF)149 void VcrossCurlF2(const Vector & x, Vector & VF)
150 {
151 Vector V; V2(x, V);
152 VF.SetSize(2);
153 VF(0) = V(1); VF(1) = - V(0); VF *= CurlF2(x);
154 }
155
VDivF2(const Vector & x,Vector & VF)156 void VDivF2(const Vector & x, Vector & VF)
157 {
158 V2(x, VF); VF *= DivF2(x);
159 }
160
MapTypeName(FiniteElement::MapType map_type)161 const std::string MapTypeName(FiniteElement::MapType map_type)
162 {
163 switch (map_type)
164 {
165 case FiniteElement::VALUE:
166 return "VALUE";
167 case FiniteElement::INTEGRAL:
168 return "INTEGRAL";
169 case FiniteElement::H_CURL:
170 return "H_CURL";
171 case FiniteElement::H_DIV:
172 return "H_DIV";
173 default:
174 return "UNKNOWN";
175 }
176 }
177
178 TEST_CASE("2D Bilinear Mass Integrators",
179 "[MixedScalarMassIntegrator]"
180 "[MixedScalarIntegrator]"
181 "[BilinearFormIntegrator]"
182 "[NonlinearFormIntegrator]")
183 {
184 int order = 2, n = 1, dim = 2;
185 double cg_rtol = 1e-14;
186 double tol = 1e-9;
187
188 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
189
190 FunctionCoefficient f2_coef(f2);
191 FunctionCoefficient q2_coef(q2);
192 FunctionCoefficient qf2_coef(qf2);
193
194 SECTION("Operators on H1")
195 {
196 H1_FECollection fec_h1(order, dim);
197 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
198
199 GridFunction f_h1(&fespace_h1); f_h1.ProjectCoefficient(f2_coef);
200
201 for (int map_type = (int)FiniteElement::VALUE;
202 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
203 {
204 SECTION("Mapping H1 to L2 (" +
205 MapTypeName((FiniteElement::MapType)map_type) + ")")
206 {
207 L2_FECollection fec_l2(order, dim,
208 BasisType::GaussLegendre,
209 (FiniteElement::MapType)map_type);
210 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
211
212 BilinearForm m_l2(&fespace_l2);
213 m_l2.AddDomainIntegrator(new MassIntegrator());
214 m_l2.Assemble();
215 m_l2.Finalize();
216
217 GridFunction g_l2(&fespace_l2);
218
219 Vector tmp_l2(fespace_l2.GetNDofs());
220
221 SECTION("Without Coefficient")
222 {
223 MixedBilinearForm blf(&fespace_h1, &fespace_l2);
224 blf.AddDomainIntegrator(new MixedScalarMassIntegrator());
225 blf.Assemble();
226 blf.Finalize();
227
228 blf.Mult(f_h1, tmp_l2); g_l2 = 0.0;
229 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
230
231 REQUIRE( g_l2.ComputeL2Error(f2_coef) < tol );
232
233 MixedBilinearForm blfw(&fespace_l2, &fespace_h1);
234 blfw.AddDomainIntegrator(new MixedScalarMassIntegrator());
235 blfw.Assemble();
236 blfw.Finalize();
237
238 SparseMatrix * blfT = Transpose(blfw.SpMat());
239 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
240 REQUIRE( diff->MaxNorm() < tol );
241 }
242 SECTION("With Coefficient")
243 {
244 MixedBilinearForm blf(&fespace_h1, &fespace_l2);
245 blf.AddDomainIntegrator(new MixedScalarMassIntegrator(q2_coef));
246 blf.Assemble();
247 blf.Finalize();
248
249 blf.Mult(f_h1, tmp_l2); g_l2 = 0.0;
250 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
251
252 REQUIRE( g_l2.ComputeL2Error(qf2_coef) < tol );
253
254 MixedBilinearForm blfw(&fespace_l2, &fespace_h1);
255 blfw.AddDomainIntegrator(new MixedScalarMassIntegrator(q2_coef));
256 blfw.Assemble();
257 blfw.Finalize();
258
259 SparseMatrix * blfT = Transpose(blfw.SpMat());
260 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
261 REQUIRE( diff->MaxNorm() < tol );
262 }
263 }
264 }
265 SECTION("Mapping H1 to H1")
266 {
267 BilinearForm m_h1(&fespace_h1);
268 m_h1.AddDomainIntegrator(new MassIntegrator());
269 m_h1.Assemble();
270 m_h1.Finalize();
271
272 GridFunction g_h1(&fespace_h1);
273
274 Vector tmp_h1(fespace_h1.GetNDofs());
275
276 SECTION("Without Coefficient")
277 {
278 MixedBilinearForm blf(&fespace_h1, &fespace_h1);
279 blf.AddDomainIntegrator(new MixedScalarMassIntegrator());
280 blf.Assemble();
281 blf.Finalize();
282
283 blf.Mult(f_h1, tmp_h1); g_h1 = 0.0;
284 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
285
286 REQUIRE( g_h1.ComputeL2Error(f2_coef) < tol );
287 }
288 SECTION("With Coefficient")
289 {
290 MixedBilinearForm blf(&fespace_h1, &fespace_h1);
291 blf.AddDomainIntegrator(new MixedScalarMassIntegrator(q2_coef));
292 blf.Assemble();
293 blf.Finalize();
294
295 blf.Mult(f_h1, tmp_h1); g_h1 = 0.0;
296 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
297
298 REQUIRE( g_h1.ComputeL2Error(qf2_coef) < tol );
299 }
300 }
301 }
302 for (int map_type_d = (int)FiniteElement::VALUE;
303 map_type_d <= (int)FiniteElement::INTEGRAL; map_type_d++)
304 {
305 SECTION("Operators on L2 (" +
306 MapTypeName((FiniteElement::MapType)map_type_d) + ")")
307 {
308 L2_FECollection fec_l2_d(order, dim,
309 BasisType::GaussLegendre,
310 (FiniteElement::MapType)map_type_d);
311 FiniteElementSpace fespace_l2_d(&mesh, &fec_l2_d);
312
313 GridFunction f_l2(&fespace_l2_d); f_l2.ProjectCoefficient(f2_coef);
314
315 for (int map_type_r = (int)FiniteElement::VALUE;
316 map_type_r <= (int)FiniteElement::INTEGRAL; map_type_r++)
317 {
318 SECTION("Mapping L2 (" +
319 MapTypeName((FiniteElement::MapType)map_type_d) +
320 ") to L2 (" +
321 MapTypeName((FiniteElement::MapType)map_type_r) + ")")
322 {
323 L2_FECollection fec_l2_r(order, dim,
324 BasisType::GaussLegendre,
325 (FiniteElement::MapType)map_type_r);
326 FiniteElementSpace fespace_l2_r(&mesh, &fec_l2_r);
327
328 BilinearForm m_l2(&fespace_l2_r);
329 m_l2.AddDomainIntegrator(new MassIntegrator());
330 m_l2.Assemble();
331 m_l2.Finalize();
332
333 GridFunction g_l2(&fespace_l2_r);
334
335 Vector tmp_l2(fespace_l2_r.GetNDofs());
336
337 SECTION("Without Coefficient")
338 {
339 MixedBilinearForm blf(&fespace_l2_d, &fespace_l2_r);
340 blf.AddDomainIntegrator(new MixedScalarMassIntegrator());
341 blf.Assemble();
342 blf.Finalize();
343
344 blf.Mult(f_l2, tmp_l2); g_l2 = 0.0;
345 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
346
347 REQUIRE( g_l2.ComputeL2Error(f2_coef) < tol );
348 }
349 SECTION("With Coefficient")
350 {
351 MixedBilinearForm blf(&fespace_l2_d, &fespace_l2_r);
352 blf.AddDomainIntegrator(new MixedScalarMassIntegrator(q2_coef));
353 blf.Assemble();
354 blf.Finalize();
355
356 blf.Mult(f_l2, tmp_l2); g_l2 = 0.0;
357 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
358
359 REQUIRE( g_l2.ComputeL2Error(qf2_coef) < tol );
360 }
361 }
362 }
363 SECTION("Mapping L2 (" +
364 MapTypeName((FiniteElement::MapType)map_type_d) +
365 ") to H1")
366 {
367 H1_FECollection fec_h1(order, dim);
368 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
369
370 BilinearForm m_h1(&fespace_h1);
371 m_h1.AddDomainIntegrator(new MassIntegrator());
372 m_h1.Assemble();
373 m_h1.Finalize();
374
375 GridFunction g_h1(&fespace_h1);
376
377 Vector tmp_h1(fespace_h1.GetNDofs());
378
379 SECTION("Without Coefficient")
380 {
381 MixedBilinearForm blf(&fespace_l2_d, &fespace_h1);
382 blf.AddDomainIntegrator(new MixedScalarMassIntegrator());
383 blf.Assemble();
384 blf.Finalize();
385
386 blf.Mult(f_l2, tmp_h1); g_h1 = 0.0;
387 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
388
389 REQUIRE( g_h1.ComputeL2Error(f2_coef) < tol );
390
391 MixedBilinearForm blfw(&fespace_h1, &fespace_l2_d);
392 blfw.AddDomainIntegrator(new MixedScalarMassIntegrator());
393 blfw.Assemble();
394 blfw.Finalize();
395
396 SparseMatrix * blfT = Transpose(blfw.SpMat());
397 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
398 REQUIRE( diff->MaxNorm() < tol );
399 }
400 SECTION("With Coefficient")
401 {
402 MixedBilinearForm blf(&fespace_l2_d, &fespace_h1);
403 blf.AddDomainIntegrator(new MixedScalarMassIntegrator(q2_coef));
404 blf.Assemble();
405 blf.Finalize();
406
407 blf.Mult(f_l2, tmp_h1); g_h1 = 0.0;
408 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
409
410 REQUIRE( g_h1.ComputeL2Error(qf2_coef) < tol );
411
412 MixedBilinearForm blfw(&fespace_h1, &fespace_l2_d);
413 blfw.AddDomainIntegrator(new MixedScalarMassIntegrator(q2_coef));
414 blfw.Assemble();
415 blfw.Finalize();
416
417 SparseMatrix * blfT = Transpose(blfw.SpMat());
418 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
419 REQUIRE( diff->MaxNorm() < tol );
420 }
421 }
422 }
423 }
424 }
425
426 TEST_CASE("2D Bilinear Vector Mass Integrators",
427 "[MixedVectorMassIntegrator]"
428 "[MixedVectorIntegrator]"
429 "[BilinearFormIntegrator]"
430 "[NonlinearFormIntegrator]")
431 {
432 int order = 2, n = 1, dim = 2;
433 double cg_rtol = 1e-14;
434 double tol = 1e-9;
435
436 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
437
438 VectorFunctionCoefficient F2_coef(dim, F2);
439 FunctionCoefficient q2_coef(q2);
440 VectorFunctionCoefficient D2_coef(dim, V2);
441 MatrixFunctionCoefficient M2_coef(dim, M2);
442 MatrixFunctionCoefficient MT2_coef(dim, MT2);
443 VectorFunctionCoefficient qF2_coef(dim, qF2);
444 VectorFunctionCoefficient DF2_coef(dim, DF2);
445 VectorFunctionCoefficient MF2_coef(dim, MF2);
446
447 SECTION("Operators on ND")
448 {
449 ND_FECollection fec_nd(order, dim);
450 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
451
452 GridFunction f_nd(&fespace_nd); f_nd.ProjectCoefficient(F2_coef);
453
454 SECTION("Mapping ND to RT")
455 {
456 {
457 // Tests requiring an RT space with same order of
458 // convergence as the ND space
459 RT_FECollection fec_rt(order - 1, dim);
460 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
461
462 BilinearForm m_rt(&fespace_rt);
463 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
464 m_rt.Assemble();
465 m_rt.Finalize();
466
467 GridFunction g_rt(&fespace_rt);
468
469 Vector tmp_rt(fespace_rt.GetNDofs());
470
471 SECTION("Without Coefficient")
472 {
473 MixedBilinearForm blf(&fespace_nd, &fespace_rt);
474 blf.AddDomainIntegrator(new MixedVectorMassIntegrator());
475 blf.Assemble();
476 blf.Finalize();
477
478 blf.Mult(f_nd, tmp_rt); g_rt = 0.0;
479 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
480
481 REQUIRE( g_rt.ComputeL2Error(F2_coef) < tol );
482
483 MixedBilinearForm blfw(&fespace_rt, &fespace_nd);
484 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator());
485 blfw.Assemble();
486 blfw.Finalize();
487
488 SparseMatrix * blfT = Transpose(blfw.SpMat());
489 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
490
491 REQUIRE( diff->MaxNorm() < tol );
492
493 delete blfT;
494 delete diff;
495
496 MixedBilinearForm blfv(&fespace_nd, &fespace_rt);
497 blfv.AddDomainIntegrator(new VectorFEMassIntegrator());
498 blfv.Assemble();
499 blfv.Finalize();
500
501 SparseMatrix * diffv = Add(1.0, blf.SpMat(), -1.0, blfv.SpMat());
502
503 REQUIRE( diffv->MaxNorm() < tol );
504
505 delete diffv;
506 }
507 }
508 {
509 // Tests requiring a higher order RT space
510 RT_FECollection fec_rt(order, dim);
511 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
512
513 BilinearForm m_rt(&fespace_rt);
514 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
515 m_rt.Assemble();
516 m_rt.Finalize();
517
518 GridFunction g_rt(&fespace_rt);
519
520 Vector tmp_rt(fespace_rt.GetNDofs());
521
522 SECTION("With Scalar Coefficient")
523 {
524 MixedBilinearForm blf(&fespace_nd, &fespace_rt);
525 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
526 blf.Assemble();
527 blf.Finalize();
528
529 blf.Mult(f_nd, tmp_rt); g_rt = 0.0;
530 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
531
532 REQUIRE( g_rt.ComputeL2Error(qF2_coef) < tol );
533
534 MixedBilinearForm blfw(&fespace_rt, &fespace_nd);
535 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
536 blfw.Assemble();
537 blfw.Finalize();
538
539 SparseMatrix * blfT = Transpose(blfw.SpMat());
540 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
541
542 REQUIRE( diff->MaxNorm() < tol );
543
544 delete blfT;
545 delete diff;
546 }
547 SECTION("With Diagonal Matrix Coefficient")
548 {
549 MixedBilinearForm blf(&fespace_nd, &fespace_rt);
550 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(D2_coef));
551 blf.Assemble();
552 blf.Finalize();
553
554 blf.Mult(f_nd, tmp_rt); g_rt = 0.0;
555 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
556
557 REQUIRE( g_rt.ComputeL2Error(DF2_coef) < tol );
558
559 MixedBilinearForm blfw(&fespace_rt, &fespace_nd);
560 blfw.AddDomainIntegrator(
561 new MixedVectorMassIntegrator(D2_coef));
562 blfw.Assemble();
563 blfw.Finalize();
564
565 SparseMatrix * blfT = Transpose(blfw.SpMat());
566 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
567
568 REQUIRE( diff->MaxNorm() < tol );
569
570 delete blfT;
571 delete diff;
572 }
573 SECTION("With Matrix Coefficient")
574 {
575 MixedBilinearForm blf(&fespace_nd, &fespace_rt);
576 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(M2_coef));
577 blf.Assemble();
578 blf.Finalize();
579
580 blf.Mult(f_nd, tmp_rt); g_rt = 0.0;
581 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
582
583 REQUIRE( g_rt.ComputeL2Error(MF2_coef) < tol );
584
585 MixedBilinearForm blfw(&fespace_rt, &fespace_nd);
586 blfw.AddDomainIntegrator(
587 new MixedVectorMassIntegrator(MT2_coef));
588 blfw.Assemble();
589 blfw.Finalize();
590
591 SparseMatrix * blfT = Transpose(blfw.SpMat());
592 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
593
594 REQUIRE( diff->MaxNorm() < tol );
595
596 delete blfT;
597 delete diff;
598 }
599 }
600 }
601 SECTION("Mapping ND to ND")
602 {
603 {
604 // Tests requiring an ND test space with same order of
605 // convergence as the ND trial space
606 BilinearForm m_nd(&fespace_nd);
607 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
608 m_nd.Assemble();
609 m_nd.Finalize();
610
611 GridFunction g_nd(&fespace_nd);
612
613 Vector tmp_nd(fespace_nd.GetNDofs());
614
615 SECTION("Without Coefficient")
616 {
617 MixedBilinearForm blf(&fespace_nd, &fespace_nd);
618 blf.AddDomainIntegrator(new MixedVectorMassIntegrator());
619 blf.Assemble();
620 blf.Finalize();
621
622 blf.Mult(f_nd, tmp_nd); g_nd = 0.0;
623 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
624
625 REQUIRE( g_nd.ComputeL2Error(F2_coef) < tol );
626
627 MixedBilinearForm blfw(&fespace_nd, &fespace_nd);
628 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator());
629 blfw.Assemble();
630 blfw.Finalize();
631
632 SparseMatrix * blfT = Transpose(blfw.SpMat());
633 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
634
635 REQUIRE( diff->MaxNorm() < tol );
636
637 delete blfT;
638 delete diff;
639
640 MixedBilinearForm blfv(&fespace_nd, &fespace_nd);
641 blfv.AddDomainIntegrator(new VectorFEMassIntegrator());
642 blfv.Assemble();
643 blfv.Finalize();
644
645 SparseMatrix * diffv = Add(1.0, blf.SpMat(), -1.0, blfv.SpMat());
646
647 REQUIRE( diffv->MaxNorm() < tol );
648
649 delete diffv;
650 }
651 }
652 {
653 // Tests requiring a higher order ND space
654 ND_FECollection fec_ndp(order+1, dim);
655 FiniteElementSpace fespace_ndp(&mesh, &fec_ndp);
656
657 BilinearForm m_ndp(&fespace_ndp);
658 m_ndp.AddDomainIntegrator(new VectorFEMassIntegrator());
659 m_ndp.Assemble();
660 m_ndp.Finalize();
661
662 GridFunction g_ndp(&fespace_ndp);
663
664 Vector tmp_ndp(fespace_ndp.GetNDofs());
665
666 SECTION("With Scalar Coefficient")
667 {
668 MixedBilinearForm blf(&fespace_nd, &fespace_ndp);
669 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
670 blf.Assemble();
671 blf.Finalize();
672
673 blf.Mult(f_nd, tmp_ndp); g_ndp = 0.0;
674 CG(m_ndp, tmp_ndp, g_ndp, 0, 200, cg_rtol * cg_rtol, 0.0);
675
676 REQUIRE( g_ndp.ComputeL2Error(qF2_coef) < tol );
677
678 MixedBilinearForm blfw(&fespace_ndp, &fespace_nd);
679 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
680 blfw.Assemble();
681 blfw.Finalize();
682
683 SparseMatrix * blfT = Transpose(blfw.SpMat());
684 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
685
686 REQUIRE( diff->MaxNorm() < tol );
687
688 delete blfT;
689 delete diff;
690 }
691 SECTION("With Diagonal Matrix Coefficient")
692 {
693 MixedBilinearForm blf(&fespace_nd, &fespace_ndp);
694 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(D2_coef));
695 blf.Assemble();
696 blf.Finalize();
697
698 blf.Mult(f_nd, tmp_ndp); g_ndp = 0.0;
699 CG(m_ndp, tmp_ndp, g_ndp, 0, 200, cg_rtol * cg_rtol, 0.0);
700
701 REQUIRE( g_ndp.ComputeL2Error(DF2_coef) < tol );
702
703 MixedBilinearForm blfw(&fespace_ndp, &fespace_nd);
704 blfw.AddDomainIntegrator(
705 new MixedVectorMassIntegrator(D2_coef));
706 blfw.Assemble();
707 blfw.Finalize();
708
709 SparseMatrix * blfT = Transpose(blfw.SpMat());
710 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
711
712 REQUIRE( diff->MaxNorm() < tol );
713
714 delete blfT;
715 delete diff;
716 }
717 SECTION("With Matrix Coefficient")
718 {
719 MixedBilinearForm blf(&fespace_nd, &fespace_ndp);
720 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(M2_coef));
721 blf.Assemble();
722 blf.Finalize();
723
724 blf.Mult(f_nd, tmp_ndp); g_ndp = 0.0;
725 CG(m_ndp, tmp_ndp, g_ndp, 0, 200, cg_rtol * cg_rtol, 0.0);
726
727 REQUIRE( g_ndp.ComputeL2Error(MF2_coef) < tol );
728
729 MixedBilinearForm blfw(&fespace_ndp, &fespace_nd);
730 blfw.AddDomainIntegrator(
731 new MixedVectorMassIntegrator(MT2_coef));
732 blfw.Assemble();
733 blfw.Finalize();
734
735 SparseMatrix * blfT = Transpose(blfw.SpMat());
736 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
737
738 REQUIRE( diff->MaxNorm() < tol );
739
740 delete blfT;
741 delete diff;
742 }
743 }
744 }
745 }
746 SECTION("Operators on RT")
747 {
748 RT_FECollection fec_rt(order - 1, dim);
749 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
750
751 GridFunction f_rt(&fespace_rt); f_rt.ProjectCoefficient(F2_coef);
752
753 SECTION("Mapping RT to ND")
754 {
755 {
756 // Tests requiring an ND test space with same order of
757 // convergence as the RT trial space
758 ND_FECollection fec_nd(order, dim);
759 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
760
761 BilinearForm m_nd(&fespace_nd);
762 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
763 m_nd.Assemble();
764 m_nd.Finalize();
765
766 GridFunction g_nd(&fespace_nd);
767
768 Vector tmp_nd(fespace_nd.GetNDofs());
769
770 SECTION("Without Coefficient")
771 {
772 MixedBilinearForm blf(&fespace_rt, &fespace_nd);
773 blf.AddDomainIntegrator(new MixedVectorMassIntegrator());
774 blf.Assemble();
775 blf.Finalize();
776
777 blf.Mult(f_rt, tmp_nd); g_nd = 0.0;
778 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
779
780 REQUIRE( g_nd.ComputeL2Error(F2_coef) < tol );
781
782 MixedBilinearForm blfw(&fespace_nd, &fespace_rt);
783 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator());
784 blfw.Assemble();
785 blfw.Finalize();
786
787 SparseMatrix * blfT = Transpose(blfw.SpMat());
788 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
789
790 REQUIRE( diff->MaxNorm() < tol );
791
792 delete blfT;
793 delete diff;
794
795 MixedBilinearForm blfv(&fespace_rt, &fespace_nd);
796 blfv.AddDomainIntegrator(new VectorFEMassIntegrator());
797 blfv.Assemble();
798 blfv.Finalize();
799
800 SparseMatrix * diffv = Add(1.0, blf.SpMat(), -1.0, blfv.SpMat());
801
802 REQUIRE( diffv->MaxNorm() < tol );
803
804 delete diffv;
805 }
806 }
807 {
808 // Tests requiring a higher order ND space
809 ND_FECollection fec_nd(order + 1, dim);
810 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
811
812 BilinearForm m_nd(&fespace_nd);
813 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
814 m_nd.Assemble();
815 m_nd.Finalize();
816
817 GridFunction g_nd(&fespace_nd);
818
819 Vector tmp_nd(fespace_nd.GetNDofs());
820
821 SECTION("With Scalar Coefficient")
822 {
823 MixedBilinearForm blf(&fespace_rt, &fespace_nd);
824 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
825 blf.Assemble();
826 blf.Finalize();
827
828 blf.Mult(f_rt, tmp_nd); g_nd = 0.0;
829 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
830
831 REQUIRE( g_nd.ComputeL2Error(qF2_coef) < tol );
832
833 MixedBilinearForm blfw(&fespace_nd, &fespace_rt);
834 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
835 blfw.Assemble();
836 blfw.Finalize();
837
838 SparseMatrix * blfT = Transpose(blfw.SpMat());
839 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
840
841 REQUIRE( diff->MaxNorm() < tol );
842
843 delete blfT;
844 delete diff;
845 }
846 SECTION("With Diagonal Matrix Coefficient")
847 {
848 MixedBilinearForm blf(&fespace_rt, &fespace_nd);
849 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(D2_coef));
850 blf.Assemble();
851 blf.Finalize();
852
853 blf.Mult(f_rt, tmp_nd); g_nd = 0.0;
854 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
855
856 REQUIRE( g_nd.ComputeL2Error(DF2_coef) < tol );
857
858 MixedBilinearForm blfw(&fespace_nd, &fespace_rt);
859 blfw.AddDomainIntegrator(
860 new MixedVectorMassIntegrator(D2_coef));
861 blfw.Assemble();
862 blfw.Finalize();
863
864 SparseMatrix * blfT = Transpose(blfw.SpMat());
865 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
866
867 REQUIRE( diff->MaxNorm() < tol );
868
869 delete blfT;
870 delete diff;
871 }
872 SECTION("With Matrix Coefficient")
873 {
874 MixedBilinearForm blf(&fespace_rt, &fespace_nd);
875 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(M2_coef));
876 blf.Assemble();
877 blf.Finalize();
878
879 blf.Mult(f_rt, tmp_nd); g_nd = 0.0;
880 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
881
882 REQUIRE( g_nd.ComputeL2Error(MF2_coef) < tol );
883
884 MixedBilinearForm blfw(&fespace_nd, &fespace_rt);
885 blfw.AddDomainIntegrator(
886 new MixedVectorMassIntegrator(MT2_coef));
887 blfw.Assemble();
888 blfw.Finalize();
889
890 SparseMatrix * blfT = Transpose(blfw.SpMat());
891 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
892
893 REQUIRE( diff->MaxNorm() < tol );
894
895 delete blfT;
896 delete diff;
897 }
898 }
899 }
900 SECTION("Mapping RT to RT")
901 {
902 {
903 // Tests requiring an RT test space with same order of
904 // convergence as the RT trial space
905 BilinearForm m_rt(&fespace_rt);
906 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
907 m_rt.Assemble();
908 m_rt.Finalize();
909
910 GridFunction g_rt(&fespace_rt);
911
912 Vector tmp_rt(fespace_rt.GetNDofs());
913
914 SECTION("Without Coefficient")
915 {
916 MixedBilinearForm blf(&fespace_rt, &fespace_rt);
917 blf.AddDomainIntegrator(new MixedVectorMassIntegrator());
918 blf.Assemble();
919 blf.Finalize();
920
921 blf.Mult(f_rt, tmp_rt); g_rt = 0.0;
922 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
923
924 REQUIRE( g_rt.ComputeL2Error(F2_coef) < tol );
925
926 MixedBilinearForm blfw(&fespace_rt, &fespace_rt);
927 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator());
928 blfw.Assemble();
929 blfw.Finalize();
930
931 SparseMatrix * blfT = Transpose(blfw.SpMat());
932 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
933
934 REQUIRE( diff->MaxNorm() < tol );
935
936 delete blfT;
937 delete diff;
938
939 MixedBilinearForm blfv(&fespace_rt, &fespace_rt);
940 blfv.AddDomainIntegrator(new VectorFEMassIntegrator());
941 blfv.Assemble();
942 blfv.Finalize();
943
944 SparseMatrix * diffv = Add(1.0, blf.SpMat(), -1.0, blfv.SpMat());
945
946 REQUIRE( diffv->MaxNorm() < tol );
947
948 delete diffv;
949 }
950 }
951 {
952 // Tests requiring a higher order RT space
953 RT_FECollection fec_rtp(order, dim);
954 FiniteElementSpace fespace_rtp(&mesh, &fec_rtp);
955
956 BilinearForm m_rtp(&fespace_rtp);
957 m_rtp.AddDomainIntegrator(new VectorFEMassIntegrator());
958 m_rtp.Assemble();
959 m_rtp.Finalize();
960
961 GridFunction g_rtp(&fespace_rtp);
962
963 Vector tmp_rtp(fespace_rtp.GetNDofs());
964
965 SECTION("With Scalar Coefficient")
966 {
967 MixedBilinearForm blf(&fespace_rt, &fespace_rtp);
968 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
969 blf.Assemble();
970 blf.Finalize();
971
972 blf.Mult(f_rt, tmp_rtp); g_rtp = 0.0;
973 CG(m_rtp, tmp_rtp, g_rtp, 0, 200, cg_rtol * cg_rtol, 0.0);
974
975 REQUIRE( g_rtp.ComputeL2Error(qF2_coef) < tol );
976
977 MixedBilinearForm blfw(&fespace_rtp, &fespace_rt);
978 blfw.AddDomainIntegrator(new MixedVectorMassIntegrator(q2_coef));
979 blfw.Assemble();
980 blfw.Finalize();
981
982 SparseMatrix * blfT = Transpose(blfw.SpMat());
983 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
984
985 REQUIRE( diff->MaxNorm() < tol );
986
987 delete blfT;
988 delete diff;
989 }
990 SECTION("With Diagonal Matrix Coefficient")
991 {
992 MixedBilinearForm blf(&fespace_rt, &fespace_rtp);
993 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(D2_coef));
994 blf.Assemble();
995 blf.Finalize();
996
997 blf.Mult(f_rt, tmp_rtp); g_rtp = 0.0;
998 CG(m_rtp, tmp_rtp, g_rtp, 0, 200, cg_rtol * cg_rtol, 0.0);
999
1000 REQUIRE( g_rtp.ComputeL2Error(DF2_coef) < tol );
1001
1002 MixedBilinearForm blfw(&fespace_rtp, &fespace_rt);
1003 blfw.AddDomainIntegrator(
1004 new MixedVectorMassIntegrator(D2_coef));
1005 blfw.Assemble();
1006 blfw.Finalize();
1007
1008 SparseMatrix * blfT = Transpose(blfw.SpMat());
1009 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
1010
1011 REQUIRE( diff->MaxNorm() < tol );
1012
1013 delete blfT;
1014 delete diff;
1015 }
1016 SECTION("With Matrix Coefficient")
1017 {
1018 MixedBilinearForm blf(&fespace_rt, &fespace_rtp);
1019 blf.AddDomainIntegrator(new MixedVectorMassIntegrator(M2_coef));
1020 blf.Assemble();
1021 blf.Finalize();
1022
1023 blf.Mult(f_rt, tmp_rtp); g_rtp = 0.0;
1024 CG(m_rtp, tmp_rtp, g_rtp, 0, 200, cg_rtol * cg_rtol, 0.0);
1025
1026 REQUIRE( g_rtp.ComputeL2Error(MF2_coef) < tol );
1027
1028 MixedBilinearForm blfw(&fespace_rtp, &fespace_rt);
1029 blfw.AddDomainIntegrator(
1030 new MixedVectorMassIntegrator(MT2_coef));
1031 blfw.Assemble();
1032 blfw.Finalize();
1033
1034 SparseMatrix * blfT = Transpose(blfw.SpMat());
1035 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
1036
1037 REQUIRE( diff->MaxNorm() < tol );
1038
1039 delete blfT;
1040 delete diff;
1041 }
1042 }
1043 }
1044 }
1045 }
1046
1047 TEST_CASE("2D Bilinear Gradient Integrator",
1048 "[MixedVectorGradientIntegrator]"
1049 "[MixedVectorIntegrator]"
1050 "[BilinearFormIntegrator]"
1051 "[NonlinearFormIntegrator]")
1052 {
1053 int order = 2, n = 1, dim = 2;
1054 double cg_rtol = 1e-14;
1055 double tol = 1e-9;
1056
1057 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1058
1059 FunctionCoefficient f2_coef(f2);
1060 FunctionCoefficient q2_coef(q2);
1061 VectorFunctionCoefficient D2_coef(dim, V2);
1062 MatrixFunctionCoefficient M2_coef(dim, M2);
1063 VectorFunctionCoefficient df2_coef(dim, Grad_f2);
1064 VectorFunctionCoefficient qdf2_coef(dim, qGrad_f2);
1065 VectorFunctionCoefficient Ddf2_coef(dim, DGrad_f2);
1066 VectorFunctionCoefficient Mdf2_coef(dim, MGrad_f2);
1067
1068 SECTION("Operators on H1")
1069 {
1070 H1_FECollection fec_h1(order, dim);
1071 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
1072
1073 GridFunction f_h1(&fespace_h1); f_h1.ProjectCoefficient(f2_coef);
1074
1075 SECTION("Mapping H1 to ND")
1076 {
1077 ND_FECollection fec_nd(order, dim);
1078 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
1079
1080 BilinearForm m_nd(&fespace_nd);
1081 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
1082 m_nd.Assemble();
1083 m_nd.Finalize();
1084
1085 GridFunction g_nd(&fespace_nd);
1086
1087 Vector tmp_nd(fespace_nd.GetNDofs());
1088
1089 SECTION("Without Coefficient")
1090 {
1091 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
1092 blf.AddDomainIntegrator(new MixedVectorGradientIntegrator());
1093 blf.Assemble();
1094 blf.Finalize();
1095 g_nd.ProjectCoefficient(df2_coef);
1096
1097 blf.Mult(f_h1, tmp_nd); g_nd = 0.0;
1098 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
1099
1100 REQUIRE( g_nd.ComputeL2Error(df2_coef) < tol );
1101 }
1102 SECTION("With Scalar Coefficient")
1103 {
1104 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
1105 blf.AddDomainIntegrator(
1106 new MixedVectorGradientIntegrator(q2_coef));
1107 blf.Assemble();
1108 blf.Finalize();
1109
1110 blf.Mult(f_h1, tmp_nd); g_nd = 0.0;
1111 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
1112
1113 REQUIRE( g_nd.ComputeL2Error(qdf2_coef) < tol );
1114 }
1115 SECTION("With Diagonal Matrix Coefficient")
1116 {
1117 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
1118 blf.AddDomainIntegrator(
1119 new MixedVectorGradientIntegrator(D2_coef));
1120 blf.Assemble();
1121 blf.Finalize();
1122
1123 blf.Mult(f_h1, tmp_nd); g_nd = 0.0;
1124 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
1125
1126 REQUIRE( g_nd.ComputeL2Error(Ddf2_coef) < tol );
1127 }
1128 SECTION("With Matrix Coefficient")
1129 {
1130 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
1131 blf.AddDomainIntegrator(
1132 new MixedVectorGradientIntegrator(M2_coef));
1133 blf.Assemble();
1134 blf.Finalize();
1135
1136 blf.Mult(f_h1, tmp_nd); g_nd = 0.0;
1137 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
1138
1139 REQUIRE( g_nd.ComputeL2Error(Mdf2_coef) < tol );
1140 }
1141 }
1142 SECTION("Mapping H1 to RT")
1143 {
1144 RT_FECollection fec_rt(order - 1, dim);
1145 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
1146
1147 BilinearForm m_rt(&fespace_rt);
1148 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
1149 m_rt.Assemble();
1150 m_rt.Finalize();
1151
1152 GridFunction g_rt(&fespace_rt);
1153
1154 Vector tmp_rt(fespace_rt.GetNDofs());
1155
1156 SECTION("Without Coefficient")
1157 {
1158 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
1159 blf.AddDomainIntegrator(new MixedVectorGradientIntegrator());
1160 blf.Assemble();
1161 blf.Finalize();
1162
1163 blf.Mult(f_h1, tmp_rt); g_rt = 0.0;
1164 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
1165
1166 REQUIRE( g_rt.ComputeL2Error(df2_coef) < tol );
1167 }
1168 SECTION("With Scalar Coefficient")
1169 {
1170 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
1171 blf.AddDomainIntegrator(
1172 new MixedVectorGradientIntegrator(q2_coef));
1173 blf.Assemble();
1174 blf.Finalize();
1175
1176 blf.Mult(f_h1, tmp_rt); g_rt = 0.0;
1177 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
1178
1179 REQUIRE( g_rt.ComputeL2Error(qdf2_coef) < tol );
1180 }
1181 SECTION("With Diagonal Matrix Coefficient")
1182 {
1183 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
1184 blf.AddDomainIntegrator(
1185 new MixedVectorGradientIntegrator(D2_coef));
1186 blf.Assemble();
1187 blf.Finalize();
1188
1189 blf.Mult(f_h1, tmp_rt); g_rt = 0.0;
1190 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
1191
1192 REQUIRE( g_rt.ComputeL2Error(Ddf2_coef) < tol );
1193 }
1194 SECTION("With Matrix Coefficient")
1195 {
1196 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
1197 blf.AddDomainIntegrator(
1198 new MixedVectorGradientIntegrator(M2_coef));
1199 blf.Assemble();
1200 blf.Finalize();
1201
1202 blf.Mult(f_h1, tmp_rt); g_rt = 0.0;
1203 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
1204
1205 REQUIRE( g_rt.ComputeL2Error(Mdf2_coef) < tol );
1206 }
1207 }
1208 }
1209 }
1210
1211 TEST_CASE("2D Bilinear Scalar Curl Integrator",
1212 "[MixedScalarCurlIntegrator]"
1213 "[MixedScalarIntegrator]"
1214 "[BilinearFormIntegrator]"
1215 "[NonlinearFormIntegrator]")
1216 {
1217 int order = 2, n = 1, dim = 2;
1218 double cg_rtol = 1e-14;
1219 double tol = 1e-9;
1220
1221 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1222
1223 VectorFunctionCoefficient F2_coef(dim, F2);
1224 FunctionCoefficient q2_coef(q2);
1225 FunctionCoefficient dF2_coef(CurlF2);
1226 FunctionCoefficient qdF2_coef(qCurlF2);
1227
1228 SECTION("Operators on ND")
1229 {
1230 ND_FECollection fec_nd(order, dim);
1231 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
1232
1233 GridFunction f_nd(&fespace_nd); f_nd.ProjectCoefficient(F2_coef);
1234
1235 for (int map_type = (int)FiniteElement::VALUE;
1236 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
1237 {
1238 SECTION("Mapping ND to L2 (" +
1239 MapTypeName((FiniteElement::MapType)map_type) + ")")
1240 {
1241 L2_FECollection fec_l2(order - 1, dim,
1242 BasisType::GaussLegendre,
1243 (FiniteElement::MapType)map_type);
1244 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
1245
1246 BilinearForm m_l2(&fespace_l2);
1247 m_l2.AddDomainIntegrator(new MassIntegrator());
1248 m_l2.Assemble();
1249 m_l2.Finalize();
1250
1251 GridFunction g_l2(&fespace_l2);
1252
1253 Vector tmp_l2(fespace_l2.GetNDofs());
1254
1255 SECTION("Without Coefficient")
1256 {
1257 MixedBilinearForm blf(&fespace_nd, &fespace_l2);
1258 blf.AddDomainIntegrator(new MixedScalarCurlIntegrator());
1259 blf.Assemble();
1260 blf.Finalize();
1261
1262 blf.Mult(f_nd, tmp_l2); g_l2 = 0.0;
1263 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
1264
1265 REQUIRE( g_l2.ComputeL2Error(dF2_coef) < tol );
1266 }
1267 SECTION("With Scalar Coefficient")
1268 {
1269 MixedBilinearForm blf(&fespace_nd, &fespace_l2);
1270 blf.AddDomainIntegrator(
1271 new MixedScalarCurlIntegrator(q2_coef));
1272 blf.Assemble();
1273 blf.Finalize();
1274
1275 blf.Mult(f_nd, tmp_l2); g_l2 = 0.0;
1276 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
1277
1278 REQUIRE( g_l2.ComputeL2Error(qdF2_coef) < tol );
1279 }
1280 }
1281 }
1282 SECTION("Mapping ND to H1")
1283 {
1284 H1_FECollection fec_h1(order, dim);
1285 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
1286
1287 BilinearForm m_h1(&fespace_h1);
1288 m_h1.AddDomainIntegrator(new MassIntegrator());
1289 m_h1.Assemble();
1290 m_h1.Finalize();
1291
1292 GridFunction g_h1(&fespace_h1);
1293
1294 Vector tmp_h1(fespace_h1.GetNDofs());
1295
1296 SECTION("Without Coefficient")
1297 {
1298 MixedBilinearForm blf(&fespace_nd, &fespace_h1);
1299 blf.AddDomainIntegrator(new MixedScalarCurlIntegrator());
1300 blf.Assemble();
1301 blf.Finalize();
1302
1303 blf.Mult(f_nd, tmp_h1); g_h1 = 0.0;
1304 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
1305
1306 REQUIRE( g_h1.ComputeL2Error(dF2_coef) < tol );
1307 }
1308 SECTION("With Scalar Coefficient")
1309 {
1310 MixedBilinearForm blf(&fespace_nd, &fespace_h1);
1311 blf.AddDomainIntegrator(
1312 new MixedScalarCurlIntegrator(q2_coef));
1313 blf.Assemble();
1314 blf.Finalize();
1315
1316 blf.Mult(f_nd, tmp_h1); g_h1 = 0.0;
1317 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
1318
1319 REQUIRE( g_h1.ComputeL2Error(qdF2_coef) < tol );
1320 }
1321 }
1322 }
1323 }
1324
1325 TEST_CASE("2D Bilinear Scalar Cross Product Gradient Integrator",
1326 "[MixedScalarCrossGradIntegrator]"
1327 "[MixedScalarVectorIntegrator]"
1328 "[BilinearFormIntegrator]"
1329 "[NonlinearFormIntegrator]")
1330 {
1331 int order = 2, n = 1, dim = 2;
1332 double cg_rtol = 1e-14;
1333 double tol = 1e-9;
1334
1335 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1336
1337 FunctionCoefficient f2_coef(f2);
1338 VectorFunctionCoefficient V2_coef(dim, V2);
1339 FunctionCoefficient Vxdf2_coef(VcrossGrad_f2);
1340
1341 SECTION("Operators on H1")
1342 {
1343 H1_FECollection fec_h1(order, dim);
1344 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
1345
1346 GridFunction f_h1(&fespace_h1); f_h1.ProjectCoefficient(f2_coef);
1347
1348 for (int map_type = (int)FiniteElement::VALUE;
1349 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
1350 {
1351 SECTION("Mapping H1 to L2 (" +
1352 MapTypeName((FiniteElement::MapType)map_type) + ")")
1353 {
1354 L2_FECollection fec_l2(order, dim,
1355 BasisType::GaussLegendre,
1356 (FiniteElement::MapType)map_type);
1357 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
1358
1359 BilinearForm m_l2(&fespace_l2);
1360 m_l2.AddDomainIntegrator(new MassIntegrator());
1361 m_l2.Assemble();
1362 m_l2.Finalize();
1363
1364 GridFunction g_l2(&fespace_l2);
1365
1366 Vector tmp_l2(fespace_l2.GetNDofs());
1367
1368 SECTION("With Vector Coefficient")
1369 {
1370 MixedBilinearForm blf(&fespace_h1, &fespace_l2);
1371 blf.AddDomainIntegrator(
1372 new MixedScalarCrossGradIntegrator(V2_coef));
1373 blf.Assemble();
1374 blf.Finalize();
1375
1376 blf.Mult(f_h1, tmp_l2); g_l2 = 0.0;
1377 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
1378
1379 REQUIRE( g_l2.ComputeL2Error(Vxdf2_coef) < tol );
1380 }
1381 }
1382 }
1383 SECTION("Mapping H1 to H1")
1384 {
1385 BilinearForm m_h1(&fespace_h1);
1386 m_h1.AddDomainIntegrator(new MassIntegrator());
1387 m_h1.Assemble();
1388 m_h1.Finalize();
1389
1390 GridFunction g_h1(&fespace_h1);
1391
1392 Vector tmp_h1(fespace_h1.GetNDofs());
1393
1394 SECTION("With Vector Coefficient")
1395 {
1396 MixedBilinearForm blf(&fespace_h1, &fespace_h1);
1397 blf.AddDomainIntegrator(
1398 new MixedScalarCrossGradIntegrator(V2_coef));
1399 blf.Assemble();
1400 blf.Finalize();
1401
1402 blf.Mult(f_h1, tmp_h1); g_h1 = 0.0;
1403 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
1404
1405 REQUIRE( g_h1.ComputeL2Error(Vxdf2_coef) < tol );
1406 }
1407 }
1408 }
1409 }
1410
1411 TEST_CASE("2D Bilinear Divergence Integrator",
1412 "[MixedScalarDivergenceIntegrator]"
1413 "[MixedScalarIntegrator]"
1414 "[BilinearFormIntegrator]"
1415 "[NonlinearFormIntegrator]")
1416 {
1417 int order = 2, n = 1, dim = 2;
1418 double cg_rtol = 1e-14;
1419 double tol = 1e-9;
1420
1421 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1422
1423 VectorFunctionCoefficient F2_coef(dim, F2);
1424 FunctionCoefficient q2_coef(q2);
1425 FunctionCoefficient dF2_coef(DivF2);
1426 FunctionCoefficient qdF2_coef(qDivF2);
1427
1428 SECTION("Operators on RT")
1429 {
1430 RT_FECollection fec_rt(order - 1, dim);
1431 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
1432
1433 GridFunction f_rt(&fespace_rt); f_rt.ProjectCoefficient(F2_coef);
1434
1435 for (int map_type = (int)FiniteElement::VALUE;
1436 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
1437 {
1438 SECTION("Mapping RT to L2 (" +
1439 MapTypeName((FiniteElement::MapType)map_type) + ")")
1440 {
1441 L2_FECollection fec_l2(order - 1, dim,
1442 BasisType::GaussLegendre,
1443 (FiniteElement::MapType)map_type);
1444 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
1445
1446 BilinearForm m_l2(&fespace_l2);
1447 m_l2.AddDomainIntegrator(new MassIntegrator());
1448 m_l2.Assemble();
1449 m_l2.Finalize();
1450
1451 GridFunction g_l2(&fespace_l2);
1452
1453 Vector tmp_l2(fespace_l2.GetNDofs());
1454
1455 SECTION("Without Coefficient")
1456 {
1457 MixedBilinearForm blf(&fespace_rt, &fespace_l2);
1458 blf.AddDomainIntegrator(new MixedScalarDivergenceIntegrator());
1459 blf.Assemble();
1460 blf.Finalize();
1461
1462 blf.Mult(f_rt, tmp_l2); g_l2 = 0.0;
1463 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
1464
1465 REQUIRE( g_l2.ComputeL2Error(dF2_coef) < tol );
1466 }
1467 SECTION("With Scalar Coefficient")
1468 {
1469 MixedBilinearForm blf(&fespace_rt, &fespace_l2);
1470 blf.AddDomainIntegrator(
1471 new MixedScalarDivergenceIntegrator(q2_coef));
1472 blf.Assemble();
1473 blf.Finalize();
1474
1475 blf.Mult(f_rt, tmp_l2); g_l2 = 0.0;
1476 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
1477
1478 REQUIRE( g_l2.ComputeL2Error(qdF2_coef) < tol );
1479 }
1480 }
1481 }
1482 SECTION("Mapping RT to H1")
1483 {
1484 H1_FECollection fec_h1(order, dim);
1485 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
1486
1487 BilinearForm m_h1(&fespace_h1);
1488 m_h1.AddDomainIntegrator(new MassIntegrator());
1489 m_h1.Assemble();
1490 m_h1.Finalize();
1491
1492 GridFunction g_h1(&fespace_h1);
1493
1494 Vector tmp_h1(fespace_h1.GetNDofs());
1495
1496 SECTION("Without Coefficient")
1497 {
1498 MixedBilinearForm blf(&fespace_rt, &fespace_h1);
1499 blf.AddDomainIntegrator(new MixedScalarDivergenceIntegrator());
1500 blf.Assemble();
1501 blf.Finalize();
1502
1503 blf.Mult(f_rt, tmp_h1); g_h1 = 0.0;
1504 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
1505
1506 REQUIRE( g_h1.ComputeL2Error(dF2_coef) < tol );
1507 }
1508 SECTION("With Scalar Coefficient")
1509 {
1510 MixedBilinearForm blf(&fespace_rt, &fespace_h1);
1511 blf.AddDomainIntegrator(
1512 new MixedScalarDivergenceIntegrator(q2_coef));
1513 blf.Assemble();
1514 blf.Finalize();
1515
1516 blf.Mult(f_rt, tmp_h1); g_h1 = 0.0;
1517 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
1518
1519 REQUIRE( g_h1.ComputeL2Error(qdF2_coef) < tol );
1520 }
1521 }
1522 }
1523 }
1524
1525 TEST_CASE("2D Bilinear Vector Divergence Integrator",
1526 "[MixedVectorDivergenceIntegrator]"
1527 "[MixedScalarVectorIntegrator]"
1528 "[BilinearFormIntegrator]"
1529 "[NonlinearFormIntegrator]")
1530 {
1531 int order = 2, n = 1, dim = 2;
1532 double cg_rtol = 1e-14;
1533 double tol = 1e-9;
1534
1535 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1536
1537 VectorFunctionCoefficient F2_coef(dim, F2);
1538 VectorFunctionCoefficient V2_coef(dim, V2);
1539 VectorFunctionCoefficient VdF2_coef(dim, VDivF2);
1540
1541 SECTION("Operators on RT")
1542 {
1543 RT_FECollection fec_rt(order - 1, dim);
1544 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
1545
1546 GridFunction f_rt(&fespace_rt); f_rt.ProjectCoefficient(F2_coef);
1547
1548 SECTION("Mapping RT to RT")
1549 {
1550 BilinearForm m_rt(&fespace_rt);
1551 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
1552 m_rt.Assemble();
1553 m_rt.Finalize();
1554
1555 GridFunction g_rt(&fespace_rt);
1556
1557 Vector tmp_rt(fespace_rt.GetNDofs());
1558
1559 SECTION("With Vector Coefficient")
1560 {
1561 MixedBilinearForm blf(&fespace_rt, &fespace_rt);
1562 blf.AddDomainIntegrator(
1563 new MixedVectorDivergenceIntegrator(V2_coef));
1564 blf.Assemble();
1565 blf.Finalize();
1566
1567 blf.Mult(f_rt, tmp_rt); g_rt = 0.0;
1568 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
1569
1570 REQUIRE( g_rt.ComputeL2Error(VdF2_coef) < tol );
1571 }
1572 }
1573 SECTION("Mapping RT to ND")
1574 {
1575 ND_FECollection fec_nd(order, dim);
1576 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
1577
1578 BilinearForm m_nd(&fespace_nd);
1579 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
1580 m_nd.Assemble();
1581 m_nd.Finalize();
1582
1583 GridFunction g_nd(&fespace_nd);
1584
1585 Vector tmp_nd(fespace_nd.GetNDofs());
1586
1587 SECTION("With Vector Coefficient")
1588 {
1589 MixedBilinearForm blf(&fespace_rt, &fespace_nd);
1590 blf.AddDomainIntegrator(
1591 new MixedVectorDivergenceIntegrator(V2_coef));
1592 blf.Assemble();
1593 blf.Finalize();
1594
1595 blf.Mult(f_rt, tmp_nd); g_nd = 0.0;
1596 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
1597
1598 REQUIRE( g_nd.ComputeL2Error(VdF2_coef) < tol );
1599 }
1600 }
1601 }
1602 }
1603
1604 TEST_CASE("2D Bilinear Vector Product Integrators",
1605 "[MixedVectorProductIntegrator]"
1606 "[MixedScalarVectorIntegrator]"
1607 "[BilinearFormIntegrator]"
1608 "[NonlinearFormIntegrator]")
1609 {
1610 int order = 2, n = 1, dim = 2;
1611 double cg_rtol = 1e-14;
1612 double tol = 1e-9;
1613
1614 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1615
1616 FunctionCoefficient f2_coef(f2);
1617 VectorFunctionCoefficient V2_coef(dim, V2);
1618 VectorFunctionCoefficient Vf2_coef(dim, Vf2);
1619
1620 SECTION("Operators on H1")
1621 {
1622 H1_FECollection fec_h1(order, dim);
1623 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
1624
1625 GridFunction f_h1(&fespace_h1); f_h1.ProjectCoefficient(f2_coef);
1626
1627 SECTION("Mapping H1 to ND")
1628 {
1629 ND_FECollection fec_nd(order + 1, dim);
1630 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
1631
1632 BilinearForm m_nd(&fespace_nd);
1633 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
1634 m_nd.Assemble();
1635 m_nd.Finalize();
1636
1637 GridFunction g_nd(&fespace_nd);
1638
1639 Vector tmp_nd(fespace_nd.GetNDofs());
1640
1641 SECTION("With Vector Coefficient")
1642 {
1643 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
1644 blf.AddDomainIntegrator(new MixedVectorProductIntegrator(V2_coef));
1645 blf.Assemble();
1646 blf.Finalize();
1647
1648 blf.Mult(f_h1, tmp_nd); g_nd = 0.0;
1649 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
1650
1651 REQUIRE( g_nd.ComputeL2Error(Vf2_coef) < tol );
1652
1653 MixedBilinearForm blfw(&fespace_nd, &fespace_h1);
1654 blfw.AddDomainIntegrator(
1655 new MixedDotProductIntegrator(V2_coef));
1656 blfw.Assemble();
1657 blfw.Finalize();
1658
1659 SparseMatrix * blfT = Transpose(blfw.SpMat());
1660 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
1661
1662 REQUIRE( diff->MaxNorm() < tol );
1663
1664 delete blfT;
1665 delete diff;
1666 }
1667 }
1668 SECTION("Mapping H1 to RT")
1669 {
1670 RT_FECollection fec_rt(order, dim);
1671 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
1672
1673 BilinearForm m_rt(&fespace_rt);
1674 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
1675 m_rt.Assemble();
1676 m_rt.Finalize();
1677
1678 GridFunction g_rt(&fespace_rt);
1679
1680 Vector tmp_rt(fespace_rt.GetNDofs());
1681
1682 SECTION("With Vector Coefficient")
1683 {
1684 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
1685 blf.AddDomainIntegrator(new MixedVectorProductIntegrator(V2_coef));
1686 blf.Assemble();
1687 blf.Finalize();
1688
1689 blf.Mult(f_h1, tmp_rt); g_rt = 0.0;
1690 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
1691
1692 REQUIRE( g_rt.ComputeL2Error(Vf2_coef) < tol );
1693
1694 MixedBilinearForm blfw(&fespace_rt, &fespace_h1);
1695 blfw.AddDomainIntegrator(
1696 new MixedDotProductIntegrator(V2_coef));
1697 blfw.Assemble();
1698 blfw.Finalize();
1699
1700 SparseMatrix * blfT = Transpose(blfw.SpMat());
1701 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
1702
1703 REQUIRE( diff->MaxNorm() < tol );
1704
1705 delete blfT;
1706 delete diff;
1707 }
1708 }
1709 }
1710 for (int map_type = (int)FiniteElement::VALUE;
1711 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
1712 {
1713 SECTION("Operators on L2 (" +
1714 MapTypeName((FiniteElement::MapType)map_type) + ")")
1715 {
1716 L2_FECollection fec_l2(order, dim,
1717 BasisType::GaussLegendre,
1718 (FiniteElement::MapType)map_type);
1719 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
1720
1721 GridFunction f_l2(&fespace_l2); f_l2.ProjectCoefficient(f2_coef);
1722
1723 SECTION("Mapping L2 (" +
1724 MapTypeName((FiniteElement::MapType)map_type) + ") to ND")
1725 {
1726 ND_FECollection fec_nd(order + 1, dim);
1727 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
1728
1729 BilinearForm m_nd(&fespace_nd);
1730 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
1731 m_nd.Assemble();
1732 m_nd.Finalize();
1733
1734 GridFunction g_nd(&fespace_nd);
1735
1736 Vector tmp_nd(fespace_nd.GetNDofs());
1737
1738 SECTION("With Vector Coefficient")
1739 {
1740 MixedBilinearForm blf(&fespace_l2, &fespace_nd);
1741 blf.AddDomainIntegrator(new MixedVectorProductIntegrator(V2_coef));
1742 blf.Assemble();
1743 blf.Finalize();
1744
1745 blf.Mult(f_l2, tmp_nd); g_nd = 0.0;
1746 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
1747
1748 REQUIRE( g_nd.ComputeL2Error(Vf2_coef) < tol );
1749
1750 MixedBilinearForm blfw(&fespace_nd, &fespace_l2);
1751 blfw.AddDomainIntegrator(
1752 new MixedDotProductIntegrator(V2_coef));
1753 blfw.Assemble();
1754 blfw.Finalize();
1755
1756 SparseMatrix * blfT = Transpose(blfw.SpMat());
1757 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
1758
1759 REQUIRE( diff->MaxNorm() < tol );
1760
1761 delete blfT;
1762 delete diff;
1763 }
1764 }
1765 SECTION("Mapping L2 (" +
1766 MapTypeName((FiniteElement::MapType)map_type) + ") to RT")
1767 {
1768 RT_FECollection fec_rt(order, dim);
1769 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
1770
1771 BilinearForm m_rt(&fespace_rt);
1772 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
1773 m_rt.Assemble();
1774 m_rt.Finalize();
1775
1776 GridFunction g_rt(&fespace_rt);
1777
1778 Vector tmp_rt(fespace_rt.GetNDofs());
1779
1780 SECTION("With Vector Coefficient")
1781 {
1782 MixedBilinearForm blf(&fespace_l2, &fespace_rt);
1783 blf.AddDomainIntegrator(new MixedVectorProductIntegrator(V2_coef));
1784 blf.Assemble();
1785 blf.Finalize();
1786
1787 blf.Mult(f_l2, tmp_rt); g_rt = 0.0;
1788 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
1789
1790 REQUIRE( g_rt.ComputeL2Error(Vf2_coef) < tol );
1791
1792 MixedBilinearForm blfw(&fespace_rt, &fespace_l2);
1793 blfw.AddDomainIntegrator(
1794 new MixedDotProductIntegrator(V2_coef));
1795 blfw.Assemble();
1796 blfw.Finalize();
1797
1798 SparseMatrix * blfT = Transpose(blfw.SpMat());
1799 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
1800
1801 REQUIRE( diff->MaxNorm() < tol );
1802
1803 delete blfT;
1804 delete diff;
1805 }
1806 }
1807 }
1808 }
1809 }
1810
1811 TEST_CASE("2D Bilinear Scalar Cross Product Integrators",
1812 "[MixedScalarCrossProductIntegrator]"
1813 "[MixedScalarVectorIntegrator]"
1814 "[BilinearFormIntegrator]"
1815 "[NonlinearFormIntegrator]")
1816 {
1817 int order = 2, n = 1, dim = 2;
1818 double cg_rtol = 1e-14;
1819 double tol = 1e-9;
1820
1821 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1822
1823 VectorFunctionCoefficient F2_coef(dim, F2);
1824 VectorFunctionCoefficient V2_coef(dim, V2);
1825 FunctionCoefficient VxF2_coef(VcrossF2);
1826
1827 SECTION("Operators on ND")
1828 {
1829 ND_FECollection fec_nd(order, dim);
1830 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
1831
1832 GridFunction f_nd(&fespace_nd); f_nd.ProjectCoefficient(F2_coef);
1833
1834 SECTION("Mapping ND to H1")
1835 {
1836 H1_FECollection fec_h1p(order + 1, dim);
1837 FiniteElementSpace fespace_h1p(&mesh, &fec_h1p);
1838
1839 BilinearForm m_h1p(&fespace_h1p);
1840 m_h1p.AddDomainIntegrator(new MassIntegrator());
1841 m_h1p.Assemble();
1842 m_h1p.Finalize();
1843
1844 GridFunction g_h1p(&fespace_h1p);
1845
1846 Vector tmp_h1p(fespace_h1p.GetNDofs());
1847
1848 SECTION("With Vector Coefficient")
1849 {
1850 MixedBilinearForm blf(&fespace_nd, &fespace_h1p);
1851 blf.AddDomainIntegrator(
1852 new MixedScalarCrossProductIntegrator(V2_coef));
1853 blf.Assemble();
1854 blf.Finalize();
1855
1856 blf.Mult(f_nd, tmp_h1p); g_h1p = 0.0;
1857 CG(m_h1p, tmp_h1p, g_h1p, 0, 200, cg_rtol * cg_rtol, 0.0);
1858
1859 REQUIRE( g_h1p.ComputeL2Error(VxF2_coef) < tol );
1860 }
1861 }
1862 for (int map_type = (int)FiniteElement::VALUE;
1863 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
1864 {
1865 SECTION("Mapping ND to L2 (" +
1866 MapTypeName((FiniteElement::MapType)map_type) + ")")
1867 {
1868 L2_FECollection fec_l2(order, dim,
1869 BasisType::GaussLegendre,
1870 (FiniteElement::MapType)map_type);
1871 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
1872
1873 BilinearForm m_l2(&fespace_l2);
1874 m_l2.AddDomainIntegrator(new MassIntegrator());
1875 m_l2.Assemble();
1876 m_l2.Finalize();
1877
1878 GridFunction g_l2(&fespace_l2);
1879
1880 Vector tmp_l2(fespace_l2.GetNDofs());
1881
1882 SECTION("With Vector Coefficient")
1883 {
1884 MixedBilinearForm blf(&fespace_nd, &fespace_l2);
1885 blf.AddDomainIntegrator(
1886 new MixedScalarCrossProductIntegrator(V2_coef));
1887 blf.Assemble();
1888 blf.Finalize();
1889
1890 blf.Mult(f_nd, tmp_l2); g_l2 = 0.0;
1891 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
1892
1893 REQUIRE( g_l2.ComputeL2Error(VxF2_coef) < tol );
1894 }
1895 }
1896 }
1897 }
1898 SECTION("Operators on RT")
1899 {
1900 RT_FECollection fec_rt(order - 1, dim);
1901 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
1902
1903 GridFunction f_rt(&fespace_rt); f_rt.ProjectCoefficient(F2_coef);
1904
1905 SECTION("Mapping RT to H1")
1906 {
1907 H1_FECollection fec_h1(order + 1, dim);
1908 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
1909
1910 BilinearForm m_h1(&fespace_h1);
1911 m_h1.AddDomainIntegrator(new MassIntegrator());
1912 m_h1.Assemble();
1913 m_h1.Finalize();
1914
1915 GridFunction g_h1(&fespace_h1);
1916
1917 Vector tmp_h1(fespace_h1.GetNDofs());
1918
1919 SECTION("With Vector Coefficient")
1920 {
1921 MixedBilinearForm blf(&fespace_rt, &fespace_h1);
1922 blf.AddDomainIntegrator(
1923 new MixedScalarCrossProductIntegrator(V2_coef));
1924 blf.Assemble();
1925 blf.Finalize();
1926
1927 blf.Mult(f_rt, tmp_h1); g_h1 = 0.0;
1928 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
1929
1930 REQUIRE( g_h1.ComputeL2Error(VxF2_coef) < tol );
1931 }
1932 }
1933 for (int map_type = (int)FiniteElement::VALUE;
1934 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
1935 {
1936 SECTION("Mapping RT to L2 (" +
1937 MapTypeName((FiniteElement::MapType)map_type) + ")")
1938 {
1939 L2_FECollection fec_l2p(order, dim,
1940 BasisType::GaussLegendre,
1941 (FiniteElement::MapType)map_type);
1942 FiniteElementSpace fespace_l2p(&mesh, &fec_l2p);
1943
1944 BilinearForm m_l2p(&fespace_l2p);
1945 m_l2p.AddDomainIntegrator(new MassIntegrator());
1946 m_l2p.Assemble();
1947 m_l2p.Finalize();
1948
1949 GridFunction g_l2p(&fespace_l2p);
1950
1951 Vector tmp_l2p(fespace_l2p.GetNDofs());
1952
1953 SECTION("With Vector Coefficient")
1954 {
1955 MixedBilinearForm blf(&fespace_rt, &fespace_l2p);
1956 blf.AddDomainIntegrator(
1957 new MixedScalarCrossProductIntegrator(V2_coef));
1958 blf.Assemble();
1959 blf.Finalize();
1960
1961 blf.Mult(f_rt, tmp_l2p); g_l2p = 0.0;
1962 CG(m_l2p, tmp_l2p, g_l2p, 0, 200, cg_rtol * cg_rtol, 0.0);
1963
1964 REQUIRE( g_l2p.ComputeL2Error(VxF2_coef) < tol );
1965 }
1966 }
1967 }
1968 }
1969 }
1970
1971 TEST_CASE("2D Bilinear Scalar Weak Cross Product Integrators",
1972 "[MixedScalarWeakCrossProductIntegrator]"
1973 "[MixedScalarVectorIntegrator]"
1974 "[BilinearFormIntegrator]"
1975 "[NonlinearFormIntegrator]")
1976 {
1977 int order = 2, n = 1, dim = 2;
1978 double cg_rtol = 1e-14;
1979 double tol = 1e-9;
1980
1981 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
1982
1983 FunctionCoefficient f2_coef(f2);
1984 VectorFunctionCoefficient V2_coef(dim, V2);
1985 VectorFunctionCoefficient Vxf2_coef(dim, Vcross_f2);
1986
1987 SECTION("Operators on H1")
1988 {
1989 H1_FECollection fec_h1(order, dim);
1990 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
1991
1992 GridFunction f_h1(&fespace_h1); f_h1.ProjectCoefficient(f2_coef);
1993
1994 SECTION("Mapping H1 to ND")
1995 {
1996 ND_FECollection fec_ndp(order + 1, dim);
1997 FiniteElementSpace fespace_ndp(&mesh, &fec_ndp);
1998
1999 BilinearForm m_ndp(&fespace_ndp);
2000 m_ndp.AddDomainIntegrator(new VectorFEMassIntegrator());
2001 m_ndp.Assemble();
2002 m_ndp.Finalize();
2003
2004 GridFunction g_ndp(&fespace_ndp);
2005
2006 Vector tmp_ndp(fespace_ndp.GetNDofs());
2007
2008 SECTION("With Vector Coefficient")
2009 {
2010 MixedBilinearForm blf(&fespace_h1, &fespace_ndp);
2011 blf.AddDomainIntegrator(
2012 new MixedScalarWeakCrossProductIntegrator(V2_coef));
2013 blf.Assemble();
2014 blf.Finalize();
2015
2016 blf.Mult(f_h1, tmp_ndp); g_ndp = 0.0;
2017 CG(m_ndp, tmp_ndp, g_ndp, 0, 200, cg_rtol * cg_rtol, 0.0);
2018
2019 REQUIRE( g_ndp.ComputeL2Error(Vxf2_coef) < tol );
2020
2021 MixedBilinearForm blfw(&fespace_ndp, &fespace_h1);
2022 blfw.AddDomainIntegrator(
2023 new MixedScalarCrossProductIntegrator(V2_coef));
2024 blfw.Assemble();
2025 blfw.Finalize();
2026
2027 SparseMatrix * blfT = Transpose(blfw.SpMat());
2028 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2029
2030 REQUIRE( diff->MaxNorm() < tol );
2031
2032 delete blfT;
2033 delete diff;
2034 }
2035 }
2036 SECTION("Mapping H1 to RT")
2037 {
2038 RT_FECollection fec_rt(order, dim);
2039 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
2040
2041 BilinearForm m_rt(&fespace_rt);
2042 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
2043 m_rt.Assemble();
2044 m_rt.Finalize();
2045
2046 GridFunction g_rt(&fespace_rt);
2047
2048 Vector tmp_rt(fespace_rt.GetNDofs());
2049
2050 SECTION("With Vector Coefficient")
2051 {
2052 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
2053 blf.AddDomainIntegrator(
2054 new MixedScalarWeakCrossProductIntegrator(V2_coef));
2055 blf.Assemble();
2056 blf.Finalize();
2057
2058 blf.Mult(f_h1, tmp_rt); g_rt = 0.0;
2059 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
2060
2061 REQUIRE( g_rt.ComputeL2Error(Vxf2_coef) < tol );
2062
2063 MixedBilinearForm blfw(&fespace_rt, &fespace_h1);
2064 blfw.AddDomainIntegrator(
2065 new MixedScalarCrossProductIntegrator(V2_coef));
2066 blfw.Assemble();
2067 blfw.Finalize();
2068
2069 SparseMatrix * blfT = Transpose(blfw.SpMat());
2070 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2071
2072 REQUIRE( diff->MaxNorm() < tol );
2073
2074 delete blfT;
2075 delete diff;
2076 }
2077 }
2078 }
2079 for (int map_type = (int)FiniteElement::VALUE;
2080 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
2081 {
2082 SECTION("Operators on L2 (" +
2083 MapTypeName((FiniteElement::MapType)map_type) + ")")
2084 {
2085 L2_FECollection fec_l2(order - 1, dim,
2086 BasisType::GaussLegendre,
2087 (FiniteElement::MapType)map_type);
2088 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
2089
2090 GridFunction f_l2(&fespace_l2); f_l2.ProjectCoefficient(f2_coef);
2091
2092 SECTION("Mapping L2 (" +
2093 MapTypeName((FiniteElement::MapType)map_type) + ") to ND")
2094 {
2095 ND_FECollection fec_ndp(order + 1, dim);
2096 FiniteElementSpace fespace_ndp(&mesh, &fec_ndp);
2097
2098 BilinearForm m_ndp(&fespace_ndp);
2099 m_ndp.AddDomainIntegrator(new VectorFEMassIntegrator());
2100 m_ndp.Assemble();
2101 m_ndp.Finalize();
2102
2103 GridFunction g_ndp(&fespace_ndp);
2104
2105 Vector tmp_ndp(fespace_ndp.GetNDofs());
2106
2107 SECTION("With Vector Coefficient")
2108 {
2109 MixedBilinearForm blf(&fespace_l2, &fespace_ndp);
2110 blf.AddDomainIntegrator(
2111 new MixedScalarWeakCrossProductIntegrator(V2_coef));
2112 blf.Assemble();
2113 blf.Finalize();
2114
2115 blf.Mult(f_l2, tmp_ndp); g_ndp = 0.0;
2116 CG(m_ndp, tmp_ndp, g_ndp, 0, 200, cg_rtol * cg_rtol, 0.0);
2117
2118 REQUIRE( g_ndp.ComputeL2Error(Vxf2_coef) < tol );
2119
2120 MixedBilinearForm blfw(&fespace_ndp, &fespace_l2);
2121 blfw.AddDomainIntegrator(
2122 new MixedScalarCrossProductIntegrator(V2_coef));
2123 blfw.Assemble();
2124 blfw.Finalize();
2125
2126 SparseMatrix * blfT = Transpose(blfw.SpMat());
2127 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2128
2129 REQUIRE( diff->MaxNorm() < tol );
2130
2131 delete blfT;
2132 delete diff;
2133 }
2134 }
2135 SECTION("Mapping L2 (" +
2136 MapTypeName((FiniteElement::MapType)map_type) + ") to RT")
2137 {
2138 RT_FECollection fec_rt(order, dim);
2139 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
2140
2141 BilinearForm m_rt(&fespace_rt);
2142 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
2143 m_rt.Assemble();
2144 m_rt.Finalize();
2145
2146 GridFunction g_rt(&fespace_rt);
2147
2148 Vector tmp_rt(fespace_rt.GetNDofs());
2149
2150 SECTION("With Vector Coefficient")
2151 {
2152 MixedBilinearForm blf(&fespace_l2, &fespace_rt);
2153 blf.AddDomainIntegrator(
2154 new MixedScalarWeakCrossProductIntegrator(V2_coef));
2155 blf.Assemble();
2156 blf.Finalize();
2157
2158 blf.Mult(f_l2, tmp_rt); g_rt = 0.0;
2159 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
2160
2161 REQUIRE( g_rt.ComputeL2Error(Vxf2_coef) < tol );
2162
2163 MixedBilinearForm blfw(&fespace_rt, &fespace_l2);
2164 blfw.AddDomainIntegrator(
2165 new MixedScalarCrossProductIntegrator(V2_coef));
2166 blfw.Assemble();
2167 blfw.Finalize();
2168
2169 SparseMatrix * blfT = Transpose(blfw.SpMat());
2170 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2171
2172 REQUIRE( diff->MaxNorm() < tol );
2173
2174 delete blfT;
2175 delete diff;
2176 }
2177 }
2178 }
2179 }
2180 }
2181
2182 TEST_CASE("2D Bilinear Scalar Weak Curl Integrators",
2183 "[MixedScalarWeakCurlIntegrator]"
2184 "[MixedScalarIntegrator]"
2185 "[BilinearFormIntegrator]"
2186 "[NonlinearFormIntegrator]")
2187 {
2188 int order = 2, n = 1, dim = 2;
2189 double tol = 1e-9;
2190
2191 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
2192
2193 FunctionCoefficient q2_coef(q2);
2194
2195 SECTION("Operators on H1")
2196 {
2197 H1_FECollection fec_h1(order, dim);
2198 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2199
2200 SECTION("Mapping H1 to ND")
2201 {
2202 ND_FECollection fec_nd(order, dim);
2203 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
2204
2205 SECTION("Without Coefficient")
2206 {
2207 MixedBilinearForm blf(&fespace_nd, &fespace_h1);
2208 blf.AddDomainIntegrator(
2209 new MixedScalarCurlIntegrator());
2210 blf.Assemble();
2211 blf.Finalize();
2212
2213 MixedBilinearForm blfw(&fespace_h1, &fespace_nd);
2214 blfw.AddDomainIntegrator(
2215 new MixedScalarWeakCurlIntegrator());
2216 blfw.Assemble();
2217 blfw.Finalize();
2218
2219 SparseMatrix * blfT = Transpose(blfw.SpMat());
2220 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
2221
2222 REQUIRE( diff->MaxNorm() < tol );
2223
2224 delete blfT;
2225 delete diff;
2226 }
2227 SECTION("With Scalar Coefficient")
2228 {
2229 MixedBilinearForm blf(&fespace_nd, &fespace_h1);
2230 blf.AddDomainIntegrator(
2231 new MixedScalarCurlIntegrator(q2_coef));
2232 blf.Assemble();
2233 blf.Finalize();
2234
2235 MixedBilinearForm blfw(&fespace_h1, &fespace_nd);
2236 blfw.AddDomainIntegrator(
2237 new MixedScalarWeakCurlIntegrator(q2_coef));
2238 blfw.Assemble();
2239 blfw.Finalize();
2240
2241 SparseMatrix * blfT = Transpose(blfw.SpMat());
2242 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
2243
2244 REQUIRE( diff->MaxNorm() < tol );
2245
2246 delete blfT;
2247 delete diff;
2248 }
2249 }
2250 }
2251 for (int map_type = (int)FiniteElement::VALUE;
2252 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
2253 {
2254 SECTION("Operators on L2 (" +
2255 MapTypeName((FiniteElement::MapType)map_type) + ")")
2256 {
2257 L2_FECollection fec_l2(order - 1, dim,
2258 BasisType::GaussLegendre,
2259 (FiniteElement::MapType)map_type);
2260 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
2261
2262 SECTION("Mapping L2 (" +
2263 MapTypeName((FiniteElement::MapType)map_type) + ") to ND")
2264 {
2265 ND_FECollection fec_nd(order, dim);
2266 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
2267
2268 SECTION("Without Coefficient")
2269 {
2270 MixedBilinearForm blf(&fespace_nd, &fespace_l2);
2271 blf.AddDomainIntegrator(
2272 new MixedScalarCurlIntegrator());
2273 blf.Assemble();
2274 blf.Finalize();
2275
2276 MixedBilinearForm blfw(&fespace_l2, &fespace_nd);
2277 blfw.AddDomainIntegrator(
2278 new MixedScalarWeakCurlIntegrator());
2279 blfw.Assemble();
2280 blfw.Finalize();
2281
2282 SparseMatrix * blfT = Transpose(blfw.SpMat());
2283 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
2284
2285 REQUIRE( diff->MaxNorm() < tol );
2286
2287 delete blfT;
2288 delete diff;
2289 }
2290 SECTION("With Scalar Coefficient")
2291 {
2292 MixedBilinearForm blf(&fespace_nd, &fespace_l2);
2293 blf.AddDomainIntegrator(
2294 new MixedScalarCurlIntegrator(q2_coef));
2295 blf.Assemble();
2296 blf.Finalize();
2297
2298 MixedBilinearForm blfw(&fespace_l2, &fespace_nd);
2299 blfw.AddDomainIntegrator(
2300 new MixedScalarWeakCurlIntegrator(q2_coef));
2301 blfw.Assemble();
2302 blfw.Finalize();
2303
2304 SparseMatrix * blfT = Transpose(blfw.SpMat());
2305 SparseMatrix * diff = Add(1.0, blf.SpMat(), -1.0, *blfT);
2306
2307 REQUIRE( diff->MaxNorm() < tol );
2308
2309 delete blfT;
2310 delete diff;
2311 }
2312 }
2313 }
2314 }
2315 }
2316
2317 TEST_CASE("2D Bilinear Scalar Weak Gradient Integrators",
2318 "[MixedScalarWeakGradientIntegrator]"
2319 "[MixedScalarIntegrator]"
2320 "[BilinearFormIntegrator]"
2321 "[NonlinearFormIntegrator]")
2322 {
2323 int order = 2, n = 1, dim = 2;
2324 double tol = 1e-9;
2325
2326 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
2327
2328 FunctionCoefficient q2_coef(q2);
2329
2330 SECTION("Operators on H1")
2331 {
2332 H1_FECollection fec_h1(order, dim);
2333 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2334
2335 SECTION("Mapping H1 to RT")
2336 {
2337 RT_FECollection fec_rt(order - 1, dim);
2338 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
2339
2340 SECTION("Without Coefficient")
2341 {
2342 MixedBilinearForm blf(&fespace_rt, &fespace_h1);
2343 blf.AddDomainIntegrator(
2344 new MixedScalarDivergenceIntegrator());
2345 blf.Assemble();
2346 blf.Finalize();
2347
2348 MixedBilinearForm blfw(&fespace_h1, &fespace_rt);
2349 blfw.AddDomainIntegrator(
2350 new MixedScalarWeakGradientIntegrator());
2351 blfw.Assemble();
2352 blfw.Finalize();
2353
2354 SparseMatrix * blfT = Transpose(blfw.SpMat());
2355 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2356
2357 REQUIRE( diff->MaxNorm() < tol );
2358
2359 delete blfT;
2360 delete diff;
2361 }
2362 SECTION("With Scalar Coefficient")
2363 {
2364 MixedBilinearForm blf(&fespace_rt, &fespace_h1);
2365 blf.AddDomainIntegrator(
2366 new MixedScalarDivergenceIntegrator(q2_coef));
2367 blf.Assemble();
2368 blf.Finalize();
2369
2370 MixedBilinearForm blfw(&fespace_h1, &fespace_rt);
2371 blfw.AddDomainIntegrator(
2372 new MixedScalarWeakGradientIntegrator(q2_coef));
2373 blfw.Assemble();
2374 blfw.Finalize();
2375
2376 SparseMatrix * blfT = Transpose(blfw.SpMat());
2377 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2378
2379 REQUIRE( diff->MaxNorm() < tol );
2380
2381 delete blfT;
2382 delete diff;
2383 }
2384 }
2385 }
2386 for (int map_type = (int)FiniteElement::VALUE;
2387 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
2388 {
2389 SECTION("Operators on L2 (" +
2390 MapTypeName((FiniteElement::MapType)map_type) + ")")
2391 {
2392 L2_FECollection fec_l2(order - 1, dim,
2393 BasisType::GaussLegendre,
2394 (FiniteElement::MapType)map_type);
2395 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
2396
2397 SECTION("Mapping L2 (" +
2398 MapTypeName((FiniteElement::MapType)map_type) + ") to RT")
2399 {
2400 RT_FECollection fec_rt(order - 1, dim);
2401 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
2402
2403 SECTION("Without Coefficient")
2404 {
2405 MixedBilinearForm blf(&fespace_rt, &fespace_l2);
2406 blf.AddDomainIntegrator(
2407 new MixedScalarDivergenceIntegrator());
2408 blf.Assemble();
2409 blf.Finalize();
2410
2411 MixedBilinearForm blfw(&fespace_l2, &fespace_rt);
2412 blfw.AddDomainIntegrator(
2413 new MixedScalarWeakGradientIntegrator());
2414 blfw.Assemble();
2415 blfw.Finalize();
2416
2417 SparseMatrix * blfT = Transpose(blfw.SpMat());
2418 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2419
2420 REQUIRE( diff->MaxNorm() < tol );
2421
2422 delete blfT;
2423 delete diff;
2424 }
2425 SECTION("With Scalar Coefficient")
2426 {
2427 MixedBilinearForm blf(&fespace_rt, &fespace_l2);
2428 blf.AddDomainIntegrator(
2429 new MixedScalarDivergenceIntegrator(q2_coef));
2430 blf.Assemble();
2431 blf.Finalize();
2432
2433 MixedBilinearForm blfw(&fespace_l2, &fespace_rt);
2434 blfw.AddDomainIntegrator(
2435 new MixedScalarWeakGradientIntegrator(q2_coef));
2436 blfw.Assemble();
2437 blfw.Finalize();
2438
2439 SparseMatrix * blfT = Transpose(blfw.SpMat());
2440 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2441
2442 REQUIRE( diff->MaxNorm() < tol );
2443
2444 delete blfT;
2445 delete diff;
2446 }
2447 }
2448 }
2449 }
2450 }
2451
2452 TEST_CASE("2D Bilinear Directional Derivative Integrators",
2453 "[MixedDirectionalDerivativeIntegrator]"
2454 "[MixedScalarVectorIntegrator]"
2455 "[BilinearFormIntegrator]"
2456 "[NonlinearFormIntegrator]")
2457 {
2458 int order = 2, n = 1, dim = 2;
2459 double cg_rtol = 1e-14;
2460 double tol = 1e-9;
2461
2462 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
2463
2464 FunctionCoefficient f2_coef(f2);
2465 VectorFunctionCoefficient V2_coef(dim, V2);
2466 FunctionCoefficient Vdf2_coef(VdotGrad_f2);
2467
2468 SECTION("Operators on H1")
2469 {
2470 H1_FECollection fec_h1(order, dim);
2471 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2472
2473 GridFunction f_h1(&fespace_h1); f_h1.ProjectCoefficient(f2_coef);
2474
2475 SECTION("Mapping H1 to ND")
2476 {
2477 BilinearForm m_h1(&fespace_h1);
2478 m_h1.AddDomainIntegrator(new MassIntegrator());
2479 m_h1.Assemble();
2480 m_h1.Finalize();
2481
2482 GridFunction g_h1(&fespace_h1);
2483
2484 Vector tmp_h1(fespace_h1.GetNDofs());
2485
2486 SECTION("With Vector Coefficient")
2487 {
2488 MixedBilinearForm blf(&fespace_h1, &fespace_h1);
2489 blf.AddDomainIntegrator(
2490 new MixedDirectionalDerivativeIntegrator(V2_coef));
2491 blf.Assemble();
2492 blf.Finalize();
2493
2494 blf.Mult(f_h1, tmp_h1); g_h1 = 0.0;
2495 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
2496
2497 REQUIRE( g_h1.ComputeL2Error(Vdf2_coef) < tol );
2498 }
2499 }
2500 for (int map_type = (int)FiniteElement::VALUE;
2501 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
2502 {
2503 SECTION("Mapping H1 to L2 (" +
2504 MapTypeName((FiniteElement::MapType)map_type) + ")")
2505 {
2506 L2_FECollection fec_l2(order - 1, dim,
2507 BasisType::GaussLegendre,
2508 (FiniteElement::MapType)map_type);
2509 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
2510
2511 BilinearForm m_l2(&fespace_l2);
2512 m_l2.AddDomainIntegrator(new MassIntegrator());
2513 m_l2.Assemble();
2514 m_l2.Finalize();
2515
2516 GridFunction g_l2(&fespace_l2);
2517
2518 Vector tmp_l2(fespace_l2.GetNDofs());
2519
2520 SECTION("With Vector Coefficient")
2521 {
2522 MixedBilinearForm blf(&fespace_h1, &fespace_l2);
2523 blf.AddDomainIntegrator(
2524 new MixedDirectionalDerivativeIntegrator(V2_coef));
2525 blf.Assemble();
2526 blf.Finalize();
2527
2528 blf.Mult(f_h1, tmp_l2); g_l2 = 0.0;
2529 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
2530
2531 REQUIRE( g_l2.ComputeL2Error(Vdf2_coef) < tol );
2532 }
2533 }
2534 }
2535 }
2536 }
2537
2538
2539 TEST_CASE("2D Bilinear Scalar Weak Divergence Integrators",
2540 "[MixedScalarWeakDivergenceIntegrator]"
2541 "[MixedScalarVectorIntegrator]"
2542 "[BilinearFormIntegrator]"
2543 "[NonlinearFormIntegrator]")
2544 {
2545 int order = 2, n = 1, dim = 2;
2546 double tol = 1e-9;
2547
2548 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
2549
2550 VectorFunctionCoefficient V2_coef(dim, V2);
2551
2552 SECTION("Operators on H1")
2553 {
2554 H1_FECollection fec_h1(order, dim);
2555 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2556
2557 SECTION("Mapping H1 to H1")
2558 {
2559 SECTION("With Vector Coefficient")
2560 {
2561 MixedBilinearForm blf(&fespace_h1, &fespace_h1);
2562 blf.AddDomainIntegrator(
2563 new MixedDirectionalDerivativeIntegrator(V2_coef));
2564 blf.Assemble();
2565 blf.Finalize();
2566
2567 MixedBilinearForm blfw(&fespace_h1, &fespace_h1);
2568 blfw.AddDomainIntegrator(
2569 new MixedScalarWeakDivergenceIntegrator(V2_coef));
2570 blfw.Assemble();
2571 blfw.Finalize();
2572
2573 SparseMatrix * blfT = Transpose(blfw.SpMat());
2574 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2575
2576 REQUIRE( diff->MaxNorm() < tol );
2577
2578 delete blfT;
2579 delete diff;
2580 }
2581 }
2582 }
2583 for (int map_type = (int)FiniteElement::VALUE;
2584 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
2585 {
2586 SECTION("Operators on L2 (" +
2587 MapTypeName((FiniteElement::MapType)map_type) + ")")
2588 {
2589 L2_FECollection fec_l2(order - 1, dim,
2590 BasisType::GaussLegendre,
2591 (FiniteElement::MapType)map_type);
2592 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
2593
2594 SECTION("Mapping L2 (" +
2595 MapTypeName((FiniteElement::MapType)map_type) + ") to H1")
2596 {
2597 H1_FECollection fec_h1(order, dim);
2598 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2599
2600 SECTION("With Vector Coefficient")
2601 {
2602 MixedBilinearForm blf(&fespace_h1, &fespace_l2);
2603 blf.AddDomainIntegrator(
2604 new MixedDirectionalDerivativeIntegrator(V2_coef));
2605 blf.Assemble();
2606 blf.Finalize();
2607
2608 MixedBilinearForm blfw(&fespace_l2, &fespace_h1);
2609 blfw.AddDomainIntegrator(
2610 new MixedScalarWeakDivergenceIntegrator(V2_coef));
2611 blfw.Assemble();
2612 blfw.Finalize();
2613
2614 SparseMatrix * blfT = Transpose(blfw.SpMat());
2615 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2616
2617 REQUIRE( diff->MaxNorm() < tol );
2618
2619 delete blfT;
2620 delete diff;
2621 }
2622 }
2623 }
2624 }
2625 }
2626
2627
2628 TEST_CASE("2D Bilinear Vector Weak Divergence Integrators",
2629 "[MixedVectorWeakDivergenceIntegrator]"
2630 "[MixedVectorIntegrator]"
2631 "[BilinearFormIntegrator]"
2632 "[NonlinearFormIntegrator]")
2633 {
2634 int order = 2, n = 1, dim = 2;
2635 double tol = 1e-9;
2636
2637 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
2638
2639 FunctionCoefficient q2_coef(q2);
2640 VectorFunctionCoefficient D2_coef(dim, V2);
2641 MatrixFunctionCoefficient M2_coef(dim, M2);
2642 MatrixFunctionCoefficient MT2_coef(dim, MT2);
2643
2644 SECTION("Operators on ND")
2645 {
2646 ND_FECollection fec_nd(order, dim);
2647 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
2648
2649 SECTION("Mapping ND to H1")
2650 {
2651 H1_FECollection fec_h1(order, dim);
2652 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2653
2654 SECTION("Without Coefficient")
2655 {
2656 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
2657 blf.AddDomainIntegrator(
2658 new MixedVectorGradientIntegrator());
2659 blf.Assemble();
2660 blf.Finalize();
2661
2662 MixedBilinearForm blfw(&fespace_nd, &fespace_h1);
2663 blfw.AddDomainIntegrator(
2664 new MixedVectorWeakDivergenceIntegrator());
2665 blfw.Assemble();
2666 blfw.Finalize();
2667
2668 SparseMatrix * blfT = Transpose(blfw.SpMat());
2669 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2670
2671 REQUIRE( diff->MaxNorm() < tol );
2672
2673 delete blfT;
2674 delete diff;
2675 }
2676 SECTION("With Scalar Coefficient")
2677 {
2678 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
2679 blf.AddDomainIntegrator(
2680 new MixedVectorGradientIntegrator(q2_coef));
2681 blf.Assemble();
2682 blf.Finalize();
2683
2684 MixedBilinearForm blfw(&fespace_nd, &fespace_h1);
2685 blfw.AddDomainIntegrator(
2686 new MixedVectorWeakDivergenceIntegrator(q2_coef));
2687 blfw.Assemble();
2688 blfw.Finalize();
2689
2690 SparseMatrix * blfT = Transpose(blfw.SpMat());
2691 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2692
2693 REQUIRE( diff->MaxNorm() < tol );
2694
2695 delete blfT;
2696 delete diff;
2697 }
2698 SECTION("With Diagonal Matrix Coefficient")
2699 {
2700 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
2701 blf.AddDomainIntegrator(
2702 new MixedVectorGradientIntegrator(D2_coef));
2703 blf.Assemble();
2704 blf.Finalize();
2705
2706 MixedBilinearForm blfw(&fespace_nd, &fespace_h1);
2707 blfw.AddDomainIntegrator(
2708 new MixedVectorWeakDivergenceIntegrator(D2_coef));
2709 blfw.Assemble();
2710 blfw.Finalize();
2711
2712 SparseMatrix * blfT = Transpose(blfw.SpMat());
2713 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2714
2715 REQUIRE( diff->MaxNorm() < tol );
2716
2717 delete blfT;
2718 delete diff;
2719 }
2720 SECTION("With Matrix Coefficient")
2721 {
2722 MixedBilinearForm blf(&fespace_h1, &fespace_nd);
2723 blf.AddDomainIntegrator(
2724 new MixedVectorGradientIntegrator(MT2_coef));
2725 blf.Assemble();
2726 blf.Finalize();
2727
2728 MixedBilinearForm blfw(&fespace_nd, &fespace_h1);
2729 blfw.AddDomainIntegrator(
2730 new MixedVectorWeakDivergenceIntegrator(M2_coef));
2731 blfw.Assemble();
2732 blfw.Finalize();
2733
2734 SparseMatrix * blfT = Transpose(blfw.SpMat());
2735 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2736
2737 REQUIRE( diff->MaxNorm() < tol );
2738
2739 delete blfT;
2740 delete diff;
2741 }
2742 }
2743 }
2744 SECTION("Operators on RT")
2745 {
2746 RT_FECollection fec_rt(order - 1, dim);
2747 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
2748
2749 SECTION("Mapping RT to H1")
2750 {
2751 H1_FECollection fec_h1(order, dim);
2752 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2753
2754 SECTION("Without Coefficient")
2755 {
2756 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
2757 blf.AddDomainIntegrator(
2758 new MixedVectorGradientIntegrator());
2759 blf.Assemble();
2760 blf.Finalize();
2761
2762 MixedBilinearForm blfw(&fespace_rt, &fespace_h1);
2763 blfw.AddDomainIntegrator(
2764 new MixedVectorWeakDivergenceIntegrator());
2765 blfw.Assemble();
2766 blfw.Finalize();
2767
2768 SparseMatrix * blfT = Transpose(blfw.SpMat());
2769 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2770
2771 REQUIRE( diff->MaxNorm() < tol );
2772
2773 delete blfT;
2774 delete diff;
2775 }
2776 SECTION("With Scalar Coefficient")
2777 {
2778 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
2779 blf.AddDomainIntegrator(
2780 new MixedVectorGradientIntegrator(q2_coef));
2781 blf.Assemble();
2782 blf.Finalize();
2783
2784 MixedBilinearForm blfw(&fespace_rt, &fespace_h1);
2785 blfw.AddDomainIntegrator(
2786 new MixedVectorWeakDivergenceIntegrator(q2_coef));
2787 blfw.Assemble();
2788 blfw.Finalize();
2789
2790 SparseMatrix * blfT = Transpose(blfw.SpMat());
2791 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2792
2793 REQUIRE( diff->MaxNorm() < tol );
2794
2795 delete blfT;
2796 delete diff;
2797 }
2798 SECTION("With Diagonal Matrix Coefficient")
2799 {
2800 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
2801 blf.AddDomainIntegrator(
2802 new MixedVectorGradientIntegrator(D2_coef));
2803 blf.Assemble();
2804 blf.Finalize();
2805
2806 MixedBilinearForm blfw(&fespace_rt, &fespace_h1);
2807 blfw.AddDomainIntegrator(
2808 new MixedVectorWeakDivergenceIntegrator(D2_coef));
2809 blfw.Assemble();
2810 blfw.Finalize();
2811
2812 SparseMatrix * blfT = Transpose(blfw.SpMat());
2813 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2814
2815 REQUIRE( diff->MaxNorm() < tol );
2816
2817 delete blfT;
2818 delete diff;
2819 }
2820 SECTION("With Matrix Coefficient")
2821 {
2822 MixedBilinearForm blf(&fespace_h1, &fespace_rt);
2823 blf.AddDomainIntegrator(
2824 new MixedVectorGradientIntegrator(MT2_coef));
2825 blf.Assemble();
2826 blf.Finalize();
2827
2828 MixedBilinearForm blfw(&fespace_rt, &fespace_h1);
2829 blfw.AddDomainIntegrator(
2830 new MixedVectorWeakDivergenceIntegrator(M2_coef));
2831 blfw.Assemble();
2832 blfw.Finalize();
2833
2834 SparseMatrix * blfT = Transpose(blfw.SpMat());
2835 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
2836
2837 REQUIRE( diff->MaxNorm() < tol );
2838
2839 delete blfT;
2840 delete diff;
2841 }
2842 }
2843 }
2844 }
2845
2846
2847 TEST_CASE("2D Bilinear Dot Product Integrators",
2848 "[MixedDotProductIntegrator]"
2849 "[MixedScalarVectorIntegrator]"
2850 "[BilinearFormIntegrator]"
2851 "[NonlinearFormIntegrator]")
2852 {
2853 int order = 2, n = 1, dim = 2;
2854 double cg_rtol = 1e-14;
2855 double tol = 1e-9;
2856
2857 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
2858
2859 VectorFunctionCoefficient F2_coef(dim, F2);
2860 VectorFunctionCoefficient V2_coef(dim, V2);
2861 FunctionCoefficient VF2_coef(VdotF2);
2862
2863 SECTION("Operators on ND")
2864 {
2865 ND_FECollection fec_nd(order, dim);
2866 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
2867
2868 GridFunction f_nd(&fespace_nd); f_nd.ProjectCoefficient(F2_coef);
2869
2870 SECTION("Mapping ND to H1")
2871 {
2872 H1_FECollection fec_h1(order, dim);
2873 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2874
2875 BilinearForm m_h1(&fespace_h1);
2876 m_h1.AddDomainIntegrator(new MassIntegrator());
2877 m_h1.Assemble();
2878 m_h1.Finalize();
2879
2880 GridFunction g_h1(&fespace_h1);
2881
2882 Vector tmp_h1(fespace_h1.GetNDofs());
2883
2884 SECTION("With Vector Coefficient")
2885 {
2886 MixedBilinearForm blf(&fespace_nd, &fespace_h1);
2887 blf.AddDomainIntegrator(
2888 new MixedDotProductIntegrator(V2_coef));
2889 blf.Assemble();
2890 blf.Finalize();
2891
2892 blf.Mult(f_nd, tmp_h1); g_h1 = 0.0;
2893 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
2894
2895 REQUIRE( g_h1.ComputeL2Error(VF2_coef) < tol );
2896 }
2897 }
2898 for (int map_type = (int)FiniteElement::VALUE;
2899 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
2900 {
2901 SECTION("Mapping ND to L2 (" +
2902 MapTypeName((FiniteElement::MapType)map_type) + ")")
2903 {
2904 L2_FECollection fec_l2(order, dim,
2905 BasisType::GaussLegendre,
2906 (FiniteElement::MapType)map_type);
2907 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
2908
2909 BilinearForm m_l2(&fespace_l2);
2910 m_l2.AddDomainIntegrator(new MassIntegrator());
2911 m_l2.Assemble();
2912 m_l2.Finalize();
2913
2914 GridFunction g_l2(&fespace_l2);
2915
2916 Vector tmp_l2(fespace_l2.GetNDofs());
2917
2918 SECTION("With Vector Coefficient")
2919 {
2920 MixedBilinearForm blf(&fespace_nd, &fespace_l2);
2921 blf.AddDomainIntegrator(
2922 new MixedDotProductIntegrator(V2_coef));
2923 blf.Assemble();
2924 blf.Finalize();
2925
2926 blf.Mult(f_nd, tmp_l2); g_l2 = 0.0;
2927 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
2928
2929 REQUIRE( g_l2.ComputeL2Error(VF2_coef) < tol );
2930 }
2931 }
2932 }
2933 }
2934 SECTION("Operators on RT")
2935 {
2936 RT_FECollection fec_rt(order - 1, dim);
2937 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
2938
2939 GridFunction f_rt(&fespace_rt); f_rt.ProjectCoefficient(F2_coef);
2940
2941 SECTION("Mapping RT to H1")
2942 {
2943 H1_FECollection fec_h1(order, dim);
2944 FiniteElementSpace fespace_h1(&mesh, &fec_h1);
2945
2946 BilinearForm m_h1(&fespace_h1);
2947 m_h1.AddDomainIntegrator(new MassIntegrator());
2948 m_h1.Assemble();
2949 m_h1.Finalize();
2950
2951 GridFunction g_h1(&fespace_h1);
2952
2953 Vector tmp_h1(fespace_h1.GetNDofs());
2954
2955 SECTION("With Vector Coefficient")
2956 {
2957 MixedBilinearForm blf(&fespace_rt, &fespace_h1);
2958 blf.AddDomainIntegrator(
2959 new MixedDotProductIntegrator(V2_coef));
2960 blf.Assemble();
2961 blf.Finalize();
2962
2963 blf.Mult(f_rt, tmp_h1); g_h1 = 0.0;
2964 CG(m_h1, tmp_h1, g_h1, 0, 200, cg_rtol * cg_rtol, 0.0);
2965
2966 REQUIRE( g_h1.ComputeL2Error(VF2_coef) < tol );
2967 }
2968 }
2969 for (int map_type = (int)FiniteElement::VALUE;
2970 map_type <= (int)FiniteElement::INTEGRAL; map_type++)
2971 {
2972 SECTION("Mapping RT to L2 (" +
2973 MapTypeName((FiniteElement::MapType)map_type) + ")")
2974 {
2975 L2_FECollection fec_l2(order, dim,
2976 BasisType::GaussLegendre,
2977 (FiniteElement::MapType)map_type);
2978 FiniteElementSpace fespace_l2(&mesh, &fec_l2);
2979
2980 BilinearForm m_l2(&fespace_l2);
2981 m_l2.AddDomainIntegrator(new MassIntegrator());
2982 m_l2.Assemble();
2983 m_l2.Finalize();
2984
2985 GridFunction g_l2(&fespace_l2);
2986
2987 Vector tmp_l2(fespace_l2.GetNDofs());
2988
2989 SECTION("With Vector Coefficient")
2990 {
2991 MixedBilinearForm blf(&fespace_rt, &fespace_l2);
2992 blf.AddDomainIntegrator(
2993 new MixedDotProductIntegrator(V2_coef));
2994 blf.Assemble();
2995 blf.Finalize();
2996
2997 blf.Mult(f_rt, tmp_l2); g_l2 = 0.0;
2998 CG(m_l2, tmp_l2, g_l2, 0, 200, cg_rtol * cg_rtol, 0.0);
2999
3000 REQUIRE( g_l2.ComputeL2Error(VF2_coef) < tol );
3001 }
3002 }
3003 }
3004 }
3005 }
3006
3007 TEST_CASE("2D Bilinear Weak Gradient Dot Product Integrators",
3008 "[MixedWeakGradDotIntegrator]"
3009 "[MixedScalarVectorIntegrator]"
3010 "[BilinearFormIntegrator]"
3011 "[NonlinearFormIntegrator]")
3012 {
3013 int order = 2, n = 1, dim = 2;
3014 double tol = 1e-9;
3015
3016 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
3017
3018 VectorFunctionCoefficient V2_coef(dim, V2);
3019
3020 SECTION("Operators on ND")
3021 {
3022 ND_FECollection fec_nd(order, dim);
3023 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
3024
3025 SECTION("Mapping ND to RT")
3026 {
3027 RT_FECollection fec_rt(order - 1, dim);
3028 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
3029
3030 SECTION("With Vector Coefficient")
3031 {
3032 MixedBilinearForm blf(&fespace_rt, &fespace_nd);
3033 blf.AddDomainIntegrator(
3034 new MixedVectorDivergenceIntegrator(V2_coef));
3035 blf.Assemble();
3036 blf.Finalize();
3037
3038 MixedBilinearForm blfw(&fespace_nd, &fespace_rt);
3039 blfw.AddDomainIntegrator(
3040 new MixedWeakGradDotIntegrator(V2_coef));
3041 blfw.Assemble();
3042 blfw.Finalize();
3043
3044 SparseMatrix * blfT = Transpose(blfw.SpMat());
3045 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
3046
3047 REQUIRE( diff->MaxNorm() < tol );
3048
3049 delete blfT;
3050 delete diff;
3051 }
3052 }
3053 }
3054 SECTION("Operators on RT")
3055 {
3056 RT_FECollection fec_rt(order - 1, dim);
3057 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
3058
3059 SECTION("Mapping RT to RT")
3060 {
3061 SECTION("With Vector Coefficient")
3062 {
3063 MixedBilinearForm blf(&fespace_rt, &fespace_rt);
3064 blf.AddDomainIntegrator(
3065 new MixedVectorDivergenceIntegrator(V2_coef));
3066 blf.Assemble();
3067 blf.Finalize();
3068
3069 MixedBilinearForm blfw(&fespace_rt, &fespace_rt);
3070 blfw.AddDomainIntegrator(
3071 new MixedWeakGradDotIntegrator(V2_coef));
3072 blfw.Assemble();
3073 blfw.Finalize();
3074
3075 SparseMatrix * blfT = Transpose(blfw.SpMat());
3076 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
3077
3078 REQUIRE( diff->MaxNorm() < tol );
3079
3080 delete blfT;
3081 delete diff;
3082 }
3083 }
3084 }
3085 }
3086
3087 TEST_CASE("2D Bilinear Scalar Cross Product Curl Integrator",
3088 "[MixedScalarCrossCurlIntegrator]"
3089 "[MixedScalarVectorIntegrator]"
3090 "[BilinearFormIntegrator]"
3091 "[NonlinearFormIntegrator]")
3092 {
3093 int order = 2, n = 1, dim = 2;
3094 double cg_rtol = 1e-14;
3095 double tol = 1e-9;
3096
3097 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
3098
3099 VectorFunctionCoefficient F2_coef(dim, F2);
3100 VectorFunctionCoefficient V2_coef(dim, V2);
3101 VectorFunctionCoefficient VxdF2_coef(dim, VcrossCurlF2);
3102
3103 SECTION("Operators on ND")
3104 {
3105 ND_FECollection fec_nd(order, dim);
3106 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
3107
3108 GridFunction f_nd(&fespace_nd); f_nd.ProjectCoefficient(F2_coef);
3109
3110 SECTION("Mapping ND to RT")
3111 {
3112 RT_FECollection fec_rt(order - 1, dim);
3113 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
3114
3115 BilinearForm m_rt(&fespace_rt);
3116 m_rt.AddDomainIntegrator(new VectorFEMassIntegrator());
3117 m_rt.Assemble();
3118 m_rt.Finalize();
3119
3120 GridFunction g_rt(&fespace_rt);
3121
3122 Vector tmp_rt(fespace_rt.GetNDofs());
3123
3124 SECTION("With Vector Coefficient")
3125 {
3126 MixedBilinearForm blf(&fespace_nd, &fespace_rt);
3127 blf.AddDomainIntegrator(
3128 new MixedScalarCrossCurlIntegrator(V2_coef));
3129 blf.Assemble();
3130 blf.Finalize();
3131
3132 blf.Mult(f_nd, tmp_rt); g_rt = 0.0;
3133 CG(m_rt, tmp_rt, g_rt, 0, 200, cg_rtol * cg_rtol, 0.0);
3134
3135 REQUIRE( g_rt.ComputeL2Error(VxdF2_coef) < tol );
3136 }
3137 }
3138 SECTION("Mapping ND to ND")
3139 {
3140 BilinearForm m_nd(&fespace_nd);
3141 m_nd.AddDomainIntegrator(new VectorFEMassIntegrator());
3142 m_nd.Assemble();
3143 m_nd.Finalize();
3144
3145 GridFunction g_nd(&fespace_nd);
3146
3147 Vector tmp_nd(fespace_nd.GetNDofs());
3148
3149 SECTION("With Vector Coefficient")
3150 {
3151 MixedBilinearForm blf(&fespace_nd, &fespace_nd);
3152 blf.AddDomainIntegrator(
3153 new MixedScalarCrossCurlIntegrator(V2_coef));
3154 blf.Assemble();
3155 blf.Finalize();
3156
3157 blf.Mult(f_nd, tmp_nd); g_nd = 0.0;
3158 CG(m_nd, tmp_nd, g_nd, 0, 200, cg_rtol * cg_rtol, 0.0);
3159
3160 REQUIRE( g_nd.ComputeL2Error(VxdF2_coef) < tol );
3161 }
3162 }
3163 }
3164 }
3165
3166 TEST_CASE("2D Bilinear Scalar Weak Curl Cross Integrators",
3167 "[MixedScalarWeakCurlCrossIntegrator]"
3168 "[MixedScalarVectorIntegrator]"
3169 "[BilinearFormIntegrator]"
3170 "[NonlinearFormIntegrator]")
3171 {
3172 int order = 2, n = 1, dim = 2;
3173 double tol = 1e-9;
3174
3175 Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
3176
3177 VectorFunctionCoefficient V2_coef(dim, V2);
3178
3179 SECTION("Operators on ND")
3180 {
3181 ND_FECollection fec_nd(order, dim);
3182 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
3183
3184 SECTION("Mapping ND to ND")
3185 {
3186 SECTION("With Vector Coefficient")
3187 {
3188 MixedBilinearForm blf(&fespace_nd, &fespace_nd);
3189 blf.AddDomainIntegrator(
3190 new MixedScalarCrossCurlIntegrator(V2_coef));
3191 blf.Assemble();
3192 blf.Finalize();
3193
3194 MixedBilinearForm blfw(&fespace_nd, &fespace_nd);
3195 blfw.AddDomainIntegrator(
3196 new MixedScalarWeakCurlCrossIntegrator(V2_coef));
3197 blfw.Assemble();
3198 blfw.Finalize();
3199
3200 SparseMatrix * blfT = Transpose(blfw.SpMat());
3201 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
3202
3203 REQUIRE( diff->MaxNorm() < tol );
3204
3205 delete blfT;
3206 delete diff;
3207 }
3208 }
3209 }
3210 SECTION("Operators on RT")
3211 {
3212 RT_FECollection fec_rt(order - 1, dim);
3213 FiniteElementSpace fespace_rt(&mesh, &fec_rt);
3214
3215 SECTION("Mapping RT to ND")
3216 {
3217 ND_FECollection fec_nd(order, dim);
3218 FiniteElementSpace fespace_nd(&mesh, &fec_nd);
3219
3220 SECTION("With Vector Coefficient")
3221 {
3222 MixedBilinearForm blf(&fespace_nd, &fespace_rt);
3223 blf.AddDomainIntegrator(
3224 new MixedScalarCrossCurlIntegrator(V2_coef));
3225 blf.Assemble();
3226 blf.Finalize();
3227
3228 MixedBilinearForm blfw(&fespace_rt, &fespace_nd);
3229 blfw.AddDomainIntegrator(
3230 new MixedScalarWeakCurlCrossIntegrator(V2_coef));
3231 blfw.Assemble();
3232 blfw.Finalize();
3233
3234 SparseMatrix * blfT = Transpose(blfw.SpMat());
3235 SparseMatrix * diff = Add(1.0, blf.SpMat(), 1.0, *blfT);
3236
3237 REQUIRE( diff->MaxNorm() < tol );
3238
3239 delete blfT;
3240 delete diff;
3241 }
3242 }
3243 }
3244 }
3245
3246 } // namespace bilininteg_2d
3247