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 #include <iostream> 16 17 using namespace mfem; 18 19 TEST_CASE("Test order of boundary integrators", 20 "[BilinearForm]") 21 { 22 // Create a simple mesh 23 int dim = 2, nx = 2, ny = 2, order = 2; 24 Element::Type e_type = Element::QUADRILATERAL; 25 Mesh mesh = Mesh::MakeCartesian2D(nx, ny, e_type); 26 27 H1_FECollection fec(order, dim); 28 FiniteElementSpace fes(&mesh, &fec); 29 30 SECTION("Order of restricted boundary integrators") 31 { 32 ConstantCoefficient one(1.0); 33 ConstantCoefficient two(2.0); 34 ConstantCoefficient three(3.0); 35 ConstantCoefficient four(4.0); 36 37 Array<int> bdr1(4); bdr1 = 0; bdr1[0] = 1; 38 Array<int> bdr2(4); bdr2 = 0; bdr2[1] = 1; 39 Array<int> bdr3(4); bdr3 = 0; bdr3[2] = 1; 40 Array<int> bdr4(4); bdr4 = 0; bdr4[3] = 1; 41 42 BilinearForm a1234(&fes); 43 a1234.AddBoundaryIntegrator(new MassIntegrator(one), bdr1); 44 a1234.AddBoundaryIntegrator(new MassIntegrator(two), bdr2); 45 a1234.AddBoundaryIntegrator(new MassIntegrator(three), bdr3); 46 a1234.AddBoundaryIntegrator(new MassIntegrator(four), bdr4); 47 a1234.Assemble(0); 48 a1234.Finalize(0); 49 50 BilinearForm a4321(&fes); 51 a4321.AddBoundaryIntegrator(new MassIntegrator(four), bdr4); 52 a4321.AddBoundaryIntegrator(new MassIntegrator(three), bdr3); 53 a4321.AddBoundaryIntegrator(new MassIntegrator(two), bdr2); 54 a4321.AddBoundaryIntegrator(new MassIntegrator(one), bdr1); 55 a4321.Assemble(0); 56 a4321.Finalize(0); 57 58 const SparseMatrix &A1234 = a1234.SpMat(); 59 const SparseMatrix &A4321 = a4321.SpMat(); 60 61 SparseMatrix *D = Add(1.0, A1234, -1.0, A4321); 62 63 REQUIRE(D->MaxNorm() == MFEM_Approx(0.0)); 64 65 delete D; 66 } 67 } 68 69 70 TEST_CASE("FormLinearSystem/SolutionScope", 71 "[BilinearForm]" 72 "[CUDA]") 73 { 74 // Create a simple mesh and FE space 75 int dim = 2, nx = 2, ny = 2, order = 2; 76 Element::Type e_type = Element::QUADRILATERAL; 77 Mesh mesh = Mesh::MakeCartesian2D(nx, ny, e_type); 78 79 H1_FECollection fec(order, dim); 80 FiniteElementSpace fes(&mesh, &fec); 81 int bdr_dof; 82 83 // Solve a PDE on the conforming mesh and FE space defined above, storing the 84 // result in 'sol'. 85 auto SolvePDE = [&](AssemblyLevel al, GridFunction &sol) __anon1fb6eaaf0102(AssemblyLevel al, GridFunction &sol) 86 { 87 // Linear form: rhs 88 ConstantCoefficient f(1.0); 89 LinearForm b(&fes); 90 b.AddDomainIntegrator(new DomainLFIntegrator(f)); 91 b.Assemble(); 92 // Bilinear form: matrix 93 BilinearForm a(&fes); 94 a.AddDomainIntegrator(new DiffusionIntegrator); 95 a.SetAssemblyLevel(al); 96 a.Assemble(); 97 // Setup b.c. 98 Array<int> ess_tdof_list; 99 REQUIRE(mesh.bdr_attributes.Max() > 0); 100 Array<int> bdr_attr_is_ess(mesh.bdr_attributes.Max()); 101 bdr_attr_is_ess = 1; 102 fes.GetEssentialTrueDofs(bdr_attr_is_ess, ess_tdof_list); 103 REQUIRE(ess_tdof_list.Size() > 0); 104 // Setup (on host) solution initial guess satisfying the desired b.c. 105 ConstantCoefficient zero(0.0); 106 sol.ProjectCoefficient(zero); // performed on host 107 // Setup the linear system 108 Vector B, X; 109 OperatorPtr A; 110 const bool copy_interior = true; // interior(sol) --> interior(X) 111 a.FormLinearSystem(ess_tdof_list, sol, b, A, X, B, copy_interior); 112 // Solve the system 113 CGSolver cg; 114 cg.SetMaxIter(2000); 115 cg.SetRelTol(1e-8); 116 cg.SetAbsTol(0.0); 117 cg.SetPrintLevel(0); 118 cg.SetOperator(*A); 119 cg.Mult(B, X); 120 // Recover the solution 121 a.RecoverFEMSolution(X, b, sol); 122 // Initialize the bdr_dof to be checked 123 ess_tdof_list.HostRead(); 124 bdr_dof = AsConst(ess_tdof_list)[0]; // here, L-dof is the same T-dof 125 }; 126 127 // Legacy full assembly 128 { 129 GridFunction sol(&fes); 130 SolvePDE(AssemblyLevel::LEGACYFULL, sol); 131 // Make sure the solution is still accessible after 'X' is destoyed 132 sol.HostRead(); 133 REQUIRE(AsConst(sol)(bdr_dof) == 0.0); 134 } 135 136 // Partial assembly 137 { 138 GridFunction sol(&fes); 139 SolvePDE(AssemblyLevel::PARTIAL, sol); 140 // Make sure the solution is still accessible after 'X' is destoyed 141 sol.HostRead(); 142 REQUIRE(AsConst(sol)(bdr_dof) == 0.0); 143 } 144 } 145