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 assemblediagonalpa
18 {
19
20 int dimension;
21
coeffFunction(const Vector & x)22 double coeffFunction(const Vector& x)
23 {
24 if (dimension == 2)
25 {
26 return sin(8.0 * M_PI * x[0]) * cos(6.0 * M_PI * x[1]) + 2.0;
27 }
28 else
29 {
30 return sin(8.0 * M_PI * x[0]) * cos(6.0 * M_PI * x[1]) *
31 sin(4.0 * M_PI * x[2]) +
32 2.0;
33 }
34 }
35
vectorCoeffFunction(const Vector & x,Vector & f)36 void vectorCoeffFunction(const Vector & x, Vector & f)
37 {
38 f = 0.0;
39 if (dimension > 1)
40 {
41 f[0] = sin(M_PI * x[1]);
42 f[1] = sin(2.5 * M_PI * x[0]);
43 }
44 if (dimension == 3)
45 {
46 f[2] = sin(6.1 * M_PI * x[2]);
47 }
48 }
49
asymmetricMatrixCoeffFunction(const Vector & x,DenseMatrix & f)50 void asymmetricMatrixCoeffFunction(const Vector & x, DenseMatrix & f)
51 {
52 f = 0.0;
53 if (dimension == 2)
54 {
55 f(0,0) = 1.1 + sin(M_PI * x[1]); // 1,1
56 f(1,0) = cos(1.3 * M_PI * x[1]); // 2,1
57 f(0,1) = cos(2.5 * M_PI * x[0]); // 1,2
58 f(1,1) = 1.1 + sin(4.9 * M_PI * x[0]); // 2,2
59 }
60 else if (dimension == 3)
61 {
62 f(0,0) = 1.1 + sin(M_PI * x[1]); // 1,1
63 f(0,1) = cos(2.5 * M_PI * x[0]); // 1,2
64 f(0,2) = sin(4.9 * M_PI * x[2]); // 1,3
65 f(1,0) = cos(M_PI * x[0]); // 2,1
66 f(1,1) = 1.1 + sin(6.1 * M_PI * x[1]); // 2,2
67 f(1,2) = cos(6.1 * M_PI * x[2]); // 2,3
68 f(2,0) = sin(1.5 * M_PI * x[1]); // 3,1
69 f(2,1) = cos(2.9 * M_PI * x[0]); // 3,2
70 f(2,2) = 1.1 + sin(6.1 * M_PI * x[2]); // 3,3
71 }
72 }
73
fullSymmetricMatrixCoeffFunction(const Vector & x,DenseMatrix & f)74 void fullSymmetricMatrixCoeffFunction(const Vector & x, DenseMatrix & f)
75 {
76 f = 0.0;
77 if (dimension == 2)
78 {
79 f(0,0) = 1.1 + sin(M_PI * x[1]); // 1,1
80 f(0,1) = cos(2.5 * M_PI * x[0]); // 1,2
81 f(1,1) = 1.1 + sin(4.9 * M_PI * x[0]); // 2,2
82 f(1,0) = f(0,1);
83 }
84 else if (dimension == 3)
85 {
86 f(0,0) = sin(M_PI * x[1]); // 1,1
87 f(0,1) = cos(2.5 * M_PI * x[0]); // 1,2
88 f(0,2) = sin(4.9 * M_PI * x[2]); // 1,3
89 f(1,1) = sin(6.1 * M_PI * x[1]); // 2,2
90 f(1,2) = cos(6.1 * M_PI * x[2]); // 2,3
91 f(2,2) = sin(6.1 * M_PI * x[2]); // 3,3
92 f(1,0) = f(0,1);
93 f(2,0) = f(0,2);
94 f(2,1) = f(1,2);
95 }
96 }
97
symmetricMatrixCoeffFunction(const Vector & x,DenseSymmetricMatrix & f)98 void symmetricMatrixCoeffFunction(const Vector & x, DenseSymmetricMatrix & f)
99 {
100 f = 0.0;
101 if (dimension == 2)
102 {
103 f(0,0) = 1.1 + sin(M_PI * x[1]); // 1,1
104 f(0,1) = cos(2.5 * M_PI * x[0]); // 1,2
105 f(1,1) = 1.1 + sin(4.9 * M_PI * x[0]); // 2,2
106 }
107 else if (dimension == 3)
108 {
109 f(0,0) = sin(M_PI * x[1]); // 1,1
110 f(0,1) = cos(2.5 * M_PI * x[0]); // 1,2
111 f(0,2) = sin(4.9 * M_PI * x[2]); // 1,3
112 f(1,1) = sin(6.1 * M_PI * x[1]); // 2,2
113 f(1,2) = cos(6.1 * M_PI * x[2]); // 2,3
114 f(2,2) = sin(6.1 * M_PI * x[2]); // 3,3
115 }
116 }
117
118 TEST_CASE("massdiag")
119 {
120 for (dimension = 2; dimension < 4; ++dimension)
121 {
122 for (int ne = 1; ne < 3; ++ne)
123 {
124 std::cout << "Testing " << dimension << "D partial assembly mass diagonal: "
125 << std::pow(ne, dimension) << " elements." << std::endl;
126 for (int order = 1; order < 5; ++order)
127 {
128 Mesh mesh;
129 if (dimension == 2)
130 {
131 mesh = Mesh::MakeCartesian2D(
132 ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
133 }
134 else
135 {
136 mesh = Mesh::MakeCartesian3D(
137 ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
138 }
139 FiniteElementCollection *h1_fec = new H1_FECollection(order, dimension);
140 FiniteElementSpace h1_fespace(&mesh, h1_fec);
141 BilinearForm paform(&h1_fespace);
142 ConstantCoefficient one(1.0);
143 paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
144 paform.AddDomainIntegrator(new MassIntegrator(one));
145 paform.Assemble();
146 Vector pa_diag(h1_fespace.GetVSize());
147 paform.AssembleDiagonal(pa_diag);
148
149 BilinearForm faform(&h1_fespace);
150 faform.AddDomainIntegrator(new MassIntegrator(one));
151 faform.Assemble();
152 faform.Finalize();
153 Vector assembly_diag(h1_fespace.GetVSize());
154 faform.SpMat().GetDiag(assembly_diag);
155
156 assembly_diag -= pa_diag;
157 double error = assembly_diag.Norml2();
158 std::cout << " order: " << order << ", error norm: " << error << std::endl;
159 REQUIRE(assembly_diag.Norml2() < 1.e-12);
160
161 delete h1_fec;
162 }
163 }
164 }
165 }
166
167 TEST_CASE("diffusiondiag")
168 {
169 for (dimension = 2; dimension < 4; ++dimension)
170 {
171 for (int ne = 1; ne < 3; ++ne)
172 {
173 std::cout << "Testing " << dimension <<
174 "D partial assembly diffusion diagonal: "
175 << std::pow(ne, dimension) << " elements." << std::endl;
176 for (int order = 1; order < 5; ++order)
177 {
178 Mesh mesh;
179 if (dimension == 2)
180 {
181 mesh = Mesh::MakeCartesian2D(
182 ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
183 }
184 else
185 {
186 mesh = Mesh::MakeCartesian3D(
187 ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
188 }
189 FiniteElementCollection *h1_fec = new H1_FECollection(order, dimension);
190 FiniteElementSpace h1_fespace(&mesh, h1_fec);
191
192 for (int coeffType = 0; coeffType < 5; ++coeffType)
193 {
194 Coefficient* coeff = nullptr;
195 VectorCoefficient* vcoeff = nullptr;
196 MatrixCoefficient* mcoeff = nullptr;
197 SymmetricMatrixCoefficient* smcoeff = nullptr;
198 if (coeffType == 0)
199 {
200 coeff = new ConstantCoefficient(12.34);
201 }
202 else if (coeffType == 1)
203 {
204 coeff = new FunctionCoefficient(&coeffFunction);
205 }
206 else if (coeffType == 2)
207 {
208 vcoeff = new VectorFunctionCoefficient(dimension, &vectorCoeffFunction);
209 }
210 else if (coeffType == 3)
211 {
212 mcoeff = new MatrixFunctionCoefficient(dimension,
213 &fullSymmetricMatrixCoeffFunction);
214 smcoeff = new SymmetricMatrixFunctionCoefficient(dimension,
215 &symmetricMatrixCoeffFunction);
216 }
217 else if (coeffType == 4)
218 {
219 mcoeff = new MatrixFunctionCoefficient(dimension,
220 &asymmetricMatrixCoeffFunction);
221 }
222
223 BilinearForm paform(&h1_fespace);
224 paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
225 BilinearForm faform(&h1_fespace);
226
227 if (coeffType >= 4)
228 {
229 paform.AddDomainIntegrator(new DiffusionIntegrator(*mcoeff));
230 faform.AddDomainIntegrator(new DiffusionIntegrator(*mcoeff));
231 }
232 else if (coeffType == 3)
233 {
234 paform.AddDomainIntegrator(new DiffusionIntegrator(*smcoeff));
235 faform.AddDomainIntegrator(new DiffusionIntegrator(*mcoeff));
236 }
237 else if (coeffType == 2)
238 {
239 paform.AddDomainIntegrator(new DiffusionIntegrator(*vcoeff));
240 faform.AddDomainIntegrator(new DiffusionIntegrator(*vcoeff));
241 }
242 else
243 {
244 paform.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
245 faform.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
246 }
247
248 paform.Assemble();
249 Vector pa_diag(h1_fespace.GetVSize());
250 paform.AssembleDiagonal(pa_diag);
251
252 faform.Assemble();
253 faform.Finalize();
254 Vector assembly_diag(h1_fespace.GetVSize());
255 faform.SpMat().GetDiag(assembly_diag);
256
257 assembly_diag -= pa_diag;
258 double error = assembly_diag.Norml2();
259 std::cout << " order: " << order << ", coefficient type "
260 << coeffType << ", error norm: " << error << std::endl;
261 REQUIRE(assembly_diag.Norml2() < 1.e-12);
262
263 delete coeff;
264 delete vcoeff;
265 delete mcoeff;
266 delete smcoeff;
267 }
268
269 delete h1_fec;
270 }
271 }
272 }
273 }
274
275 template <typename INTEGRATOR>
test_vdiagpa(int dim,int order)276 double test_vdiagpa(int dim, int order)
277 {
278 Mesh mesh;
279 if (dim == 2)
280 {
281 mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, 0, 1.0, 1.0);
282 }
283 else if (dim == 3)
284 {
285 mesh = Mesh::MakeCartesian3D(2, 2, 2, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
286 }
287
288 H1_FECollection fec(order, dim);
289 FiniteElementSpace fes(&mesh, &fec, dim);
290
291 BilinearForm form(&fes);
292 form.SetAssemblyLevel(AssemblyLevel::PARTIAL);
293 form.AddDomainIntegrator(new INTEGRATOR);
294 form.Assemble();
295
296 Vector diag(fes.GetVSize());
297 form.AssembleDiagonal(diag);
298
299 BilinearForm form_full(&fes);
300 form_full.AddDomainIntegrator(new INTEGRATOR);
301 form_full.Assemble();
302 form_full.Finalize();
303
304 Vector diag_full(fes.GetVSize());
305 form_full.SpMat().GetDiag(diag_full);
306
307 diag_full -= diag;
308
309 return diag_full.Norml2();
310 }
311
312 TEST_CASE("Vector Mass Diagonal PA", "[PartialAssembly], [AssembleDiagonal]")
313 {
314 SECTION("2D")
315 {
316 REQUIRE(test_vdiagpa<VectorMassIntegrator>(2,
317 2) == MFEM_Approx(0.0));
318
319 REQUIRE(test_vdiagpa<VectorMassIntegrator>(2,
320 3) == MFEM_Approx(0.0));
321 }
322
323 SECTION("3D")
324 {
325 REQUIRE(test_vdiagpa<VectorMassIntegrator>(3,
326 2) == MFEM_Approx(0.0));
327
328 REQUIRE(test_vdiagpa<VectorMassIntegrator>(3,
329 3) == MFEM_Approx(0.0));
330 }
331 }
332
333 TEST_CASE("Vector Diffusion Diagonal PA",
334 "[PartialAssembly], [AssembleDiagonal]")
335 {
336 SECTION("2D")
337 {
338 REQUIRE(
339 test_vdiagpa<VectorDiffusionIntegrator>(2,
340 2) == MFEM_Approx(0.0));
341
342 REQUIRE(test_vdiagpa<VectorDiffusionIntegrator>(2,
343 3) == MFEM_Approx(0.0));
344 }
345
346 SECTION("3D")
347 {
348 REQUIRE(test_vdiagpa<VectorDiffusionIntegrator>(3,
349 2) == MFEM_Approx(0.0));
350
351 REQUIRE(test_vdiagpa<VectorDiffusionIntegrator>(3,
352 3) == MFEM_Approx(0.0));
353 }
354 }
355
356 TEST_CASE("Hcurl/Hdiv diagonal PA",
357 "[CUDA]")
358 {
359 for (dimension = 2; dimension < 4; ++dimension)
360 {
361 for (int coeffType = 0; coeffType < 5; ++coeffType)
362 {
363 const int numSpaces = (coeffType == 0) ? 2 : 1;
364
365 Coefficient* coeff = nullptr;
366 DiagonalMatrixCoefficient* dcoeff = nullptr;
367 MatrixCoefficient* mcoeff = nullptr;
368 SymmetricMatrixCoefficient* smcoeff = nullptr;
369 if (coeffType == 0)
370 {
371 coeff = new ConstantCoefficient(12.34);
372 }
373 else if (coeffType == 1)
374 {
375 coeff = new FunctionCoefficient(&coeffFunction);
376 }
377 else if (coeffType == 2)
378 {
379 dcoeff = new VectorFunctionCoefficient(dimension, &vectorCoeffFunction);
380 }
381 else if (coeffType == 3)
382 {
383 mcoeff = new MatrixFunctionCoefficient(dimension,
384 &fullSymmetricMatrixCoeffFunction);
385 smcoeff = new SymmetricMatrixFunctionCoefficient(dimension,
386 &symmetricMatrixCoeffFunction);
387 }
388 else if (coeffType == 4)
389 {
390 mcoeff = new MatrixFunctionCoefficient(dimension,
391 &asymmetricMatrixCoeffFunction);
392 }
393
394 enum Spaces {Hcurl, Hdiv};
395
396 for (int spaceType = 0; spaceType < numSpaces; ++spaceType)
397 {
398 const int numIntegrators = (dimension == 3 || coeffType < 2) ? 2 : 1;
399 for (int integrator = 0; integrator < numIntegrators; ++integrator)
400 {
401 for (int ne = 1; ne < 3; ++ne)
402 {
403 if (spaceType == Hcurl)
404 std::cout << "Testing " << dimension <<
405 "D partial assembly H(curl) diagonal for integrator " << integrator
406 << " and coeffType " << coeffType << ": "
407 << std::pow(ne, dimension) << " elements." << std::endl;
408 else
409 std::cout << "Testing " << dimension <<
410 "D partial assembly H(div) diagonal for integrator " << integrator
411 << " and coeffType " << coeffType << ": "
412 << std::pow(ne, dimension) << " elements." << std::endl;
413
414 for (int order = 1; order < 4; ++order)
415 {
416 Mesh mesh;
417 if (dimension == 2)
418 {
419 mesh = Mesh::MakeCartesian2D(
420 ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
421 }
422 else
423 {
424 mesh = Mesh::MakeCartesian3D(
425 ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
426 }
427
428 FiniteElementCollection* fec = (spaceType == Hcurl) ?
429 (FiniteElementCollection*) new ND_FECollection(order, dimension) :
430 (FiniteElementCollection*) new RT_FECollection(order, dimension);
431
432 FiniteElementSpace fespace(&mesh, fec);
433 BilinearForm paform(&fespace);
434 BilinearForm faform(&fespace);
435 paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
436 if (integrator == 0)
437 {
438 if (coeffType >= 4)
439 {
440 paform.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
441 faform.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
442 }
443 else if (coeffType == 3)
444 {
445 paform.AddDomainIntegrator(new VectorFEMassIntegrator(*smcoeff));
446 faform.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
447 }
448 else if (coeffType == 2)
449 {
450 paform.AddDomainIntegrator(new VectorFEMassIntegrator(*dcoeff));
451 faform.AddDomainIntegrator(new VectorFEMassIntegrator(*dcoeff));
452 }
453 else
454 {
455 paform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
456 faform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
457 }
458 }
459 else
460 {
461 if (spaceType == Hcurl)
462 {
463 const FiniteElement *fel = fespace.GetFE(0);
464 const IntegrationRule *intRule = &MassIntegrator::GetRule(*fel, *fel,
465 *mesh.GetElementTransformation(0));
466
467 if (coeffType >= 4)
468 {
469 paform.AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, intRule));
470 faform.AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, intRule));
471 }
472 else if (coeffType == 3)
473 {
474 paform.AddDomainIntegrator(new CurlCurlIntegrator(*smcoeff, intRule));
475 faform.AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, intRule));
476 }
477 else if (coeffType == 2)
478 {
479 paform.AddDomainIntegrator(new CurlCurlIntegrator(*dcoeff, intRule));
480 faform.AddDomainIntegrator(new CurlCurlIntegrator(*dcoeff, intRule));
481 }
482 else
483 {
484 paform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff, intRule));
485 faform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff, intRule));
486 }
487 }
488 else
489 {
490 paform.AddDomainIntegrator(new DivDivIntegrator(*coeff));
491 faform.AddDomainIntegrator(new DivDivIntegrator(*coeff));
492 }
493 }
494 paform.Assemble();
495 Vector pa_diag(fespace.GetVSize());
496 paform.AssembleDiagonal(pa_diag);
497
498 faform.Assemble();
499 faform.Finalize();
500 Vector assembly_diag(fespace.GetVSize());
501 faform.SpMat().GetDiag(assembly_diag);
502
503 assembly_diag -= pa_diag;
504 double error = assembly_diag.Norml2();
505 std::cout << " order: " << order << ", error norm: " << error << std::endl;
506 REQUIRE(assembly_diag.Norml2() < 1.e-11);
507
508 delete fec;
509 }
510 } // ne
511 } // integrator
512 } // spaceType
513
514 delete coeff;
515 delete dcoeff;
516 delete mcoeff;
517 delete smcoeff;
518 } // coeffType
519 } // dimension
520 }
521
522 } // namespace assemblediagonalpa
523